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

Tue May 23 13:23:49 CDT 2017

> 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
> 
> 
> 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 <mailto: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.
> 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 ?
> 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.
> 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 <mailto:stefano.zampini at gmail.com>>
> À: "petsc-maint" <knepley at gmail.com <mailto:knepley at gmail.com>>
> Cc: "petsc-dev" <petsc-dev at mcs.anl.gov <mailto:petsc-dev at mcs.anl.gov>>, "PETSc users list" <petsc-users at mcs.anl.gov <mailto:petsc-users at mcs.anl.gov>>, "Franck Houssen" <franck.houssen at inria.fr <mailto: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 <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 <mailto:knepley at gmail.com>> ha scritto:
> On Sun, May 21, 2017 at 11:23 AM, Franck Houssen <franck.houssen at inria.fr <mailto: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/MatGetSubMatrix.html>
> http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatGetSubMatrices.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 <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/ <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/20170523/2b2916aa/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: matISLocalMat.cpp
Type: application/octet-stream
Size: 3788 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20170523/2b2916aa/attachment.obj>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20170523/2b2916aa/attachment-0001.html>