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

[Franck Houssen]

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 

----- Mail original -----

> 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/94f74525/attachment.html>