<div dir="ltr"><div dir="ltr"><font face="arial, sans-serif">Many thanks for the hint. As for your last two sentences, I am in a similar situation. So, l</font><font face="arial, sans-serif">et me mention some last thoughts, in case they help locate the issue</font>. If not, I will embark on the debugging adventure with Windows debuggers.<div><font face="arial, sans-serif"><br></font></div><div><font face="arial, sans-serif">Consider again the 4x4 matrix (represented by the IS [0,1,2,3], with rows 0,1 on rank 0 and 2,3 on rank 1). If I create the </font><span style="font-family:arial,sans-serif">non-inflated IS [0,1,2], [3], with </span><span style="font-family:arial,sans-serif">inflated IS [0,1,2],
[1,2,3], no problems arise.</span></div><div><font face="arial, sans-serif">Here is a comparison with the non-working case:</font></div><div><ul><li><font face="arial, sans-serif">works: non-inflated IS: [0,1,2], [3], inflated IS: [0,1,2], [1,2,3]</font></li><li><font face="arial, sans-serif">doesn't work: non-inflated IS: [0,1], [2,3], inflated IS: [0,1], [1,2,3] </font></li></ul><div><font face="arial, sans-serif">One difference is that in the working case, the non-inflating subdomain [0,1,2], is not on a single process.</font></div></div><div><font face="arial, sans-serif"><br></font></div><div><font face="arial, sans-serif">Moreover, is it possible that the issue has something to do with </font><span style="font-family:monospace">PetscObjectReference</span><font face="arial, sans-serif"> as seen in </font>/src/ksp/pc/impls/gasm/gasm.c line 1732, PCGASMCreateSubdomains2D?</div><div><br></div></div><div>Thanks again<br></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">Il giorno gio 18 mag 2023 alle ore 04:26 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> Yikes. Such huge numbers usually come from integer overflow or memory corruption. <div><br></div><div>The code to decide on the memory that needs allocating is straightforward</div><div><br></div><div><div>PetscErrorCode MatCreateSubMatrices_MPIAIJ(Mat C, PetscInt ismax, const IS isrow[], const IS iscol[], MatReuse scall, Mat *submat[])</div><div>{</div><div> PetscInt nmax, nstages = 0, i, pos, max_no, nrow, ncol, in[2], out[2];</div><div> PetscBool rowflag, colflag, wantallmatrix = PETSC_FALSE;</div><div> Mat_SeqAIJ *subc;</div><div> Mat_SubSppt *smat;</div><div><br></div><div> PetscFunctionBegin;</div><div> /* Check for special case: each processor has a single IS */</div><div> if (C->submat_singleis) { /* flag is set in PCSetUp_ASM() to skip MPI_Allreduce() */</div><div> PetscCall(MatCreateSubMatrices_MPIAIJ_SingleIS(C, ismax, isrow, iscol, scall, submat));</div><div> C->submat_singleis = PETSC_FALSE; /* resume its default value in case C will be used for non-single IS */</div><div> PetscFunctionReturn(PETSC_SUCCESS);</div><div> }</div><div><br></div><div> /* Collect global wantallmatrix and nstages */</div><div> if (!C->cmap->N) nmax = 20 * 1000000 / sizeof(PetscInt);</div><div> else nmax = 20 * 1000000 / (C->cmap->N * sizeof(PetscInt));</div><div> if (!nmax) nmax = 1;</div><div><br></div><div> if (scall == MAT_INITIAL_MATRIX) {</div><div> /* Collect global wantallmatrix and nstages */</div><div> if (ismax == 1 && C->rmap->N == C->cmap->N) {</div><div> PetscCall(ISIdentity(*isrow, &rowflag));</div><div> PetscCall(ISIdentity(*iscol, &colflag));</div><div> PetscCall(ISGetLocalSize(*isrow, &nrow));</div><div> PetscCall(ISGetLocalSize(*iscol, &ncol));</div><div> if (rowflag && colflag && nrow == C->rmap->N && ncol == C->cmap->N) {</div><div> wantallmatrix = PETSC_TRUE;</div><div><br></div><div> PetscCall(PetscOptionsGetBool(((PetscObject)C)->options, ((PetscObject)C)->prefix, "-use_fast_submatrix", &wantallmatrix, NULL));</div><div> }</div><div> }</div><div><br></div><div> /* Determine the number of stages through which submatrices are done</div><div> Each stage will extract nmax submatrices.</div><div> nmax is determined by the matrix column dimension.</div><div> If the original matrix has 20M columns, only one submatrix per stage is allowed, etc.</div><div> */</div><div> nstages = ismax / nmax + ((ismax % nmax) ? 1 : 0); /* local nstages */</div><div><br></div><div> in[0] = -1 * (PetscInt)wantallmatrix;</div><div> in[1] = nstages;</div><div> PetscCall(MPIU_Allreduce(in, out, 2, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)C)));</div><div> wantallmatrix = (PetscBool)(-out[0]);</div><div> nstages = out[1]; /* Make sure every processor loops through the global nstages */</div><div><br></div><div> } else { /* MAT_REUSE_MATRIX */</div><div> if (ismax) {</div><div> subc = (Mat_SeqAIJ *)(*submat)[0]->data;</div><div> smat = subc->submatis1;</div><div> } else { /* (*submat)[0] is a dummy matrix */</div><div> smat = (Mat_SubSppt *)(*submat)[0]->data;</div><div> }</div><div> if (!smat) {</div><div> /* smat is not generated by MatCreateSubMatrix_MPIAIJ_All(...,MAT_INITIAL_MATRIX,...) */</div><div> wantallmatrix = PETSC_TRUE;</div><div> } else if (smat->singleis) {</div><div> PetscCall(MatCreateSubMatrices_MPIAIJ_SingleIS(C, ismax, isrow, iscol, scall, submat));</div><div> PetscFunctionReturn(PETSC_SUCCESS);</div><div> } else {</div><div> nstages = smat->nstages;</div><div> }</div><div> }</div><div><br></div><div> if (wantallmatrix) {</div><div> PetscCall(MatCreateSubMatrix_MPIAIJ_All(C, MAT_GET_VALUES, scall, submat));</div><div> PetscFunctionReturn(PETSC_SUCCESS);</div><div> }</div><div><br></div><div> /* Allocate memory to hold all the submatrices and dummy submatrices */</div><div> if (scall == MAT_INITIAL_MATRIX) PetscCall(PetscCalloc1(ismax + nstages, submat));</div><div><br></div><div> for (i = 0, pos = 0; i < nstages; i++) {</div><div> if (pos + nmax <= ismax) max_no = nmax;</div><div> else if (pos >= ismax) max_no = 0;</div><div> else max_no = ismax - pos;</div><div><br></div><div> PetscCall(MatCreateSubMatrices_MPIAIJ_Local(C, max_no, isrow + pos, iscol + pos, scall, *submat + pos));</div><div> if (!max_no) {</div><div> if (scall == MAT_INITIAL_MATRIX) { /* submat[pos] is a dummy matrix */</div><div> smat = (Mat_SubSppt *)(*submat)[pos]->data;</div><div> smat->nstages = nstages;</div><div> }</div><div> pos++; /* advance to next dummy matrix if any */</div><div> } else pos += max_no;</div><div> }</div><div><br></div><div> if (ismax && scall == MAT_INITIAL_MATRIX) {</div><div> /* save nstages for reuse */</div><div> subc = (Mat_SeqAIJ *)(*submat)[0]->data;</div><div> smat = subc->submatis1;</div><div> smat->nstages = nstages;</div><div> }</div><div> PetscFunctionReturn(PETSC_SUCCESS);</div><div>}</div><div><br></div><div>The easiest way to debug would be to put a breakpoint in MatCreateSubMatrices_MPIAIJ on MPI rank zero and next through the subroutine to see where the crazy number appears that gets passed down in the line if (scall == MAT_INITIAL_MATRIX) PetscCall(PetscCalloc1(ismax + nstages, submat)); where either ismax or nstages has a crazy value.</div><div><br></div><div>If you are using the GNU compilers you can use the command line options -start_in_debugger noxterm -debugger_ranks 0 to start the Gnu debugger. If you are using the Microsoft Windows compilers you will need to use their debugger, I don't know how to do that (and shudder at the thought :-).</div><div><br><blockquote type="cite"><div>On May 17, 2023, at 7:51 PM, Leonardo Mutti <<a href="mailto:leonardo.mutti01@universitadipavia.it" target="_blank">leonardo.mutti01@universitadipavia.it</a>> wrote:</div><br><div><div dir="ltr"><table cellpadding="0" id="m_7983612177498689413gmail-undefined" style="width:1266.4px;table-layout:fixed;border-spacing:0px"><tbody><tr><td style="width:1266.4px;vertical-align:top;height:100px"><div id="m_7983612177498689413gmail-:wr" style="border-radius:1px;box-sizing:border-box;padding:0px 0px 16px;border:0px;margin:0px"><div></div></div></td></tr></tbody></table>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.<div>What I am doing wrong and what should be done differently, this is still unclear to me.</div><div><br></div><div>The simplified code runs on 2 processors, I built a 4x4 matrix.</div><div>The subdomains are now given by [0,1] and [2,3], with [2,3] inflating to [0,1,2,3].</div><div><br></div><div><div>Thank you again.</div></div><div><br></div><div>Error:</div><div><br></div><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"><div><font face="monospace">[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------</font></div><div><font face="monospace">[0]PETSC ERROR: Out of memory. This could be due to allocating</font></div><div><font face="monospace">[0]PETSC ERROR: too large an object or bleeding by not properly</font></div><div><font face="monospace">[0]PETSC ERROR: destroying unneeded objects.</font></div><div><font face="monospace">[0]PETSC ERROR: Memory allocated 0 Memory used by process 0</font></div><div><font face="monospace">[0]PETSC ERROR: Try running with -malloc_dump or -malloc_view for info.</font></div><div><font face="monospace">[0]PETSC ERROR: Memory requested 9437902811936987136</font></div><div><font face="monospace">[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!</font></div><div><font face="monospace">[0]PETSC ERROR: Option left: name:-pc_gasm_view_subdomains value: 1 source: code</font></div><div><font face="monospace">[0]PETSC ERROR: Option left: name:-sub_ksp_type value: gmres source: code</font></div><div><font face="monospace">[0]PETSC ERROR: Option left: name:-sub_pc_type value: none source: code</font></div><div><font face="monospace">[0]PETSC ERROR: See <a href="https://petsc.org/release/faq/" target="_blank">https://petsc.org/release/faq/</a> for trouble shooting.</font></div><div><font face="monospace">[0]PETSC ERROR: Petsc Development GIT revision: v3.19.1-234-g43977f8d16 GIT Date: 2023-05-08 14:50:03 +0000</font></div><div><font face="monospace">[...]</font></div><div><font face="monospace">[0]PETSC ERROR: #1 PetscMallocAlign() at ...\Sources\Git\Test-FV\PETSC-~1\src\sys\memory\mal.c:66</font></div><div><font face="monospace">[0]PETSC ERROR: #2 PetscMallocA() at ...\Sources\Git\Test-FV\PETSC-~1\src\sys\memory\mal.c:411</font></div><div><font face="monospace">[0]PETSC ERROR: #3 MatCreateSubMatrices_MPIAIJ() at ...\Sources\Git\Test-FV\PETSC-~1\src\mat\impls\aij\mpi\mpiov.c:2025</font></div><div><font face="monospace">[0]PETSC ERROR: #4 MatCreateSubMatricesMPI_MPIXAIJ() at ...\Sources\Git\Test-FV\PETSC-~1\src\mat\impls\aij\mpi\mpiov.c:3136</font></div><div><font face="monospace">[0]PETSC ERROR: #5 MatCreateSubMatricesMPI_MPIAIJ() at ...\Sources\Git\Test-FV\PETSC-~1\src\mat\impls\aij\mpi\mpiov.c:3208</font></div><div><font face="monospace">[0]PETSC ERROR: #6 MatCreateSubMatricesMPI() at ...\Sources\Git\Test-FV\PETSC-~1\src\mat\INTERF~1\matrix.c:7071</font></div><div><font face="monospace">[0]PETSC ERROR: #7 PCSetUp_GASM() at ...\Sources\Git\Test-FV\PETSC-~1\src\ksp\pc\impls\gasm\gasm.c:556</font></div><div><font face="monospace">[0]PETSC ERROR: #8 PCSetUp() at ...\Sources\Git\Test-FV\PETSC-~1\src\ksp\pc\INTERF~1\precon.c:994</font></div><div><font face="monospace">[0]PETSC ERROR: #9 KSPSetUp() at ...\Sources\Git\Test-FV\PETSC-~1\src\ksp\ksp\INTERF~1\itfunc.c:406</font></div></blockquote><div><font face="arial, sans-serif"><br></font></div><div><font face="arial, sans-serif">Code, the important bit is after </font><span style="color:rgb(0,255,0);font-family:monospace">! GASM, SETTING SUBDOMAINS</span><font face="arial, sans-serif">:</font></div><div><span style="color:rgb(0,255,0)"><font face="arial, sans-serif"><br></font></span></div><div><font face="monospace"> <font color="#cccccc"> Mat :: A<br> Vec :: b<br> PetscInt :: M,N_blocks,block_size,NSub,I<br> PetscErrorCode :: ierr<br> PetscScalar :: v<br> KSP :: ksp<br> PC :: pc<br> IS :: subdomains_IS(2), inflated_IS(2)<br> PetscInt :: NMPI,MYRANK,IERMPI<br> INTEGER :: start<br> <br> call PetscInitialize(PETSC_NULL_CHARACTER, ierr)<br> call PetscLogDefaultBegin(ierr)<br> call MPI_COMM_SIZE(MPI_COMM_WORLD, NMPI, IERMPI)<br> CALL MPI_COMM_RANK(MPI_COMM_WORLD, MYRANK,IERMPI) <br> <br> N_blocks = 2<br> block_size = 2<br> M = N_blocks * block_size<br> <br><br> ! INTRO: create matrix and right hand side, create IS </font></font></div><div><font face="monospace" color="#cccccc"><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> call VecCreate(PETSC_COMM_WORLD,b,ierr)<br> call VecSetSizes(b, PETSC_DECIDE, M,ierr)<br> call VecSetFromOptions(b,ierr)<br><br> DO I=(MYRANK*(M/NMPI)),((MYRANK+1)*(M/NMPI)-1)<br> <br> ! Set matrix<br> v=1<br> call MatSetValue(A, I, I, v, INSERT_VALUES, ierr)<br> IF (I-block_size .GE. 0) THEN<br> v=-1<br> call MatSetValue(A, I, I-block_size, v, INSERT_VALUES, ierr)<br> ENDIF<br> ! Set rhs<br> v = I<br> call VecSetValue(b,I,v, INSERT_VALUES,ierr)<br> <br> END DO <br> <br> call MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr)<br> call MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr)<br> call VecAssemblyBegin(b,ierr)<br> call VecAssemblyEnd(b,ierr)<br><br> ! FIRST KSP/PC SETUP</font></div><div><font face="monospace"><font color="#cccccc"><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 PCSetType(pc, PCGASM, ierr)<br><br></font><br><font color="#00ff00"> ! GASM, SETTING SUBDOMAINS<br></font> <br> if (myrank == 0) then<br> call ISCreateStride(PETSC_COMM_SELF, 2, 0, 1, subdomains_IS(1), ierr)<br> call ISCreateStride(PETSC_COMM_WORLD, 0, 0, 1, subdomains_IS(2), ierr)<br> call ISCreateStride(PETSC_COMM_SELF, 2, 0, 1, inflated_IS(1), ierr)<br> call ISCreateStride(PETSC_COMM_WORLD, 2, 0, 1, inflated_IS(2), ierr)<br> start = 1<br> NSub = 2<br> else<br> call ISCreateStride(PETSC_COMM_WORLD, 2, 2, 1, subdomains_IS(2), ierr)<br> call ISCreateStride(PETSC_COMM_WORLD, 2, 2, 1, inflated_IS(2), ierr)<br> start = 2<br> NSub = 1<br> endif<br><br> call PCGASMSetSubdomains(pc,NSub,subdomains_IS(start:2),inflated_IS(start:2),ierr)<br> call PCGASMDestroySubdomains(NSub,subdomains_IS(start:2),inflated_IS(start:2),ierr)<br> <br><font color="#cccccc"> ! GASM: SET SUBSOLVERS<br><br> call PetscOptionsSetValue(PETSC_NULL_OPTIONS,"-sub_ksp_type", "gmres", ierr)<br> call PetscOptionsSetValue(PETSC_NULL_OPTIONS,"-sub_pc_type", "none", ierr) <br> call PetscOptionsSetValue(PETSC_NULL_OPTIONS,"-pc_gasm_view_subdomains", "1", ierr)<br> <br> call KSPSetUp(ksp, ierr)<br> call PCSetUp(pc, ierr)<br> call KSPSetFromOptions(ksp, ierr)<br> call PCSetFromOptions(pc, ierr)<br> call KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD, ierr)<br><br> call MatDestroy(A, ierr)<br> call PetscFinalize(ierr) </font> </font><br></div><div><font face="monospace"><br></font></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">Il giorno mer 17 mag 2023 alle ore 21:55 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><br><div><br><blockquote type="cite"><div>On May 17, 2023, at 11:10 AM, Leonardo Mutti <<a href="mailto:leonardo.mutti01@universitadipavia.it" target="_blank">leonardo.mutti01@universitadipavia.it</a>> wrote:</div><br><div><div dir="ltr"><div><font face="arial, sans-serif">Dear developers, let me kindly ask for your help again. </font></div><div><font face="arial, sans-serif">In the following snippet, </font><font face="arial, sans-serif">a bi-diagonal matrix A is set up. It measures 8x8 blocks, each block is 2x2 elements. </font><span style="font-family:arial,sans-serif">I would like to create the correct IS objects for PCGASM. </span></div><div><font face="arial, sans-serif">The non-overlapping IS should be: [<b>0,1</b>], [<b>2,3</b>],[<b>4,5</b>], ..., [<b>14,15</b>]. </font><span style="font-family:arial,sans-serif">The overlapping IS should be: [</span><b style="font-family:arial,sans-serif">0,1</b><span style="font-family:arial,sans-serif">], [0,1,</span><b style="font-family:arial,sans-serif">2,3</b><span style="font-family:arial,sans-serif">], [2,3,</span><b style="font-family:arial,sans-serif">4,5</b><span style="font-family:arial,sans-serif">], ..., [12,13,</span><b style="font-family:arial,sans-serif">14,15</b><span style="font-family:arial,sans-serif">]</span></div><div><font face="arial, sans-serif">I am running the code with 4 processors. </font><font face="arial, sans-serif">For some reason, after calling </font><span style="font-family:monospace">PCGASMDestroySubdomains </span><font face="arial, sans-serif">the code crashes with </font><font face="monospace">severe (157): Program Exception - access violation</font><font face="arial, sans-serif">. A visual inspection of the indices using </font><font face="monospace">ISView </font><font face="arial, sans-serif">looks good.</font></div></div></div></blockquote><div><br></div> 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.</div><div><br></div><div><br><blockquote type="cite"><div><div dir="ltr"><div><font face="arial, sans-serif">Thanks again,</font></div><div><font face="arial, sans-serif">Leonardo</font></div><div dir="ltr"><font face="monospace"><font color="#eeeeee"><br></font></font></div><div dir="ltr"><font face="monospace"><font color="#eeeeee"> Mat :: A<br> Vec :: b<br> PetscInt :: M,N_blocks,block_size,I,J,NSub,converged_reason,srank,erank,color,subcomm<br> PetscMPIInt :: size<br> PetscErrorCode :: ierr<br> PetscScalar :: v<br> KSP :: ksp<br> PC :: pc<br> IS,ALLOCATABLE :: subdomains_IS(:), inflated_IS(:)<br> PetscInt :: NMPI,MYRANK,IERMPI<br> INTEGER :: IS_counter, is_start, is_end<br> <br> call PetscInitialize(PETSC_NULL_CHARACTER, ierr)<br> call PetscLogDefaultBegin(ierr)<br> call MPI_COMM_SIZE(MPI_COMM_WORLD, NMPI, IERMPI)<br> CALL MPI_COMM_RANK(MPI_COMM_WORLD, MYRANK,IERMPI) <br> <br> N_blocks = 8<br> block_size = 2<br> M = N_blocks * block_size<br> <br> ALLOCATE(subdomains_IS(N_blocks))<br> ALLOCATE(inflated_IS(N_blocks))<br> <br> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!<br> ! ASSUMPTION: no block spans more than one rank (the inflated blocks can)<br> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! <br> <br> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!<br> ! INTRO: create matrix and right hand side<br> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!<br> <br> ! How many inflated blocks span more than one rank? NMPI-1 !<br> <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> call VecCreate(PETSC_COMM_WORLD,b,ierr)<br> call VecSetSizes(b, PETSC_DECIDE, M,ierr)<br> call VecSetFromOptions(b,ierr)<br><br> DO I=(MYRANK*(M/NMPI)),((MYRANK+1)*(M/NMPI)-1)<br> <br> ! Set matrix<br> v=1<br> call MatSetValue(A, I, I, v, INSERT_VALUES, ierr)<br> IF (I-block_size .GE. 0) THEN<br> v=-1<br> call MatSetValue(A, I, I-block_size, v, INSERT_VALUES, ierr)<br> ENDIF<br> ! Set rhs<br> v = I<br> call VecSetValue(b,I,v, INSERT_VALUES,ierr)<br> <br> END DO <br> <br> call MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr)<br> call MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr)<br> call VecAssemblyBegin(b,ierr)<br> call VecAssemblyEnd(b,ierr)<br><br> <br> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!<br> ! FIRST KSP/PC SETUP<br> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!<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 PCSetType(pc, PCGASM, ierr)</font></font><div><font face="monospace"><br><br><font color="#00ff00"> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!<br> !! GASM, SETTING SUBDOMAINS<br> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!<br></font></font><span style="font-family:monospace"><br></span></div><div><span style="font-family:monospace"> DO IS_COUNTER=1,N_blocks</span><font face="monospace"><br></font></div><div><font face="monospace"> <br> srank = MAX(((IS_COUNTER-2)*block_size)/(M/NMPI),0) <font color="#00ff00">! start rank reached by inflated block</font><br> erank = MIN(((IS_COUNTER-1)*block_size)/(M/NMPI),NMPI-1)<font color="#00ff00"> ! end rank reached by inflated block. Coincides with rank containing non-inflated block</font></font></div><div><font face="monospace"><br></font></div></div><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"><div><div><font face="monospace"> <font color="#00ff00">! Create subcomms</font></font></div></div><div><div><font face="monospace"> color = MPI_UNDEFINED</font></div></div><div><div><font face="monospace"> IF (myrank == srank .or. myrank == erank) THEN</font></div></div><div><div><font face="monospace"> color = 1</font></div></div><div><div><font face="monospace"> ENDIF</font></div></div><div><div><font face="monospace"> call MPI_Comm_split(MPI_COMM_WORLD,color,MYRANK,subcomm,ierr) </font></div></div></blockquote><div dir="ltr"><div><font face="monospace"><br></font></div><div><font face="monospace"> <font color="#00ff00">! Create IS</font><br> IF (srank .EQ. erank) THEN </font>
<span style="font-family:monospace"><font color="#00ff00">! Block and overlap are on the same rank</font></span><font face="monospace"><br> IF (MYRANK .EQ. srank) THEN <br> call ISCreateStride(PETSC_COMM_SELF,block_size,(IS_COUNTER-1)*block_size,1,subdomains_IS(IS_COUNTER),ierr)<br> IF (IS_COUNTER .EQ. 1) THEN <font color="#00ff00"> ! the first block is not inflated</font><br> </font>
<span style="font-family:monospace">call ISCreateStride(PETSC_COMM_SELF,block_size,(IS_COUNTER-1)*block_size,1,inflated_IS(IS_COUNTER),ierr)</span><font face="monospace"><br> ELSE<br> call ISCreateStride(PETSC_COMM_SELF,2*block_size,(IS_COUNTER-2)*block_size,1,inflated_IS(IS_COUNTER),ierr)<br> ENDIF<br> ENDIF<br> else</font><font color="#00ff00"> <span style="font-family:monospace">! Block and overlap not on the same rank</span></font><font face="monospace"><br> if (myrank == erank) then<font color="#00ff00"> ! the block</font><br> call </font>
<span style="font-family:monospace">ISCreateStride</span><font face="monospace">(subcomm,block_size,(IS_COUNTER-1)*block_size,1,subdomains_IS(IS_COUNTER),ierr)<br> call </font>
<span style="font-family:monospace">ISCreateStride</span><font face="monospace">(subcomm,block_size,(IS_COUNTER-1)*block_size,1,inflated_IS(IS_COUNTER),ierr)<br> endif<br> if (myrank == srank) then<font color="#00ff00"> ! the overlap</font><br> call </font>
<span style="font-family:monospace">ISCreateStride</span><font face="monospace">(subcomm,block_size,(IS_COUNTER-2)*block_size,1,inflated_IS(IS_COUNTER),ierr)<br> call </font>
<span style="font-family:monospace">ISCreateStride</span><font face="monospace">(subcomm,0,(IS_COUNTER-1)*block_size,1,subdomains_IS(IS_COUNTER),ierr)<br> endif<br> endif <br> <br> call MPI_Comm_free(subcomm, ierr)<br> END DO<br> <br><font color="#00ff00"> ! Set the domains/subdomains</font></font><div><span style="font-family:monospace"> NSub = N_blocks/NMPI</span><br style="font-family:monospace"><font face="monospace"> is_start = 1 + myrank * NSub<br> is_end = min(is_start + NSub, N_blocks)<br> if (myrank + 1 < NMPI) then <br> NSub = NSub + 1<br> endif <br> <br> call PCGASMSetSubdomains(pc,NSub,subdomains_IS(is_start:is_end),inflated_IS(is_start:is_end),ierr)<br> call PCGASMDestroySubdomains(NSub,subdomains_IS(is_start:is_end),inflated_IS(is_start:is_end),ierr)<br> <br><font color="#eeeeee"> call PetscOptionsSetValue(PETSC_NULL_OPTIONS,"-sub_ksp_type", "gmres", ierr)<br> call PetscOptionsSetValue(PETSC_NULL_OPTIONS,"-sub_pc_type", "none", ierr) <br> call PetscOptionsSetValue(PETSC_NULL_OPTIONS,"-pc_gasm_view_subdomains", "1", ierr)<br> <br> call KSPSetUp(ksp, ierr)<br> call PCSetUp(pc, ierr)<br> call KSPSetFromOptions(ksp, ierr)<br> call PCSetFromOptions(pc, ierr)<br> <br> call KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD, ierr)<br></font> <br></font></div></div></div><font face="monospace"><br></font><div class="gmail_quote"><div dir="ltr" class="gmail_attr"><font face="monospace">Il giorno mer 10 mag 2023 alle ore 03:02 Barry Smith <<a href="mailto:bsmith@petsc.dev" target="_blank">bsmith@petsc.dev</a>> ha scritto:<br></font></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div><font face="monospace"><br></font><div><font face="monospace"><br></font><blockquote type="cite"><div><font face="monospace">On May 9, 2023, at 4:58 PM, LEONARDO MUTTI <<a href="mailto:leonardo.mutti01@universitadipavia.it" target="_blank">leonardo.mutti01@universitadipavia.it</a>> wrote:</font></div><font face="monospace"><br></font><div><div dir="auto"><div><font face="monospace">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). </font></div><div dir="auto"><font face="monospace">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?</font></div></div></div></blockquote><div><font face="monospace"><br></font></div><font face="monospace">I believe so</font></div><div><font face="monospace"><br></font><blockquote type="cite"><div><div dir="auto"><div dir="auto"><font face="monospace">Thanks,</font></div><div dir="auto"><font face="monospace">Leonardo</font></div><div dir="auto"><font face="monospace"><br></font><div class="gmail_quote" dir="auto"><div dir="ltr" class="gmail_attr"><font face="monospace">Il mar 9 mag 2023, 20:45 Barry Smith <<a href="mailto:bsmith@petsc.dev" target="_blank">bsmith@petsc.dev</a>> ha scritto:<br></font></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><font face="monospace"><br></font></div><font face="monospace"> 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.</font><div><font face="monospace"><br></font></div><div><font face="monospace"> 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</font></div><div><font face="monospace"><br></font></div><div><font face="monospace"> I do not understand why your code has </font></div><div><font face="monospace"><br></font></div><div><blockquote type="cite"><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div class="gmail_quote"><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><font face="monospace">indices_first_domain = [0,1,8,9] ! corresponds to diag(1,1)</font></div></blockquote></div></div></blockquote></div></div></div></blockquote></div></blockquote><font face="monospace"><br></font></div><div><font face="monospace">it should just be 0</font></div><div><font face="monospace"><br></font></div><div><font face="monospace"><br></font><div><font face="monospace"><br></font></div><div><font face="monospace"><br></font><div><font face="monospace"><br></font><blockquote type="cite"><div><font face="monospace">On May 9, 2023, at 12:44 PM, LEONARDO MUTTI <<a href="mailto:leonardo.mutti01@universitadipavia.it" rel="noreferrer" target="_blank">leonardo.mutti01@universitadipavia.it</a>> wrote:</font></div><font face="monospace"><br></font><div><div dir="ltr"><font face="monospace">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).</font></div><font face="monospace"><br></font><div class="gmail_quote"><div dir="ltr" class="gmail_attr"><font face="monospace">Il giorno mar 9 mag 2023 alle ore 18:31 LEONARDO MUTTI <<a href="mailto:leonardo.mutti01@universitadipavia.it" rel="noreferrer" target="_blank">leonardo.mutti01@universitadipavia.it</a>> ha scritto:<br></font></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><font face="monospace"><br><br></font><div class="gmail_quote"><div dir="ltr" class="gmail_attr"><font face="monospace">---------- Forwarded message ---------<br>Da: <strong class="gmail_sendername" dir="auto">LEONARDO MUTTI</strong> <span dir="auto"><<a href="mailto:leonardo.mutti01@universitadipavia.it" rel="noreferrer" target="_blank">leonardo.mutti01@universitadipavia.it</a>></span><br>Date: mar 9 mag 2023 alle ore 18:29<br>Subject: Re: [petsc-users] Understanding index sets for PCGASM<br>To: Matthew Knepley <<a href="mailto:knepley@gmail.com" rel="noreferrer" target="_blank">knepley@gmail.com</a>><br></font></div><font face="monospace"><br><br></font><div dir="ltr"><font face="monospace">Thank you for your answer, but I am still confused, sorry. </font><div><font face="monospace">Consider <a href="https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/ksp/tests/ex71f.F90" rel="noreferrer" target="_blank">https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/ksp/tests/ex71f.F90</a> on one processor.</font></div><div><font face="monospace">Let M=12 for the sake of simplicity, i.e. we deal with a 12x12 2D grid, hence, a 144x144 matrix.</font></div><div><font face="monospace">Let NSubx = 3, so that on the grid we do 3 vertical and 3 horizontal subdivisions.</font></div><div><font face="monospace">We should obtain 9 subdomains that are grids of 4x4 nodes each, thus corresponding to 9 submatrices of size 16x16. </font></div><div><div><font face="monospace">In my run I obtain NSub = 9 (great) and subdomain_IS(i), i=1,...,9, reads:</font></div></div><div><font face="monospace"><br></font></div><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"><div><font face="monospace"><i>IS Object: 1 MPI process</i><br></font></div><div><i><font face="monospace"> type: general</font></i></div><div><i><font face="monospace">Number of indices in set 16</font></i></div><div><i><font face="monospace">0 0</font></i></div><div><i><font face="monospace">1 1</font></i></div><div><i><font face="monospace">2 2</font></i></div><div><i><font face="monospace">3 3</font></i></div><div><i><font face="monospace">4 12</font></i></div><div><i><font face="monospace">5 13</font></i></div><div><i><font face="monospace">6 14</font></i></div><div><i><font face="monospace">7 15</font></i></div><div><i><font face="monospace">8 24</font></i></div><div><i><font face="monospace">9 25</font></i></div><div><i><font face="monospace">10 26</font></i></div><div><i><font face="monospace">11 27</font></i></div><div><i><font face="monospace">12 36</font></i></div><div><i><font face="monospace">13 37</font></i></div><div><i><font face="monospace">14 38</font></i></div><div><i><font face="monospace">15 39</font></i></div><div><i><font face="monospace">IS Object: 1 MPI process</font></i></div><div><i><font face="monospace"> type: general</font></i></div><div><i><font face="monospace">Number of indices in set 16</font></i></div><div><i><font face="monospace">0 4</font></i></div><div><i><font face="monospace">1 5</font></i></div><div><i><font face="monospace">2 6</font></i></div><div><i><font face="monospace">3 7</font></i></div><div><i><font face="monospace">4 16</font></i></div><div><i><font face="monospace">5 17</font></i></div><div><i><font face="monospace">6 18</font></i></div><div><i><font face="monospace">7 19</font></i></div><div><i><font face="monospace">8 28</font></i></div><div><i><font face="monospace">9 29</font></i></div><div><i><font face="monospace">10 30</font></i></div><div><i><font face="monospace">11 31</font></i></div><div><i><font face="monospace">12 40</font></i></div><div><i><font face="monospace">13 41</font></i></div><div><i><font face="monospace">14 42</font></i></div><div><i><font face="monospace">15 43</font></i></div><div><i><font face="monospace">IS Object: 1 MPI process</font></i></div><div><i><font face="monospace"> type: general</font></i></div><div><i><font face="monospace">Number of indices in set 16</font></i></div><div><i><font face="monospace">0 8</font></i></div><div><i><font face="monospace">1 9</font></i></div><div><i><font face="monospace">2 10</font></i></div><div><i><font face="monospace">3 11</font></i></div><div><i><font face="monospace">4 20</font></i></div><div><i><font face="monospace">5 21</font></i></div><div><i><font face="monospace">6 22</font></i></div><div><i><font face="monospace">7 23</font></i></div><div><i><font face="monospace">8 32</font></i></div><div><i><font face="monospace">9 33</font></i></div><div><i><font face="monospace">10 34</font></i></div><div><i><font face="monospace">11 35</font></i></div><div><i><font face="monospace">12 44</font></i></div><div><i><font face="monospace">13 45</font></i></div><div><i><font face="monospace">14 46</font></i></div><div><i><font face="monospace">15 47</font></i></div><div><i><font face="monospace">IS Object: 1 MPI process</font></i></div><div><i><font face="monospace"> type: general</font></i></div><div><i><font face="monospace">Number of indices in set 16</font></i></div><div><i><font face="monospace">0 48</font></i></div><div><i><font face="monospace">1 49</font></i></div><div><i><font face="monospace">2 50</font></i></div><div><i><font face="monospace">3 51</font></i></div><div><i><font face="monospace">4 60</font></i></div><div><i><font face="monospace">5 61</font></i></div><div><i><font face="monospace">6 62</font></i></div><div><i><font face="monospace">7 63</font></i></div><div><i><font face="monospace">8 72</font></i></div><div><i><font face="monospace">9 73</font></i></div><div><i><font face="monospace">10 74</font></i></div><div><i><font face="monospace">11 75</font></i></div><div><i><font face="monospace">12 84</font></i></div><div><i><font face="monospace">13 85</font></i></div><div><i><font face="monospace">14 86</font></i></div><div><i><font face="monospace">15 87</font></i></div><div><i><font face="monospace">IS Object: 1 MPI process</font></i></div><div><i><font face="monospace"> type: general</font></i></div><div><i><font face="monospace">Number of indices in set 16</font></i></div><div><i><font face="monospace">0 52</font></i></div><div><i><font face="monospace">1 53</font></i></div><div><i><font face="monospace">2 54</font></i></div><div><i><font face="monospace">3 55</font></i></div><div><i><font face="monospace">4 64</font></i></div><div><i><font face="monospace">5 65</font></i></div><div><i><font face="monospace">6 66</font></i></div><div><i><font face="monospace">7 67</font></i></div><div><i><font face="monospace">8 76</font></i></div><div><i><font face="monospace">9 77</font></i></div><div><i><font face="monospace">10 78</font></i></div><div><i><font face="monospace">11 79</font></i></div><div><i><font face="monospace">12 88</font></i></div><div><i><font face="monospace">13 89</font></i></div><div><i><font face="monospace">14 90</font></i></div><div><i><font face="monospace">15 91</font></i></div><div><i><font face="monospace">IS Object: 1 MPI process</font></i></div><div><i><font face="monospace"> type: general</font></i></div><div><i><font face="monospace">Number of indices in set 16</font></i></div><div><i><font face="monospace">0 56</font></i></div><div><i><font face="monospace">1 57</font></i></div><div><i><font face="monospace">2 58</font></i></div><div><i><font face="monospace">3 59</font></i></div><div><i><font face="monospace">4 68</font></i></div><div><i><font face="monospace">5 69</font></i></div><div><i><font face="monospace">6 70</font></i></div><div><i><font face="monospace">7 71</font></i></div><div><i><font face="monospace">8 80</font></i></div><div><i><font face="monospace">9 81</font></i></div><div><i><font face="monospace">10 82</font></i></div><div><i><font face="monospace">11 83</font></i></div><div><i><font face="monospace">12 92</font></i></div><div><i><font face="monospace">13 93</font></i></div><div><i><font face="monospace">14 94</font></i></div><div><i><font face="monospace">15 95</font></i></div><div><i><font face="monospace">IS Object: 1 MPI process</font></i></div><div><i><font face="monospace"> type: general</font></i></div><div><i><font face="monospace">Number of indices in set 16</font></i></div><div><i><font face="monospace">0 96</font></i></div><div><i><font face="monospace">1 97</font></i></div><div><i><font face="monospace">2 98</font></i></div><div><i><font face="monospace">3 99</font></i></div><div><i><font face="monospace">4 108</font></i></div><div><i><font face="monospace">5 109</font></i></div><div><i><font face="monospace">6 110</font></i></div><div><i><font face="monospace">7 111</font></i></div><div><i><font face="monospace">8 120</font></i></div><div><i><font face="monospace">9 121</font></i></div><div><i><font face="monospace">10 122</font></i></div><div><i><font face="monospace">11 123</font></i></div><div><i><font face="monospace">12 132</font></i></div><div><i><font face="monospace">13 133</font></i></div><div><i><font face="monospace">14 134</font></i></div><div><i><font face="monospace">15 135</font></i></div><div><i><font face="monospace">IS Object: 1 MPI process</font></i></div><div><i><font face="monospace"> type: general</font></i></div><div><i><font face="monospace">Number of indices in set 16</font></i></div><div><i><font face="monospace">0 100</font></i></div><div><i><font face="monospace">1 101</font></i></div><div><i><font face="monospace">2 102</font></i></div><div><i><font face="monospace">3 103</font></i></div><div><i><font face="monospace">4 112</font></i></div><div><i><font face="monospace">5 113</font></i></div><div><i><font face="monospace">6 114</font></i></div><div><i><font face="monospace">7 115</font></i></div><div><i><font face="monospace">8 124</font></i></div><div><i><font face="monospace">9 125</font></i></div><div><i><font face="monospace">10 126</font></i></div><div><i><font face="monospace">11 127</font></i></div><div><i><font face="monospace">12 136</font></i></div><div><i><font face="monospace">13 137</font></i></div><div><i><font face="monospace">14 138</font></i></div><div><i><font face="monospace">15 139</font></i></div><div><i><font face="monospace">IS Object: 1 MPI process</font></i></div><div><i><font face="monospace"> type: general</font></i></div><div><i><font face="monospace">Number of indices in set 16</font></i></div><div><i><font face="monospace">0 104</font></i></div><div><i><font face="monospace">1 105</font></i></div><div><i><font face="monospace">2 106</font></i></div><div><i><font face="monospace">3 107</font></i></div><div><i><font face="monospace">4 116</font></i></div><div><i><font face="monospace">5 117</font></i></div><div><i><font face="monospace">6 118</font></i></div><div><i><font face="monospace">7 119</font></i></div><div><i><font face="monospace">8 128</font></i></div><div><i><font face="monospace">9 129</font></i></div><div><i><font face="monospace">10 130</font></i></div><div><i><font face="monospace">11 131</font></i></div><div><i><font face="monospace">12 140</font></i></div><div><i><font face="monospace">13 141</font></i></div><div><i><font face="monospace">14 142</font></i></div><div><i><font face="monospace">15 143</font></i></div><div><font face="monospace"><br></font></div></blockquote><div><font face="monospace">As you said, no number here reaches 144. </font></div><div><font face="monospace">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.</font></div><div><font face="monospace">It would really help me if you could briefly explain how the output above encodes the subdivision into subdomains.<br></font></div><div><font face="monospace">Many thanks again,<br></font></div><div><font face="monospace">Leonardo<br></font><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"><div><font face="monospace"><br></font></div></blockquote></div></div><font face="monospace"><br></font><div class="gmail_quote"><div dir="ltr" class="gmail_attr"><font face="monospace">Il giorno mar 9 mag 2023 alle ore 16:24 Matthew Knepley <<a href="mailto:knepley@gmail.com" rel="noreferrer" target="_blank">knepley@gmail.com</a>> ha scritto:<br></font></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div dir="ltr"><font face="monospace">On Tue, May 9, 2023 at 10:05 AM LEONARDO MUTTI <<a href="mailto:leonardo.mutti01@universitadipavia.it" rel="noreferrer" target="_blank">leonardo.mutti01@universitadipavia.it</a>> wrote:<br></font></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div><font face="monospace">Great thanks! I can now successfully run <a href="https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/ksp/tests/ex71f.F90" rel="noreferrer" target="_blank">https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/ksp/tests/ex71f.F90</a>.</font></div><div><font face="monospace"><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.</font></div><div><font face="monospace"><br></font></div><div><font face="monospace">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).</font></div><div><font face="monospace"><br></font></div><div><font face="monospace">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.</font></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)<br></font></div><div><font face="monospace"><br></font></div><div><font face="monospace">This code is failing in multiple points. At call PCSetUp(pc, ierr) it produces:</font></div><div><font face="monospace"><br></font></div><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"><div><i><font color="#ff0000" face="monospace">[0]PETSC ERROR: Argument out of range</font></i></div><div><i><font color="#ff0000" face="monospace">[0]PETSC ERROR: Scatter indices in ix are out of range</font></i></div><div><i><font color="#ff0000" face="monospace">...</font></i></div><div><i><font color="#ff0000" face="monospace">[0]PETSC ERROR: #1 VecScatterCreate() at ***\src\vec\is\sf\INTERF~1\vscat.c:736</font></i></div><div><i><font color="#ff0000" face="monospace">[0]PETSC ERROR: #2 PCSetUp_GASM() at ***\src\ksp\pc\impls\gasm\gasm.c:433</font></i></div><div><i><font color="#ff0000" face="monospace">[0]PETSC ERROR: #3 PCSetUp() at ***\src\ksp\pc\INTERF~1\precon.c:994</font></i></div><div><i><font face="monospace"><br></font></i></div></blockquote><div><font face="monospace">And at call KSPSolve(ksp,b,x,ierr) it produces:</font></div><div><font face="monospace"><br></font></div><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"><div><i><font color="#ff0000" face="monospace">forrtl: severe (157): Program Exception - access violation</font></i></div></blockquote><div><font face="monospace"><br></font></div><div><font face="monospace">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" rel="noreferrer" target="_blank">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.</font></div></div></blockquote><div><font face="monospace"><br></font></div><div><font face="monospace">This is not correct, I believe. The indices are row/col indices, not indices into dense blocks, so for</font></div><div><font face="monospace">your example, they are all in [0, 8].</font></div><div><font face="monospace"><br></font></div><div><font face="monospace"> Thanks,</font></div><div><font face="monospace"><br></font></div><div><font face="monospace"> Matt</font></div><div><font face="monospace"> </font></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div><font face="monospace">Note that each submatrix does not represent some physical subdomain, the subdivision is just at the algebraic level.</font></div><div><font face="monospace">I thus have the following questions:</font></div><div><ul><li><font face="monospace">is this the correct way of creating the IS objects, given my objective at the beginning of the email? Is the ordering correct?</font></li><li><font face="monospace">what am I doing wrong that is generating the above errors?</font></li></ul><div><font face="monospace">Thanks for the patience and the time.</font></div><div><font face="monospace">Best, <br>Leonardo</font></div></div><div class="gmail_quote"><div dir="ltr" class="gmail_attr"><font face="monospace"><br></font></div><div dir="ltr" class="gmail_attr"><font face="monospace">Il giorno ven 5 mag 2023 alle ore 18:43 Barry Smith <<a href="mailto:bsmith@petsc.dev" rel="noreferrer" target="_blank">bsmith@petsc.dev</a>> ha scritto:<br></font></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><font face="monospace"><br></font></div><font face="monospace"> Added in <b style="color:rgb(200,20,201);font-size:14px">barry/2023-05-04/add-pcgasm-set-subdomains </b>see also <a href="https://gitlab.com/petsc/petsc/-/merge_requests/6419" rel="noreferrer" target="_blank">https://gitlab.com/petsc/petsc/-/merge_requests/6419</a></font><div><font face="monospace"><br></font></div><div><font face="monospace"> Barry</font></div><div><font face="monospace"><br></font><div><font face="monospace"><br></font><blockquote type="cite"><div><font face="monospace">On May 4, 2023, at 11:23 AM, LEONARDO MUTTI <<a href="mailto:leonardo.mutti01@universitadipavia.it" rel="noreferrer" target="_blank">leonardo.mutti01@universitadipavia.it</a>> wrote:</font></div><font face="monospace"><br></font><div><div dir="ltr"><font face="monospace">Thank you for the help. </font><div><font face="monospace">Adding to my example:</font><div><i><font face="monospace"> call PCGASMSetSubdomains(pc,NSub, subdomains_IS, inflated_IS,ierr)<br> call PCGASMDestroySubdomains(NSub,subdomains_IS,inflated_IS,ierr)<br></font></i></div><div><font face="monospace">results in:<br></font></div><div><div><i><font face="monospace"> Error LNK2019 unresolved external symbol PCGASMDESTROYSUBDOMAINS referenced in function ... <br></font></i></div><div><i><font face="monospace"> Error LNK2019 unresolved external symbol PCGASMSETSUBDOMAINS referenced in function ... <br></font></i></div><div><font face="monospace">I'm not sure if the interfaces are missing or if I have a compilation problem.</font></div></div><div><font face="monospace">Thank you again. </font></div><div><font face="monospace">Best,<br>Leonardo</font></div></div></div><font face="monospace"><br></font><div class="gmail_quote"><div dir="ltr" class="gmail_attr"><font face="monospace">Il giorno sab 29 apr 2023 alle ore 20:30 Barry Smith <<a href="mailto:bsmith@petsc.dev" rel="noreferrer" target="_blank">bsmith@petsc.dev</a>> ha scritto:<br></font></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><font face="monospace"><br></font></div><font face="monospace"> Thank you for the test code. I have a fix in the branch <span style="color:rgb(51,50,56);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" rel="noreferrer" target="_blank">barry/2023-04-29/fix-pcasmcreatesubdomains2d</a> with merge request <a href="https://gitlab.com/petsc/petsc/-/merge_requests/6394" rel="noreferrer" target="_blank">https://gitlab.com/petsc/petsc/-/merge_requests/6394</a></font><div><font face="monospace"><br></font></div><div><font face="monospace"> The functions did not have proper Fortran stubs and interfaces so I had to provide them manually in the new branch.<br></font><div><font face="monospace"><br></font></div><div><font face="monospace"> Use</font></div><div><font face="monospace"><br></font></div><div><font face="monospace"> git fetch</font></div><div><font face="monospace"> 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" rel="noreferrer" target="_blank">barry/2023-04-29/fix-pcasmcreatesubdomains2d</a></font></div><div><font face="monospace"> ./configure etc</font></div><div><font face="monospace"><br></font></div><div><font face="monospace"> 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.</font></div><div><font face="monospace"><br></font></div><div><font face="monospace"> Please let us know if you have any later questions.</font></div><div><font face="monospace"><br></font></div><div><font face="monospace"> Barry</font></div><div><font face="monospace"><br></font></div><div><font face="monospace"><br></font></div><div><font face="monospace"><br></font></div><div><div><font face="monospace"><br></font><blockquote type="cite"><div><font face="monospace">On Apr 28, 2023, at 12:07 PM, LEONARDO MUTTI <<a href="mailto:leonardo.mutti01@universitadipavia.it" rel="noreferrer" target="_blank">leonardo.mutti01@universitadipavia.it</a>> wrote:</font></div><font face="monospace"><br></font><div><div dir="ltr"><font face="monospace">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:</font><div><font face="monospace"><br></font></div><div><font face="monospace">#include <petsc/finclude/petscmat.h><br></font></div><div><font face="monospace">#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)</font></div><div><font face="monospace"> <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> </font></div><div><font face="monospace">!-----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)</font></div><div><font face="monospace"><br></font></div><div><font face="monospace">Running this on one processor, I get NSub = 4.</font></div><div><font face="monospace">If PCASM and PCASMCreateSubdomains2D are used instead, I get NSub = 16 as expected.</font></div><div><font face="monospace">Moreover, I get in the end "forrtl: severe (157): Program Exception - access violation". So:</font></div><div><font face="monospace">1) why do I get two different results with ASM, and GASM?</font></div><div><font face="monospace">2) why do I get access violation and how can I solve this?</font></div><div><font face="monospace">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></font></div><div><font face="monospace"><br></font></div><div><font face="monospace"> 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></font></div><div><font face="monospace">Thus:</font></div><div><font face="monospace">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?</font></div><div><font face="monospace"><br></font></div><div><font face="monospace">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.</font></div><div><font face="monospace"><br></font></div><div><font face="monospace">Thanks in advance,</font></div><div><font face="monospace">Leonardo</font></div></div>
</div></blockquote></div><font face="monospace"><br></font></div></div></div></blockquote></div>
</div></blockquote></div><font face="monospace"><br></font></div></div></blockquote></div></div>
</blockquote></div><font face="monospace"><br clear="all"></font><div><font face="monospace"><br></font></div><font face="monospace"><span>-- </span><br></font><div dir="ltr"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div><font face="monospace">What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</font></div><div><font face="monospace"><br></font></div><div><font face="monospace"><a href="http://www.cse.buffalo.edu/~knepley/" rel="noreferrer" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></font></div></div></div></div></div></div></div></div>
</blockquote></div>
</div></div>
</blockquote></div>
</div></blockquote></div><font face="monospace"><br></font></div></div></div></blockquote></div></div></div>
</div></blockquote></div><br></div></blockquote></div></div>
</div></blockquote></div><br></div></blockquote></div>
</div></blockquote></div><br></div></div></blockquote></div></div>