The code I sent compile and run at my side with petsc-3.7.6 (on debian/testing with gcc-6.3). The code you sent does not compile at my side. Anyway, no big deal.
The modification you propose as far as I understand is to replace "ISCreateGeneral(PETSC_COMM_WORLD" with "ISCreateGeneral(PETSC_COMM_SELF" : still not working at my side (empty dirichlet local matrix). I will try to get that with a MPI matrix (that would contain same data that MatIS : that's what I tried to avoid as this doubles allocations - anyway, no big deal). Franck ----- Mail original ----- > De: "Stefano Zampini" <[email protected]> > À: "Franck Houssen" <[email protected]> > Cc: "petsc-dev" <[email protected]>, "PETSc users list" > <[email protected]>, "petsc-maint" <[email protected]> > Envoyé: Mardi 23 Mai 2017 20:23:49 > Objet: Re: [petsc-dev] Using PETSc MatIS, how to get local matrix (= one > domain) before and after assembly ? > > On May 23, 2017, at 6:34 PM, Franck Houssen < [email protected] > > > wrote: > > > OK. I am supposed to destroy the matrix returned by MatISGetMPIXAIJ ? > > Yes > > Also, my example still not get the final assembled local matrix (the > > MatCreateSubMatrix returns an empty matrix) but as far as I understand my > > (global) index set is OK: what did I miss ? > > I really doubt you can use the example you have sent. It doesn’t compile, as > MatCreateSubMatrix needs an extra argument. > Attached a modified version that does what I guess is what you are looking > for (sequential Dirichlet problems on the subdomains). > > Franck > > > ----- Mail original ----- > > > > De: "Stefano Zampini" < [email protected] > > > > > > > À: "Franck Houssen" < [email protected] > > > > > > > Cc: "petsc-dev" < [email protected] >, "PETSc users list" < > > > [email protected] >, "petsc-maint" < [email protected] > > > > > > > Envoyé: Mardi 23 Mai 2017 13:16:18 > > > > > > Objet: Re: [petsc-dev] Using PETSc MatIS, how to get local matrix (= one > > > domain) before and after assembly ? > > > > > > MatISGetMPIXAIJ is collective, as it assembles the global operator. To > > > get > > > the matrices you are looking for, you should call MatCreateSubMatrix on > > > the > > > assembled global operator, with the global indices representing the > > > subdomain problem. Each process needs to call both functions > > > > > > Stefano > > > > > > Il 23 Mag 2017 11:41, "Franck Houssen" < [email protected] > ha > > > scritto: > > > > > > > I have a 3x3 global matrix made of two overlapping 2x2 local matrix (= > > > > diagonal with 1.). Each local matrix correspond to one domain (each > > > > domain > > > > is delegated to one MPI proc, so, I have 2 MPI procs because I have 2 > > > > domains). > > > > > > > > > > This is the simplest possible example: I have two 2x2 (local) diag > > > > matrix > > > > that overlap so that the global matrix built from them is 1, 2, 1 on > > > > the > > > > diagonal (local contributions add up in the middle). > > > > > > > > > > Now, I need for each MPI proc to get the assembled local matrix > > > > (sometimes > > > > called the dirichlet matrix) : this is a local matrix (sequential - not > > > > distributed with MPI) that accounts for contribution of neighboring > > > > domains > > > > (MPI proc). > > > > > > > > > > How to get the local assembled matrix ? MatGetLocalSubMatrix does not > > > > work > > > > (throw error - see example attached). MatGetSubMatrix returns a MPI > > > > distributed matrix, not a local (sequential) one. > > > > > > > > > > 1. My understanding is that MatISGetMPIXAIJ should return a local > > > > matrix > > > > (sequential AIJ matrix) : the MPI in the name recall that you get the > > > > assembled matrix (with contributions from the shared border) from the > > > > other > > > > MPI processus. Correct ? In my simple example, I replaced > > > > MatGetLocalSubMatrix with MatISGetMPIXAIJ : I get a deadlock which was > > > > surprising to me... Is MatISGetMPIXAIJ a collective call ? > > > > > > > > > > 2. Supposing this is a collective call (and that point 1 is not > > > > correct), > > > > I > > > > ride up MatISGetMPIXAIJ before the "if (rank > 0)" : I don't deadlock > > > > now, > > > > but it seems I get a global matrix which is not the assembled local > > > > matrix > > > > I > > > > am looking for. > > > > > > > > > > 3. I am supposed to destroy the matrix returned by MatISGetMPIXAIJ ? (I > > > > believe yes - not sure as AFAIU wording should associate Destroy > > > > methods > > > > to > > > > Create methods) > > > > > > > > > > Franck > > > > > > > > > > The git diff illustrate modifications I tried to add to the initial > > > > file > > > > attached to this thread: > > > > > > > > > > --- a/matISLocalMat.cpp > > > > > > > > > > +++ b/matISLocalMat.cpp > > > > > > > > > > @@ -31,6 +31,8 @@ int main(int argc,char **argv) { > > > > > > > > > > MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(A, > > > > MAT_FINAL_ASSEMBLY); > > > > > > > > > > MatView(A, PETSC_VIEWER_STDOUT_WORLD); > > > > PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD); // Diag: 1, 2, 1 > > > > > > > > > > + Mat assembledLocalMat; > > > > > > > > > > + MatISGetMPIXAIJ(A, MAT_INITIAL_MATRIX, &assembledLocalMat); > > > > > > > > > > if (rank > 0) { // Do not pollute stdout: print only 1 proc > > > > > > > > > > std::cout << std::endl << "non assembled local matrix:" << std::endl << > > > > std::endl; > > > > > > > > > > Mat nonAssembledLocalMat; > > > > > > > > > > @@ -38,11 +40,10 @@ int main(int argc,char **argv) { > > > > > > > > > > MatView(nonAssembledLocalMat, PETSC_VIEWER_STDOUT_SELF); // Diag: 1, 1 > > > > > > > > > > std::cout << std::endl << "assembled local matrix:" << std::endl << > > > > std::endl; > > > > > > > > > > - Mat assembledLocalMat; > > > > > > > > > > - IS is; ISCreateGeneral(PETSC_COMM_SELF, localSize, localIdx, > > > > PETSC_COPY_VALUES, &is); > > > > > > > > > > - MatGetLocalSubMatrix(A, is, is, &assembledLocalMat); // KO ?!... > > > > > > > > > > - MatView(assembledLocalMat, PETSC_VIEWER_STDOUT_SELF); // Would like > > > > to > > > > get > > > > => Diag: 2, 1 > > > > > > > > > > + //IS is; ISCreateGeneral(PETSC_COMM_SELF, localSize, localIdx, > > > > PETSC_COPY_VALUES, &is); > > > > > > > > > > + //MatGetLocalSubMatrix(A, is, is, &assembledLocalMat); // KO ?!... > > > > > > > > > > } > > > > > > > > > > + MatView(assembledLocalMat, PETSC_VIEWER_STDOUT_WORLD); // Would like > > > > to > > > > get > > > > => Diag: 2, 1 > > > > > > > > > > > De: "Stefano Zampini" < [email protected] > > > > > > > > > > > > > > > > À: "petsc-maint" < [email protected] > > > > > > > > > > > > > > > > Cc: "petsc-dev" < [email protected] >, "PETSc users list" < > > > > > [email protected] >, "Franck Houssen" < [email protected] > > > > > > > > > > > > > > > > > > > > > Envoyé: Dimanche 21 Mai 2017 22:51:34 > > > > > > > > > > > > > > > Objet: Re: [petsc-dev] Using PETSc MatIS, how to get local matrix (= > > > > > one > > > > > domain) before and after assembly ? > > > > > > > > > > > > > > > To assemble the operator in aij format, use > > > > > > > > > > > > > > > MatISGetMPIXAIJ > > > > > > > > > > > > > > > http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatISGetMPIXAIJ.html > > > > > > > > > > > > > > > Il 21 Mag 2017 18:43, "Matthew Knepley" < [email protected] > ha > > > > > scritto: > > > > > > > > > > > > > > > > On Sun, May 21, 2017 at 11:23 AM, Franck Houssen < > > > > > > [email protected] > > > > > > > > > > > > > wrote: > > > > > > > > > > > > > > > > > > > > > > I have a 3x3 global matrix is built (diag: 1, 2, 1): it's made of > > > > > > > 2 > > > > > > > overlapping 2x2 local matrix (diag: 1, 1). > > > > > > > > > > > > > > > > > > > > > > > > > > > > Getting non assembled local matrix is OK with MatISGetLocalMat. > > > > > > > > > > > > > > > > > > > > > > > > > > > > How to get assembled local matrix (initial local matrix + > > > > > > > neigbhor > > > > > > > contributions on the borders) ? (expected result is diag: 2, 1) > > > > > > > > > > > > > > > > > > > > > > > > > > > You can always use > > > > > > > > > > > > > > > > > > > > > http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatGetSubMatrix.html > > > > > > > > > > > > > > > > > > > > > http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatGetSubMatrices.html > > > > > > > > > > > > > > > > > > > > > to get copies, but if you just want to build things, you can use > > > > > > > > > > > > > > > > > > > > > http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatGetLocalSubMatrix.html > > > > > > > > > > > > > > > > > > > > > Thanks, > > > > > > > > > > > > > > > > > > > > > Matt > > > > > > > > > > > > > > > > > > > > > > Franck > > > > > > > > > > > > > > > > > > > > > > > > > > > -- > > > > > > > > > > > > > > > > > > > > > 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 > > > > > > > > > > > > > > > > > > > > > http://www.caam.rice.edu/~mk51/ > > > > > > > > > > > > > > > > > <matISLocalMat.cpp> >
