Partial typo: I expect 9x(16+16) numbers to be stored in subdomain_IS : # subdomains x (row indices of the submatrix + col indices of the submatrix).
Il giorno mar 9 mag 2023 alle ore 18:31 LEONARDO MUTTI < leonardo.mutt...@universitadipavia.it> ha scritto: > > > ---------- Forwarded message --------- > Da: LEONARDO MUTTI <leonardo.mutt...@universitadipavia.it> > Date: mar 9 mag 2023 alle ore 18:29 > Subject: Re: [petsc-users] Understanding index sets for PCGASM > To: Matthew Knepley <knep...@gmail.com> > > > Thank you for your answer, but I am still confused, sorry. > Consider > https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/ksp/tests/ex71f.F90 on > one processor. > Let M=12 for the sake of simplicity, i.e. we deal with a 12x12 2D grid, > hence, a 144x144 matrix. > Let NSubx = 3, so that on the grid we do 3 vertical and 3 horizontal > subdivisions. > We should obtain 9 subdomains that are grids of 4x4 nodes each, thus > corresponding to 9 submatrices of size 16x16. > In my run I obtain NSub = 9 (great) and subdomain_IS(i), i=1,...,9, reads: > > *IS Object: 1 MPI process* > * type: general* > *Number of indices in set 16* > *0 0* > *1 1* > *2 2* > *3 3* > *4 12* > *5 13* > *6 14* > *7 15* > *8 24* > *9 25* > *10 26* > *11 27* > *12 36* > *13 37* > *14 38* > *15 39* > *IS Object: 1 MPI process* > * type: general* > *Number of indices in set 16* > *0 4* > *1 5* > *2 6* > *3 7* > *4 16* > *5 17* > *6 18* > *7 19* > *8 28* > *9 29* > *10 30* > *11 31* > *12 40* > *13 41* > *14 42* > *15 43* > *IS Object: 1 MPI process* > * type: general* > *Number of indices in set 16* > *0 8* > *1 9* > *2 10* > *3 11* > *4 20* > *5 21* > *6 22* > *7 23* > *8 32* > *9 33* > *10 34* > *11 35* > *12 44* > *13 45* > *14 46* > *15 47* > *IS Object: 1 MPI process* > * type: general* > *Number of indices in set 16* > *0 48* > *1 49* > *2 50* > *3 51* > *4 60* > *5 61* > *6 62* > *7 63* > *8 72* > *9 73* > *10 74* > *11 75* > *12 84* > *13 85* > *14 86* > *15 87* > *IS Object: 1 MPI process* > * type: general* > *Number of indices in set 16* > *0 52* > *1 53* > *2 54* > *3 55* > *4 64* > *5 65* > *6 66* > *7 67* > *8 76* > *9 77* > *10 78* > *11 79* > *12 88* > *13 89* > *14 90* > *15 91* > *IS Object: 1 MPI process* > * type: general* > *Number of indices in set 16* > *0 56* > *1 57* > *2 58* > *3 59* > *4 68* > *5 69* > *6 70* > *7 71* > *8 80* > *9 81* > *10 82* > *11 83* > *12 92* > *13 93* > *14 94* > *15 95* > *IS Object: 1 MPI process* > * type: general* > *Number of indices in set 16* > *0 96* > *1 97* > *2 98* > *3 99* > *4 108* > *5 109* > *6 110* > *7 111* > *8 120* > *9 121* > *10 122* > *11 123* > *12 132* > *13 133* > *14 134* > *15 135* > *IS Object: 1 MPI process* > * type: general* > *Number of indices in set 16* > *0 100* > *1 101* > *2 102* > *3 103* > *4 112* > *5 113* > *6 114* > *7 115* > *8 124* > *9 125* > *10 126* > *11 127* > *12 136* > *13 137* > *14 138* > *15 139* > *IS Object: 1 MPI process* > * type: general* > *Number of indices in set 16* > *0 104* > *1 105* > *2 106* > *3 107* > *4 116* > *5 117* > *6 118* > *7 119* > *8 128* > *9 129* > *10 130* > *11 131* > *12 140* > *13 141* > *14 142* > *15 143* > > As you said, no number here reaches 144. > But the number stored in subdomain_IS are 9x16= #subdomains x 16, whereas > I would expect, also given your latest reply, 9x16x16x2=#subdomains x > submatrix height x submatrix width x length of a (row,column) pair. > It would really help me if you could briefly explain how the output above > encodes the subdivision into subdomains. > Many thanks again, > Leonardo > > > > Il giorno mar 9 mag 2023 alle ore 16:24 Matthew Knepley <knep...@gmail.com> > ha scritto: > >> On Tue, May 9, 2023 at 10:05 AM LEONARDO MUTTI < >> leonardo.mutt...@universitadipavia.it> wrote: >> >>> Great thanks! I can now successfully run >>> https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/ksp/tests/ex71f.F90. >>> >>> 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. >>> >>> 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). >>> >>> 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. >>> >>> #include <petsc/finclude/petscmat.h> >>> #include <petsc/finclude/petscksp.h> >>> #include <petsc/finclude/petscpc.h> >>> USE petscmat >>> USE petscksp >>> USE petscpc >>> USE MPI >>> >>> Mat :: A >>> Vec :: b, x >>> PetscInt :: M, I, J, ISLen, NSub >>> PetscMPIInt :: size >>> PetscErrorCode :: ierr >>> PetscScalar :: v >>> KSP :: ksp >>> PC :: pc >>> IS :: subdomains_IS(2), inflated_IS(2) >>> PetscInt,DIMENSION(4) :: indices_first_domain >>> PetscInt,DIMENSION(36) :: indices_second_domain >>> >>> call PetscInitialize(PETSC_NULL_CHARACTER, ierr) >>> call MPI_Comm_size(PETSC_COMM_WORLD, size, ierr) >>> >>> >>> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! >>> ! INTRO: create matrix and right hand side >>> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! >>> >>> WRITE(*,*) "Assembling A,b" >>> >>> M = 8 >>> call MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, >>> & M, M, PETSC_DEFAULT_INTEGER, PETSC_NULL_INTEGER, >>> & PETSC_DEFAULT_INTEGER, PETSC_NULL_INTEGER,A, ierr) >>> DO I=1,M >>> DO J=1,M >>> IF ((I .EQ. J) .AND. (I .LE. 2 )) THEN >>> v = 1 >>> ELSE IF ((I .EQ. J) .AND. (I .GT. 2 )) THEN >>> v = 2 >>> ELSE >>> v = 0 >>> ENDIF >>> call MatSetValue(A, I-1, J-1, v, INSERT_VALUES, ierr) >>> END DO >>> END DO >>> >>> call MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr) >>> call MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr) >>> >>> call VecCreate(PETSC_COMM_WORLD,b,ierr) >>> call VecSetSizes(b, PETSC_DECIDE, M,ierr) >>> call VecSetFromOptions(b,ierr) >>> >>> do I=1,M >>> v = 0.5 >>> call VecSetValue(b,I-1,v, INSERT_VALUES,ierr) >>> end do >>> >>> call VecAssemblyBegin(b,ierr) >>> call VecAssemblyEnd(b,ierr) >>> >>> >>> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! >>> ! FIRST KSP/PC SETUP >>> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! >>> >>> WRITE(*,*) "KSP/PC first setup" >>> >>> call KSPCreate(PETSC_COMM_WORLD, ksp, ierr) >>> call KSPSetOperators(ksp, A, A, ierr) >>> call KSPSetType(ksp, 'preonly', ierr) >>> call KSPGetPC(ksp, pc, ierr) >>> call KSPSetUp(ksp, ierr) >>> call PCSetType(pc, PCGASM, ierr) >>> >>> >>> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! >>> ! GASM, SETTING SUBDOMAINS >>> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! >>> >>> WRITE(*,*) "Setting GASM subdomains" >>> >>> ! Let's create the subdomain IS and inflated_IS >>> ! They are equal if no overlap is present >>> ! They are 1: 0,1,8,9 >>> ! 2: 10,...,15,18,...,23,...,58,...,63 >>> >>> indices_first_domain = [0,1,8,9] ! corresponds to diag(1,1) >>> do I=0,5 >>> do J=0,5 >>> indices_second_domain(I*6+1+J) = 18 + J + 8*I ! corresponds >>> to diag(2,2,...,2) >>> !WRITE(*,*) I*6+1+J, 18 + J + 8*I >>> end do >>> end do >>> >>> ! Convert into IS >>> ISLen = 4 >>> call ISCreateGeneral(PETSC_COMM_WORLD,ISLen,indices_first_domain, >>> & PETSC_COPY_VALUES, subdomains_IS(1), ierr) >>> call ISCreateGeneral(PETSC_COMM_WORLD,ISLen,indices_first_domain, >>> & PETSC_COPY_VALUES, inflated_IS(1), ierr) >>> ISLen = 36 >>> call ISCreateGeneral(PETSC_COMM_WORLD,ISLen,indices_second_domain, >>> & PETSC_COPY_VALUES, subdomains_IS(2), ierr) >>> call ISCreateGeneral(PETSC_COMM_WORLD,ISLen,indices_second_domain, >>> & PETSC_COPY_VALUES, inflated_IS(2), ierr) >>> >>> NSub = 2 >>> call PCGASMSetSubdomains(pc,NSub, >>> & subdomains_IS,inflated_IS,ierr) >>> call PCGASMDestroySubdomains(NSub, >>> & subdomains_IS,inflated_IS,ierr) >>> >>> >>> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! >>> ! GASM: SET SUBSOLVERS >>> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! >>> >>> WRITE(*,*) "Setting subsolvers for GASM" >>> >>> call PCSetUp(pc, ierr) ! should I add this? >>> >>> call PetscOptionsSetValue(PETSC_NULL_OPTIONS, >>> & "-sub_pc_type", "lu", ierr) >>> call PetscOptionsSetValue(PETSC_NULL_OPTIONS, >>> & "-sub_ksp_type", "preonly", ierr) >>> >>> call KSPSetFromOptions(ksp, ierr) >>> call PCSetFromOptions(pc, ierr) >>> >>> >>> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! >>> ! DUMMY SOLUTION: DID IT WORK? >>> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! >>> >>> WRITE(*,*) "Solve" >>> >>> call VecDuplicate(b,x,ierr) >>> call KSPSolve(ksp,b,x,ierr) >>> >>> call MatDestroy(A, ierr) >>> call KSPDestroy(ksp, ierr) >>> call PetscFinalize(ierr) >>> >>> This code is failing in multiple points. At call PCSetUp(pc, ierr) it >>> produces: >>> >>> *[0]PETSC ERROR: Argument out of range* >>> *[0]PETSC ERROR: Scatter indices in ix are out of range* >>> *...* >>> *[0]PETSC ERROR: #1 VecScatterCreate() at >>> ***\src\vec\is\sf\INTERF~1\vscat.c:736* >>> *[0]PETSC ERROR: #2 PCSetUp_GASM() at >>> ***\src\ksp\pc\impls\gasm\gasm.c:433* >>> *[0]PETSC ERROR: #3 PCSetUp() at ***\src\ksp\pc\INTERF~1\precon.c:994* >>> >>> And at call KSPSolve(ksp,b,x,ierr) it produces: >>> >>> *forrtl: severe (157): Program Exception - access violation* >>> >>> >>> The index sets are setup coherently with the outputs of e.g. >>> https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/ksp/tests/output/ex71f_1.out: >>> in particular each element of the matrix A corresponds to a number from 0 >>> to 63. >>> >> >> This is not correct, I believe. The indices are row/col indices, not >> indices into dense blocks, so for >> your example, they are all in [0, 8]. >> >> Thanks, >> >> Matt >> >> >>> Note that each submatrix does not represent some physical subdomain, the >>> subdivision is just at the algebraic level. >>> I thus have the following questions: >>> >>> - is this the correct way of creating the IS objects, given my >>> objective at the beginning of the email? Is the ordering correct? >>> - what am I doing wrong that is generating the above errors? >>> >>> Thanks for the patience and the time. >>> Best, >>> Leonardo >>> >>> Il giorno ven 5 mag 2023 alle ore 18:43 Barry Smith <bsm...@petsc.dev> >>> ha scritto: >>> >>>> >>>> Added in *barry/2023-05-04/add-pcgasm-set-subdomains *see also >>>> https://gitlab.com/petsc/petsc/-/merge_requests/6419 >>>> >>>> Barry >>>> >>>> >>>> On May 4, 2023, at 11:23 AM, LEONARDO MUTTI < >>>> leonardo.mutt...@universitadipavia.it> wrote: >>>> >>>> Thank you for the help. >>>> Adding to my example: >>>> >>>> >>>> * call PCGASMSetSubdomains(pc,NSub, subdomains_IS, >>>> inflated_IS,ierr) call >>>> PCGASMDestroySubdomains(NSub,subdomains_IS,inflated_IS,ierr)* >>>> results in: >>>> >>>> * Error LNK2019 unresolved external symbol PCGASMDESTROYSUBDOMAINS >>>> referenced in function ... * >>>> >>>> * Error LNK2019 unresolved external symbol PCGASMSETSUBDOMAINS >>>> referenced in function ... * >>>> I'm not sure if the interfaces are missing or if I have a compilation >>>> problem. >>>> Thank you again. >>>> Best, >>>> Leonardo >>>> >>>> Il giorno sab 29 apr 2023 alle ore 20:30 Barry Smith <bsm...@petsc.dev> >>>> ha scritto: >>>> >>>>> >>>>> Thank you for the test code. I have a fix in the branch >>>>> barry/2023-04-29/fix-pcasmcreatesubdomains2d >>>>> <https://gitlab.com/petsc/petsc/-/commits/barry/2023-04-29/fix-pcasmcreatesubdomains2d> >>>>> with >>>>> merge request https://gitlab.com/petsc/petsc/-/merge_requests/6394 >>>>> >>>>> The functions did not have proper Fortran stubs and interfaces so I >>>>> had to provide them manually in the new branch. >>>>> >>>>> Use >>>>> >>>>> git fetch >>>>> git checkout barry/2023-04-29/fix-pcasmcreatesubdomains2d >>>>> <https://gitlab.com/petsc/petsc/-/commits/barry/2023-04-29/fix-pcasmcreatesubdomains2d> >>>>> ./configure etc >>>>> >>>>> 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. >>>>> >>>>> Please let us know if you have any later questions. >>>>> >>>>> Barry >>>>> >>>>> >>>>> >>>>> >>>>> On Apr 28, 2023, at 12:07 PM, LEONARDO MUTTI < >>>>> leonardo.mutt...@universitadipavia.it> wrote: >>>>> >>>>> 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: >>>>> >>>>> #include <petsc/finclude/petscmat.h> >>>>> #include <petsc/finclude/petscksp.h> >>>>> #include <petsc/finclude/petscpc.h> >>>>> USE petscmat >>>>> USE petscksp >>>>> USE petscpc >>>>> >>>>> Mat :: A >>>>> PetscInt :: M, NSubx, dof, overlap, NSub >>>>> INTEGER :: I,J >>>>> PetscErrorCode :: ierr >>>>> PetscScalar :: v >>>>> KSP :: ksp >>>>> PC :: pc >>>>> IS :: subdomains_IS, inflated_IS >>>>> >>>>> call PetscInitialize(PETSC_NULL_CHARACTER , ierr) >>>>> >>>>> !-----Create a dummy matrix >>>>> M = 16 >>>>> call MatCreateAIJ(MPI_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, >>>>> & M, M, >>>>> & PETSC_DEFAULT_INTEGER, PETSC_NULL_INTEGER, >>>>> & PETSC_DEFAULT_INTEGER, PETSC_NULL_INTEGER, >>>>> & A, ierr) >>>>> >>>>> DO I=1,M >>>>> DO J=1,M >>>>> v = I*J >>>>> CALL MatSetValue (A,I-1,J-1,v, >>>>> & INSERT_VALUES , ierr) >>>>> END DO >>>>> END DO >>>>> >>>>> call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY , ierr) >>>>> call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY , ierr) >>>>> >>>>> !-----Create KSP and PC >>>>> call KSPCreate(PETSC_COMM_WORLD,ksp, ierr) >>>>> call KSPSetOperators(ksp,A,A, ierr) >>>>> call KSPSetType(ksp,"bcgs",ierr) >>>>> call KSPGetPC(ksp,pc,ierr) >>>>> call KSPSetUp(ksp, ierr) >>>>> call PCSetType(pc,PCGASM, ierr) >>>>> call PCSetUp(pc , ierr) >>>>> >>>>> !-----GASM setup >>>>> NSubx = 4 >>>>> dof = 1 >>>>> overlap = 0 >>>>> >>>>> call PCGASMCreateSubdomains2D(pc, >>>>> & M, M, >>>>> & NSubx, NSubx, >>>>> & dof, overlap, >>>>> & NSub, subdomains_IS, inflated_IS, ierr) >>>>> >>>>> call ISView(subdomains_IS, PETSC_VIEWER_STDOUT_WORLD, ierr) >>>>> >>>>> call KSPDestroy(ksp, ierr) >>>>> call PetscFinalize(ierr) >>>>> >>>>> Running this on one processor, I get NSub = 4. >>>>> If PCASM and PCASMCreateSubdomains2D are used instead, I get NSub = 16 >>>>> as expected. >>>>> Moreover, I get in the end "forrtl: severe (157): Program Exception - >>>>> access violation". So: >>>>> 1) why do I get two different results with ASM, and GASM? >>>>> 2) why do I get access violation and how can I solve this? >>>>> 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: >>>>> >>>>> subroutine PCGASMCreateSubdomains2D(a,b,c,d,e,f,g,h,i,j,z) >>>>> import tPC,tIS >>>>> PC a ! PC >>>>> PetscInt b ! PetscInt >>>>> PetscInt c ! PetscInt >>>>> PetscInt d ! PetscInt >>>>> PetscInt e ! PetscInt >>>>> PetscInt f ! PetscInt >>>>> PetscInt g ! PetscInt >>>>> PetscInt h ! PetscInt >>>>> IS i ! IS >>>>> IS j ! IS >>>>> PetscErrorCode z >>>>> end subroutine PCGASMCreateSubdomains2D >>>>> Thus: >>>>> 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? >>>>> >>>>> 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. >>>>> >>>>> Thanks in advance, >>>>> Leonardo >>>>> >>>>> >>>>> >>>> >> >> -- >> What most experimenters take for granted before they begin their >> experiments is infinitely more interesting than any results to which their >> experiments lead. >> -- Norbert Wiener >> >> https://www.cse.buffalo.edu/~knepley/ >> <http://www.cse.buffalo.edu/~knepley/> >> >