<div dir="ltr"><div>Great thanks! I can now successfully run <a href="https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/ksp/tests/ex71f.F90">https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/ksp/tests/ex71f.F90</a>.</div><div><br>Going forward with my experiments, let me post a new code snippet (very similar to
ex71f.F90) that I cannot get to work, probably I must be setting up the IS objects incorrectly.</div><div><br></div><div>I have an 8x8 matrix A=diag(1,1,2,2,...,2) and a vector b=(0.5,...,0.5). We have only one processor, and I want to solve Ax=b using GASM. In particular, KSP is set to preonly, GASM is the preconditioner and it uses on each submatrix an lu direct solver (sub_ksp = preonly, sub_pc = lu).</div><div><br></div><div>For the GASM algorithm, I divide A into diag(1,1) and diag(2,2,...,2). For simplicity I set 0 overlap. Now I want to use GASM to solve Ax=b. The code follows.</div><div><font face="monospace"><br>#include <petsc/finclude/petscmat.h><br>#include <petsc/finclude/petscksp.h><br>#include <petsc/finclude/petscpc.h><br> USE petscmat<br> USE petscksp<br> USE petscpc<br> USE MPI<br><br></font></div><div><font face="monospace"> Mat :: A<br> Vec :: b, x<br> PetscInt :: M, I, J, ISLen, NSub<br> PetscMPIInt :: size<br> PetscErrorCode :: ierr<br> PetscScalar :: v<br> KSP :: ksp<br> PC :: pc<br> IS :: subdomains_IS(2), inflated_IS(2)<br> PetscInt,DIMENSION(4) :: indices_first_domain<br> PetscInt,DIMENSION(36) :: indices_second_domain<br> <br> call PetscInitialize(PETSC_NULL_CHARACTER, ierr)<br> call MPI_Comm_size(PETSC_COMM_WORLD, size, ierr)<br><br><br> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!<br> ! INTRO: create matrix and right hand side<br> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!<br> <br> WRITE(*,*) "Assembling A,b"<br><br> M = 8<br> call MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE,<br> & M, M, PETSC_DEFAULT_INTEGER, PETSC_NULL_INTEGER,<br> & PETSC_DEFAULT_INTEGER, PETSC_NULL_INTEGER,A, ierr)<br> DO I=1,M<br> DO J=1,M<br> IF ((I .EQ. J) .AND. (I .LE. 2 )) THEN<br> v = 1<br> ELSE IF ((I .EQ. J) .AND. (I .GT. 2 )) THEN<br> v = 2<br> ELSE<br> v = 0<br> ENDIF<br> call MatSetValue(A, I-1, J-1, v, INSERT_VALUES, ierr)<br> END DO<br> END DO<br><br> call MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr)<br> call MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr)<br> <br> call VecCreate(PETSC_COMM_WORLD,b,ierr)<br> call VecSetSizes(b, PETSC_DECIDE, M,ierr)<br> call VecSetFromOptions(b,ierr)<br> <br> do I=1,M<br> v = 0.5<br> call VecSetValue(b,I-1,v, INSERT_VALUES,ierr)<br> end do<br> <br> call VecAssemblyBegin(b,ierr)<br> call VecAssemblyEnd(b,ierr)<br> </font></div><div><font face="monospace"> <br> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!<br> ! FIRST KSP/PC SETUP<br> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!<br> <br> WRITE(*,*) "KSP/PC first setup"<br> <br> call KSPCreate(PETSC_COMM_WORLD, ksp, ierr)<br> call KSPSetOperators(ksp, A, A, ierr)<br> call KSPSetType(ksp, 'preonly', ierr)<br> call KSPGetPC(ksp, pc, ierr)<br> call KSPSetUp(ksp, ierr)<br> call PCSetType(pc, PCGASM, ierr)<br> </font></div><div><font face="monospace"> <br> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!<br> ! GASM, SETTING SUBDOMAINS<br> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!<br> <br> WRITE(*,*) "Setting GASM subdomains"<br> <br> ! Let's create the subdomain IS and inflated_IS<br> ! They are equal if no overlap is present<br> ! They are 1: 0,1,8,9<br> ! 2: 10,...,15,18,...,23,...,58,...,63 <br> <br> indices_first_domain = [0,1,8,9] ! corresponds to diag(1,1)<br> do I=0,5<br> do J=0,5<br> indices_second_domain(I*6+1+J) = 18 + J + 8*I ! corresponds to diag(2,2,...,2)<br> !WRITE(*,*) I*6+1+J, 18 + J + 8*I<br> end do<br> end do<br> <br> ! Convert into IS<br> ISLen = 4<br> call ISCreateGeneral(PETSC_COMM_WORLD,ISLen,indices_first_domain,<br> & PETSC_COPY_VALUES, subdomains_IS(1), ierr)<br> call ISCreateGeneral(PETSC_COMM_WORLD,ISLen,indices_first_domain,<br> & PETSC_COPY_VALUES, inflated_IS(1), ierr)<br> ISLen = 36<br> call ISCreateGeneral(PETSC_COMM_WORLD,ISLen,indices_second_domain,<br> & PETSC_COPY_VALUES, subdomains_IS(2), ierr)<br> call ISCreateGeneral(PETSC_COMM_WORLD,ISLen,indices_second_domain,<br> & PETSC_COPY_VALUES, inflated_IS(2), ierr)<br><br> NSub = 2<br> call PCGASMSetSubdomains(pc,NSub,<br> & subdomains_IS,inflated_IS,ierr)<br> call PCGASMDestroySubdomains(NSub,<br> & subdomains_IS,inflated_IS,ierr)<br> <br> <br> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!<br> ! GASM: SET SUBSOLVERS<br> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! <br> <br> WRITE(*,*) "Setting subsolvers for GASM"<br> <br> call PCSetUp(pc, ierr) ! should I add this?<br> <br> call PetscOptionsSetValue(PETSC_NULL_OPTIONS,<br> & "-sub_pc_type", "lu", ierr)<br> call PetscOptionsSetValue(PETSC_NULL_OPTIONS,<br> & "-sub_ksp_type", "preonly", ierr)<br> <br> call KSPSetFromOptions(ksp, ierr)<br> call PCSetFromOptions(pc, ierr)<br> </font></div><div><font face="monospace"> <br> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!<br> ! DUMMY SOLUTION: DID IT WORK?<br> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! <br> <br> WRITE(*,*) "Solve"<br> <br> call VecDuplicate(b,x,ierr)<br> call KSPSolve(ksp,b,x,ierr)<br> <br> call MatDestroy(A, ierr)<br> call KSPDestroy(ksp, ierr)<br> call PetscFinalize(ierr)</font><br></div><div><br></div><div>This code is failing in multiple points. At <font face="monospace">call PCSetUp(pc, ierr)</font> it produces:</div><div><br></div><blockquote style="margin:0 0 0 40px;border:none;padding:0px"><div><i><font color="#ff0000">[0]PETSC ERROR: Argument out of range</font></i></div><div><i><font color="#ff0000">[0]PETSC ERROR: Scatter indices in ix are out of range</font></i></div><div><i><font color="#ff0000">...</font></i></div><div><i><font color="#ff0000">[0]PETSC ERROR: #1 VecScatterCreate() at ***\src\vec\is\sf\INTERF~1\vscat.c:736</font></i></div><div><i><font color="#ff0000">[0]PETSC ERROR: #2 PCSetUp_GASM() at ***\src\ksp\pc\impls\gasm\gasm.c:433</font></i></div><div><i><font color="#ff0000">[0]PETSC ERROR: #3 PCSetUp() at ***\src\ksp\pc\INTERF~1\precon.c:994</font></i></div><div><i><br></i></div></blockquote><div>And at <span style="font-family:monospace">call KSPSolve(ksp,b,x,ierr) </span><font face="arial, sans-serif">it produces:</font></div><div><font face="arial, sans-serif"><br></font></div><blockquote style="margin:0 0 0 40px;border:none;padding:0px"><div><i><font color="#ff0000">forrtl: severe (157): Program Exception - access violation</font></i></div></blockquote><div><br></div><div>The index sets are setup coherently with the outputs of e.g. <a href="https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/ksp/tests/output/ex71f_1.out">https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/ksp/tests/output/ex71f_1.out</a>: in particular each element of the matrix A corresponds to a number from 0 to 63.</div><div>Note that each submatrix does not represent some physical subdomain, the subdivision is just at the algebraic level.</div><div>I thus have the following questions:</div><div><ul><li>is this the correct way of creating the IS objects, given my objective at the beginning of the email? Is the ordering correct?</li><li>what am I doing wrong that is generating the above errors?</li></ul><div>Thanks for the patience and the time.</div><div>Best, <br>Leonardo</div></div><div class="gmail_quote"><div dir="ltr" class="gmail_attr"><br></div><div dir="ltr" class="gmail_attr">Il giorno ven 5 mag 2023 alle ore 18:43 Barry Smith <<a href="mailto:bsmith@petsc.dev">bsmith@petsc.dev</a>> ha scritto:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div><div><br></div> Added in <b style="color:rgb(200,20,201);font-family:Menlo;font-size:14px">barry/2023-05-04/add-pcgasm-set-subdomains </b>see also <a href="https://gitlab.com/petsc/petsc/-/merge_requests/6419" target="_blank">https://gitlab.com/petsc/petsc/-/merge_requests/6419</a><div><br></div><div> Barry</div><div><br><div><br><blockquote type="cite"><div>On May 4, 2023, at 11:23 AM, LEONARDO MUTTI <<a href="mailto:leonardo.mutti01@universitadipavia.it" target="_blank">leonardo.mutti01@universitadipavia.it</a>> wrote:</div><br><div><div dir="ltr">Thank you for the help. <div>Adding to my example:<div><i> call PCGASMSetSubdomains(pc,NSub, subdomains_IS, inflated_IS,ierr)<br> call PCGASMDestroySubdomains(NSub,subdomains_IS,inflated_IS,ierr)<br></i></div><div>results in:<br></div><div><div><i> Error LNK2019 unresolved external symbol PCGASMDESTROYSUBDOMAINS referenced in function ... <br></i></div><div><i> Error LNK2019 unresolved external symbol PCGASMSETSUBDOMAINS referenced in function ... <br></i></div><div>I'm not sure if the interfaces are missing or if I have a compilation problem.</div></div><div>Thank you again. </div><div>Best,<br>Leonardo</div></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">Il giorno sab 29 apr 2023 alle ore 20:30 Barry Smith <<a href="mailto:bsmith@petsc.dev" target="_blank">bsmith@petsc.dev</a>> ha scritto:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div><div><br></div> Thank you for the test code. I have a fix in the branch <span style="color:rgb(51,50,56);font-family:"GitLab Sans",-apple-system,BlinkMacSystemFont,"Segoe UI",Roboto,"Noto Sans",Ubuntu,Cantarell,"Helvetica Neue",sans-serif,"Apple Color Emoji","Segoe UI Emoji","Segoe UI Symbol","Noto Color Emoji";font-size:14px;background-color:rgb(251,250,253)"> </span><a href="https://gitlab.com/petsc/petsc/-/commits/barry/2023-04-29/fix-pcasmcreatesubdomains2d" style="box-sizing:border-box;font-variant-ligatures:none;color:rgb(31,117,203);text-decoration:none;font-size:13.3px" target="_blank">barry/2023-04-29/fix-pcasmcreatesubdomains2d</a> with merge request <a href="https://gitlab.com/petsc/petsc/-/merge_requests/6394" target="_blank">https://gitlab.com/petsc/petsc/-/merge_requests/6394</a><div><br></div><div> The functions did not have proper Fortran stubs and interfaces so I had to provide them manually in the new branch.<br><div><br></div><div> Use</div><div><br></div><div> git fetch</div><div> git checkout <a href="https://gitlab.com/petsc/petsc/-/commits/barry/2023-04-29/fix-pcasmcreatesubdomains2d" style="box-sizing:border-box;font-variant-ligatures:none;color:rgb(31,117,203);text-decoration:none;font-size:13.3px" target="_blank">barry/2023-04-29/fix-pcasmcreatesubdomains2d</a></div><div> ./configure etc</div><div><br></div><div> Your now working test code is in src/ksp/ksp/tests/ex71f.F90 I had to change things slightly and I updated the error handling for the latest version.</div><div><br></div><div> Please let us know if you have any later questions.</div><div><br></div><div> Barry</div><div><br></div><div><br></div><div><br></div><div><div><br><blockquote type="cite"><div>On Apr 28, 2023, at 12:07 PM, LEONARDO MUTTI <<a href="mailto:leonardo.mutti01@universitadipavia.it" target="_blank">leonardo.mutti01@universitadipavia.it</a>> wrote:</div><br><div><div dir="ltr">Hello. I am having a hard time understanding the index sets to feed PCGASMSetSubdomains, and I am working in Fortran (as a PETSc novice). To get more intuition on how the IS objects behave I tried the following minimal (non) working example, which should tile a 16x16 matrix into 16 square, non-overlapping submatrices:<div><br></div><div>#include <petsc/finclude/petscmat.h><br></div><div>#include <petsc/finclude/petscksp.h><br>#include <petsc/finclude/petscpc.h><br> USE petscmat<br> USE petscksp<br> USE petscpc<br> <br> Mat :: A<br> PetscInt :: M, NSubx, dof, overlap, NSub<br> INTEGER :: I,J<br> PetscErrorCode :: ierr<br> PetscScalar :: v<br> KSP :: ksp<br> PC :: pc<br> IS :: subdomains_IS, inflated_IS<br> <br> call PetscInitialize(PETSC_NULL_CHARACTER , ierr)</div><div> <br>!-----Create a dummy matrix<br> M = 16<br> call MatCreateAIJ(MPI_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE,<br> & M, M,<br> & PETSC_DEFAULT_INTEGER, PETSC_NULL_INTEGER,<br> & PETSC_DEFAULT_INTEGER, PETSC_NULL_INTEGER,<br> & A, ierr)<br> <br> DO I=1,M<br> DO J=1,M<br> v = I*J<br> CALL MatSetValue (A,I-1,J-1,v,<br> & INSERT_VALUES , ierr)<br> END DO<br> END DO<br> <br> call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY , ierr)<br> call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY , ierr)<br> </div><div>!-----Create KSP and PC<br> call KSPCreate(PETSC_COMM_WORLD,ksp, ierr)<br> call KSPSetOperators(ksp,A,A, ierr)<br> call KSPSetType(ksp,"bcgs",ierr)<br> call KSPGetPC(ksp,pc,ierr)<br> call KSPSetUp(ksp, ierr)<br> call PCSetType(pc,PCGASM, ierr)<br> call PCSetUp(pc , ierr)<br> <br>!-----GASM setup<br> NSubx = 4<br> dof = 1<br> overlap = 0<br> <br> call PCGASMCreateSubdomains2D(pc, <br> & M, M,<br> & NSubx, NSubx,<br> & dof, overlap,<br> & NSub, subdomains_IS, inflated_IS, ierr)<br><br> call ISView(subdomains_IS, PETSC_VIEWER_STDOUT_WORLD, ierr)<br> <br> call KSPDestroy(ksp, ierr) <br> call PetscFinalize(ierr)</div><div><br></div><div>Running this on one processor, I get NSub = 4.</div><div>If PCASM and PCASMCreateSubdomains2D are used instead, I get NSub = 16 as expected.</div><div>Moreover, I get in the end "forrtl: severe (157): Program Exception - access violation". So:</div><div>1) why do I get two different results with ASM, and GASM?</div><div>2) why do I get access violation and how can I solve this?</div><div>In fact, in C, subdomains_IS, inflated_IS should pointers to IS objects. As I see on the Fortran interface, the arguments to PCGASMCreateSubdomains2D are IS objects:<br></div><div><br></div><div> subroutine PCGASMCreateSubdomains2D(a,b,c,d,e,f,g,h,i,j,z)<br> import tPC,tIS<br> PC a ! PC<br> PetscInt b ! PetscInt<br> PetscInt c ! PetscInt<br> PetscInt d ! PetscInt<br> PetscInt e ! PetscInt<br> PetscInt f ! PetscInt<br> PetscInt g ! PetscInt<br> PetscInt h ! PetscInt<br> IS i ! IS<br> IS j ! IS<br> PetscErrorCode z<br> end subroutine PCGASMCreateSubdomains2D<br></div><div>Thus:</div><div>3) what should be inside e.g., subdomains_IS? I expect it to contain, for every created subdomain, the list of rows and columns defining the subblock in the matrix, am I right?</div><div><br></div><div>Context: I have a block-tridiagonal system arising from space-time finite elements, and I want to solve it with GMRES+PCGASM preconditioner, where each overlapping submatrix is on the diagonal and of size 3x3 blocks (and spanning multiple processes). This is PETSc 3.17.1 on Windows.</div><div><br></div><div>Thanks in advance,</div><div>Leonardo</div></div>
</div></blockquote></div><br></div></div></div></blockquote></div>
</div></blockquote></div><br></div></div></blockquote></div></div>