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

Franck Houssen franck.houssen at inria.fr
Wed May 24 04:46:01 CDT 2017 

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" <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-dev/attachments/20170524/56b55fea/attachment.html>

More information about the petsc-dev mailing list