Thanks for the reply. Even without Valgrind (which I can't use since I'm on Windows), by further simplifying the example, I was able to have PETSc display a more informative message. What I am doing wrong and what should be done differently, this is still unclear to me.
The simplified code runs on 2 processors, I built a 4x4 matrix. The subdomains are now given by [0,1] and [2,3], with [2,3] inflating to [0,1,2,3]. Thank you again. Error: [0]PETSC ERROR: --------------------- Error Message -------------------------------------------------------------- [0]PETSC ERROR: Out of memory. This could be due to allocating [0]PETSC ERROR: too large an object or bleeding by not properly [0]PETSC ERROR: destroying unneeded objects. [0]PETSC ERROR: Memory allocated 0 Memory used by process 0 [0]PETSC ERROR: Try running with -malloc_dump or -malloc_view for info. [0]PETSC ERROR: Memory requested 9437902811936987136 [0]PETSC ERROR: WARNING! There are option(s) set that were not used! Could be the program crashed before they were used or a spelling mistake, etc! [0]PETSC ERROR: Option left: name:-pc_gasm_view_subdomains value: 1 source: code [0]PETSC ERROR: Option left: name:-sub_ksp_type value: gmres source: code [0]PETSC ERROR: Option left: name:-sub_pc_type value: none source: code [0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting. [0]PETSC ERROR: Petsc Development GIT revision: v3.19.1-234-g43977f8d16 GIT Date: 2023-05-08 14:50:03 +0000 [...] [0]PETSC ERROR: #1 PetscMallocAlign() at ...\Sources\Git\Test-FV\PETSC-~1\src\sys\memory\mal.c:66 [0]PETSC ERROR: #2 PetscMallocA() at ...\Sources\Git\Test-FV\PETSC-~1\src\sys\memory\mal.c:411 [0]PETSC ERROR: #3 MatCreateSubMatrices_MPIAIJ() at ...\Sources\Git\Test-FV\PETSC-~1\src\mat\impls\aij\mpi\mpiov.c:2025 [0]PETSC ERROR: #4 MatCreateSubMatricesMPI_MPIXAIJ() at ...\Sources\Git\Test-FV\PETSC-~1\src\mat\impls\aij\mpi\mpiov.c:3136 [0]PETSC ERROR: #5 MatCreateSubMatricesMPI_MPIAIJ() at ...\Sources\Git\Test-FV\PETSC-~1\src\mat\impls\aij\mpi\mpiov.c:3208 [0]PETSC ERROR: #6 MatCreateSubMatricesMPI() at ...\Sources\Git\Test-FV\PETSC-~1\src\mat\INTERF~1\matrix.c:7071 [0]PETSC ERROR: #7 PCSetUp_GASM() at ...\Sources\Git\Test-FV\PETSC-~1\src\ksp\pc\impls\gasm\gasm.c:556 [0]PETSC ERROR: #8 PCSetUp() at ...\Sources\Git\Test-FV\PETSC-~1\src\ksp\pc\INTERF~1\precon.c:994 [0]PETSC ERROR: #9 KSPSetUp() at ...\Sources\Git\Test-FV\PETSC-~1\src\ksp\ksp\INTERF~1\itfunc.c:406 Code, the important bit is after ! GASM, SETTING SUBDOMAINS: Mat :: A Vec :: b PetscInt :: M,N_blocks,block_size,NSub,I PetscErrorCode :: ierr PetscScalar :: v KSP :: ksp PC :: pc IS :: subdomains_IS(2), inflated_IS(2) PetscInt :: NMPI,MYRANK,IERMPI INTEGER :: start call PetscInitialize(PETSC_NULL_CHARACTER, ierr) call PetscLogDefaultBegin(ierr) call MPI_COMM_SIZE(MPI_COMM_WORLD, NMPI, IERMPI) CALL MPI_COMM_RANK(MPI_COMM_WORLD, MYRANK,IERMPI) N_blocks = 2 block_size = 2 M = N_blocks * block_size ! INTRO: create matrix and right hand side, create IS 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) call VecCreate(PETSC_COMM_WORLD,b,ierr) call VecSetSizes(b, PETSC_DECIDE, M,ierr) call VecSetFromOptions(b,ierr) DO I=(MYRANK*(M/NMPI)),((MYRANK+1)*(M/NMPI)-1) ! Set matrix v=1 call MatSetValue(A, I, I, v, INSERT_VALUES, ierr) IF (I-block_size .GE. 0) THEN v=-1 call MatSetValue(A, I, I-block_size, v, INSERT_VALUES, ierr) ENDIF ! Set rhs v = I call VecSetValue(b,I,v, INSERT_VALUES,ierr) END DO call MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr) call MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr) call VecAssemblyBegin(b,ierr) call VecAssemblyEnd(b,ierr) ! FIRST KSP/PC 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 PCSetType(pc, PCGASM, ierr) ! GASM, SETTING SUBDOMAINS if (myrank == 0) then call ISCreateStride(PETSC_COMM_SELF, 2, 0, 1, subdomains_IS(1), ierr) call ISCreateStride(PETSC_COMM_WORLD, 0, 0, 1, subdomains_IS(2), ierr) call ISCreateStride(PETSC_COMM_SELF, 2, 0, 1, inflated_IS(1), ierr) call ISCreateStride(PETSC_COMM_WORLD, 2, 0, 1, inflated_IS(2), ierr) start = 1 NSub = 2 else call ISCreateStride(PETSC_COMM_WORLD, 2, 2, 1, subdomains_IS(2), ierr) call ISCreateStride(PETSC_COMM_WORLD, 2, 2, 1, inflated_IS(2), ierr) start = 2 NSub = 1 endif call PCGASMSetSubdomains(pc,NSub,subdomains_IS(start:2),inflated_IS(start:2),ierr) call PCGASMDestroySubdomains(NSub,subdomains_IS(start:2),inflated_IS(start:2),ierr) ! GASM: SET SUBSOLVERS call PetscOptionsSetValue(PETSC_NULL_OPTIONS,"-sub_ksp_type", "gmres", ierr) call PetscOptionsSetValue(PETSC_NULL_OPTIONS,"-sub_pc_type", "none", ierr) call PetscOptionsSetValue(PETSC_NULL_OPTIONS,"-pc_gasm_view_subdomains", "1", ierr) call KSPSetUp(ksp, ierr) call PCSetUp(pc, ierr) call KSPSetFromOptions(ksp, ierr) call PCSetFromOptions(pc, ierr) call KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD, ierr) call MatDestroy(A, ierr) call PetscFinalize(ierr) Il giorno mer 17 mag 2023 alle ore 21:55 Barry Smith <[email protected]> ha scritto: > > > On May 17, 2023, at 11:10 AM, Leonardo Mutti < > [email protected]> wrote: > > Dear developers, let me kindly ask for your help again. > In the following snippet, a bi-diagonal matrix A is set up. It measures > 8x8 blocks, each block is 2x2 elements. I would like to create the > correct IS objects for PCGASM. > The non-overlapping IS should be: [*0,1*], [*2,3*],[*4,5*], ..., [*14,15* > ]. The overlapping IS should be: [*0,1*], [0,1,*2,3*], [2,3,*4,5*], ..., > [12,13,*14,15*] > I am running the code with 4 processors. For some reason, after calling > PCGASMDestroySubdomains > the code crashes with severe (157): Program Exception - access violation. > A visual inspection of the indices using ISView looks good. > > > Likely memory corruption or use of an object or an array that was > already freed. Best to use Valgrind to find the exact location of the mess. > > > Thanks again, > Leonardo > > Mat :: A > Vec :: b > PetscInt :: > M,N_blocks,block_size,I,J,NSub,converged_reason,srank,erank,color,subcomm > PetscMPIInt :: size > PetscErrorCode :: ierr > PetscScalar :: v > KSP :: ksp > PC :: pc > IS,ALLOCATABLE :: subdomains_IS(:), inflated_IS(:) > PetscInt :: NMPI,MYRANK,IERMPI > INTEGER :: IS_counter, is_start, is_end > > call PetscInitialize(PETSC_NULL_CHARACTER, ierr) > call PetscLogDefaultBegin(ierr) > call MPI_COMM_SIZE(MPI_COMM_WORLD, NMPI, IERMPI) > CALL MPI_COMM_RANK(MPI_COMM_WORLD, MYRANK,IERMPI) > > N_blocks = 8 > block_size = 2 > M = N_blocks * block_size > > ALLOCATE(subdomains_IS(N_blocks)) > ALLOCATE(inflated_IS(N_blocks)) > > !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! > ! ASSUMPTION: no block spans more than one rank (the inflated blocks > can) > !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! > > !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! > ! INTRO: create matrix and right hand side > !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! > > ! How many inflated blocks span more than one rank? NMPI-1 ! > > 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) > call VecCreate(PETSC_COMM_WORLD,b,ierr) > call VecSetSizes(b, PETSC_DECIDE, M,ierr) > call VecSetFromOptions(b,ierr) > > DO I=(MYRANK*(M/NMPI)),((MYRANK+1)*(M/NMPI)-1) > > ! Set matrix > v=1 > call MatSetValue(A, I, I, v, INSERT_VALUES, ierr) > IF (I-block_size .GE. 0) THEN > v=-1 > call MatSetValue(A, I, I-block_size, v, INSERT_VALUES, ierr) > ENDIF > ! Set rhs > v = I > call VecSetValue(b,I,v, INSERT_VALUES,ierr) > > END DO > > call MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr) > call MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr) > call VecAssemblyBegin(b,ierr) > call VecAssemblyEnd(b,ierr) > > > !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! > ! FIRST KSP/PC 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 PCSetType(pc, PCGASM, ierr) > > > !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! > !! GASM, SETTING SUBDOMAINS > !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! > > DO IS_COUNTER=1,N_blocks > > srank = MAX(((IS_COUNTER-2)*block_size)/(M/NMPI),0) ! start rank > reached by inflated block > erank = MIN(((IS_COUNTER-1)*block_size)/(M/NMPI),NMPI-1) ! end > rank reached by inflated block. Coincides with rank containing non-inflated > block > > ! Create subcomms > color = MPI_UNDEFINED > IF (myrank == srank .or. myrank == erank) THEN > color = 1 > ENDIF > call MPI_Comm_split(MPI_COMM_WORLD,color,MYRANK,subcomm,ierr) > > > ! Create IS > IF (srank .EQ. erank) THEN ! Block and overlap are on the same > rank > IF (MYRANK .EQ. srank) THEN > call > ISCreateStride(PETSC_COMM_SELF,block_size,(IS_COUNTER-1)*block_size,1,subdomains_IS(IS_COUNTER),ierr) > IF (IS_COUNTER .EQ. 1) THEN ! the first block is not > inflated > call > ISCreateStride(PETSC_COMM_SELF,block_size,(IS_COUNTER-1)*block_size,1,inflated_IS(IS_COUNTER),ierr) > ELSE > call > ISCreateStride(PETSC_COMM_SELF,2*block_size,(IS_COUNTER-2)*block_size,1,inflated_IS(IS_COUNTER),ierr) > ENDIF > ENDIF > else ! Block and overlap not on the same rank > if (myrank == erank) then ! the block > call ISCreateStride > (subcomm,block_size,(IS_COUNTER-1)*block_size,1,subdomains_IS(IS_COUNTER),ierr) > call ISCreateStride > (subcomm,block_size,(IS_COUNTER-1)*block_size,1,inflated_IS(IS_COUNTER),ierr) > endif > if (myrank == srank) then ! the overlap > call ISCreateStride > (subcomm,block_size,(IS_COUNTER-2)*block_size,1,inflated_IS(IS_COUNTER),ierr) > call ISCreateStride > (subcomm,0,(IS_COUNTER-1)*block_size,1,subdomains_IS(IS_COUNTER),ierr) > endif > endif > > call MPI_Comm_free(subcomm, ierr) > END DO > > ! Set the domains/subdomains > NSub = N_blocks/NMPI > is_start = 1 + myrank * NSub > is_end = min(is_start + NSub, N_blocks) > if (myrank + 1 < NMPI) then > NSub = NSub + 1 > endif > > call > PCGASMSetSubdomains(pc,NSub,subdomains_IS(is_start:is_end),inflated_IS(is_start:is_end),ierr) > call > PCGASMDestroySubdomains(NSub,subdomains_IS(is_start:is_end),inflated_IS(is_start:is_end),ierr) > > call PetscOptionsSetValue(PETSC_NULL_OPTIONS,"-sub_ksp_type", > "gmres", ierr) > call PetscOptionsSetValue(PETSC_NULL_OPTIONS,"-sub_pc_type", "none", > ierr) > call > PetscOptionsSetValue(PETSC_NULL_OPTIONS,"-pc_gasm_view_subdomains", "1", > ierr) > > call KSPSetUp(ksp, ierr) > call PCSetUp(pc, ierr) > call KSPSetFromOptions(ksp, ierr) > call PCSetFromOptions(pc, ierr) > > call KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD, ierr) > > > Il giorno mer 10 mag 2023 alle ore 03:02 Barry Smith <[email protected]> > ha scritto: > >> >> >> On May 9, 2023, at 4:58 PM, LEONARDO MUTTI < >> [email protected]> wrote: >> >> In my notation diag(1,1) means a diagonal 2x2 matrix with 1,1 on the >> diagonal, submatrix in the 8x8 diagonal matrix diag(1,1,2,2,...,2). >> Am I then correct that the IS representing diag(1,1) is 0,1, and that >> diag(2,2,...,2) is represented by 2,3,4,5,6,7? >> >> >> I believe so >> >> Thanks, >> Leonardo >> >> Il mar 9 mag 2023, 20:45 Barry Smith <[email protected]> ha scritto: >> >>> >>> It is simplier than you are making it out to be. Each IS[] is a list >>> of rows (and columns) in the sub (domain) matrix. In your case with the >>> matrix of 144 by 144 the indices will go from 0 to 143. >>> >>> In your simple Fortran code you have a completely different problem. A >>> matrix with 8 rows and columns. In that case if you want the first IS to >>> represent just the first row (and column) in the matrix then it should >>> contain only 0. The second submatrix which is all rows (but the first) >>> should have 1,2,3,4,5,6,7 >>> >>> I do not understand why your code has >>> >>> indices_first_domain = [0,1,8,9] ! corresponds to diag(1,1) >>>>>> >>>>> >>> it should just be 0 >>> >>> >>> >>> >>> >>> On May 9, 2023, at 12:44 PM, LEONARDO MUTTI < >>> [email protected]> wrote: >>> >>> 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 < >>> [email protected]> ha scritto: >>> >>>> >>>> >>>> ---------- Forwarded message --------- >>>> Da: LEONARDO MUTTI <[email protected]> >>>> Date: mar 9 mag 2023 alle ore 18:29 >>>> Subject: Re: [petsc-users] Understanding index sets for PCGASM >>>> To: Matthew Knepley <[email protected]> >>>> >>>> >>>> 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 < >>>> [email protected]> ha scritto: >>>> >>>>> On Tue, May 9, 2023 at 10:05 AM LEONARDO MUTTI < >>>>> [email protected]> 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 <[email protected]> >>>>>> 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 < >>>>>>> [email protected]> 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 < >>>>>>> [email protected]> 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 < >>>>>>>> [email protected]> 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/> >>>>> >>>> >>> >> >
