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

Tue May 23 11:34:53 CDT 2017

OK. I am supposed to destroy the matrix returned by MatISGetMPIXAIJ ? 
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 ? 

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/
> > > 
> > 
> 
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20170523/68af16d0/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: matISLocalMat.cpp
Type: text/x-c++src
Size: 3685 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20170523/68af16d0/attachment.cpp>