[petsc-users] [petsc-dev] Using PETSc MatIS, how to get local matrix (= one domain) before and after assembly ?

Wed May 24 08:11:28 CDT 2017

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" <stefano.zampini at gmail.com>
> À: "Franck Houssen" <franck.houssen at inria.fr>
> Cc: "petsc-dev" <petsc-dev at mcs.anl.gov>, "PETSc users list"
> <petsc-users at mcs.anl.gov>, "petsc-maint" <knepley at gmail.com>
> 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 < franck.houssen at inria.fr >
> > 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" < stefano.zampini at gmail.com >
> > 
> 
> > > À: "Franck Houssen" < franck.houssen at inria.fr >
> > 
> 
> > > Cc: "petsc-dev" < petsc-dev at mcs.anl.gov >, "PETSc users list" <
> > > petsc-users at mcs.anl.gov >, "petsc-maint" < knepley at gmail.com >
> > 
> 
> > > 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 < franck.houssen at inria.fr >
> > > > 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" < stefano.zampini at gmail.com >
> > > > 
> > > 
> > 
> 
> > > > > À: "Franck Houssen" < franck.houssen at inria.fr >
> > > > 
> > > 
> > 
> 
> > > > > Cc: "petsc-dev" < petsc-dev at mcs.anl.gov >, "PETSc users list" <
> > > > > petsc-users at mcs.anl.gov >, "petsc-maint" < knepley at gmail.com >
> > > > 
> > > 
> > 
> 
> > > > > 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" < franck.houssen at inria.fr > 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" < stefano.zampini at gmail.com >
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> 
> > > > > > > À: "petsc-maint" < knepley at gmail.com >
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> 
> > > > > > > Cc: "petsc-dev" < petsc-dev at mcs.anl.gov >, "PETSc users list" <
> > > > > > > petsc-users at mcs.anl.gov >, "Franck Houssen" <
> > > > > > > franck.houssen at inria.fr
> > > > > > > >
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> 
> > > > > > > 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" < knepley at gmail.com > ha
> > > > > > > scritto:
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> 

> > > > > > > > On Sun, May 21, 2017 at 11:23 AM, Franck Houssen <
> > > > > > > > franck.houssen at inria.fr
> > > > > > > > >
> > > > > > > > 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>
> > > 
> > 
> 
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20170524/22a01549/attachment-0001.html>