---------- 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/> >
