OK, this is working now ! As the API changed between the latest stable and the master branch, I was actually not using the correct method.
Thanks Stefano, 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é: Mercredi 24 Mai 2017 13:42:10 > Objet: Re: [petsc-dev] Using PETSc MatIS, how to get local matrix (= one > domain) before and after assembly ? > > On May 24, 2017, at 11:46 AM, Franck Houssen < [email protected] > > > wrote: > > > 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. > > MatGetSubMatrix/MatGetSubMatrices have been renamed to > MatCreateSubMatrix/MatCreateSubMatrices in petsc-dev. I thought you were > using the master branch and not the latest release. Sorry for the confusion. > To compile the code I have sent, just rename MatCreateSubMatrices with > MatGetSubMatrices and it should work. > http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatGetSubMatrices.html#MatGetSubMatrices > http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatGetSubMatrix.html#MatGetSubMatrix > > 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). > > In the code, you are already extracting submatrices from MPIAIJ format, not > from MATIS. Attached a code that compiles and runs with petsc-3.7.6 > > 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> > > > > > >
