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" <[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/
> > >
> >
>
// Using PETSc MatIS, how to get local matrix (= one domain) before and after assembly (neumann/dirichlet).
//
// A 3x3 global matrix is built (diag: 1, 2, 1): it's made of 2 overlapping 2x2 local matrix (diag: 1, 1).
// Get non assembled local matrix : OK with MatISGetLocalMat.
//
// How to get assembled local matrix (initial local matrix + neigbhor contributions on the borders) ?
// Need to use MatGetLocalSubMatrix ? Not working ?! Expected result is Diag: 2, 1.
//
// ~> g++ -o matISLocalMat.exe matISLocalMat.cpp -lpetsc -lm; mpirun -n 2 matISLocalMat.exe
#include <iostream>
#include "petsc.h"
int main(int argc,char **argv) {
PetscInitialize(&argc, &argv, NULL, NULL);
int size = 0; MPI_Comm_size(MPI_COMM_WORLD, &size); if (size != 2) return 1;
int rank = 0; MPI_Comm_rank(MPI_COMM_WORLD, &rank);
PetscInt localSize = 2, globalSize = 3;
PetscInt localIdx[2] = {0, 0};
if (rank == 0) {localIdx[0] = 0; localIdx[1] = 1;}
else {localIdx[0] = 1; localIdx[1] = 2;}
ISLocalToGlobalMapping map;
ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD, 1, localSize, localIdx, PETSC_COPY_VALUES, &map);
Mat A;
// MatIS means that the matrix is not assembled. The easiest way to think of this is that processes do not have
// to hold full rows. One process can hold part of row i, and another processes can hold another part. However, there are still
// the same number of global rows. The local size here is not the size of the local IS block, since that is a property only of MatIS. It is the
// size of the local piece of the vector you multiply. This allows PETSc to understand the parallel layout of the Vec, and how it
// matched the Mat.
// Conclusion: it's crucial to use PETSC_DECIDE in place of localSize in your call to MatCreateIS.
MatCreateIS(PETSC_COMM_WORLD, 1, PETSC_DECIDE, PETSC_DECIDE, globalSize, globalSize, map, map, &A);
ISLocalToGlobalMappingDestroy(&map);
MatISSetPreallocation(A, localSize, NULL, localSize, NULL);
PetscScalar localVal[4] = {1., 0., 0., 1.};
MatSetValues(A, 2, localIdx, 2, localIdx, localVal, ADD_VALUES);
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
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;
MatISGetLocalMat(A, &nonAssembledLocalMat);
MatView(nonAssembledLocalMat, PETSC_VIEWER_STDOUT_SELF); // Diag: 1, 1
PetscViewerFlush(PETSC_VIEWER_STDOUT_SELF);
}
MPI_Barrier(MPI_COMM_WORLD);
if (rank == 0) std::cout << std::endl << "assembled global matrix:" << std::endl << std::endl;
MPI_Barrier(MPI_COMM_WORLD);
Mat assembledGlobalMat;
MatISGetMPIXAIJ(A, MAT_INITIAL_MATRIX, &assembledGlobalMat); // Collective call.
MatView(assembledGlobalMat, PETSC_VIEWER_STDOUT_WORLD);
PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);
if (rank == 0) std::cout << std::endl << "index set to extract assembled local matrix:" << std::endl << std::endl;
MPI_Barrier(MPI_COMM_WORLD);
IS is; ISCreateGeneral(PETSC_COMM_WORLD, localSize, localIdx, PETSC_COPY_VALUES, &is);
ISView(is, PETSC_VIEWER_STDOUT_WORLD);
PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);
if (rank == 0) std::cout << std::endl << "assembled local matrix:" << std::endl << std::endl;
Mat assembledLocalMat;
MatCreateSubMatrix(assembledGlobalMat, is, is, &assembledLocalMat); // Collective call.
MatView(assembledLocalMat, PETSC_VIEWER_STDOUT_WORLD);
PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);
PetscFinalize();
return 0;
}
