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

[Stefano Zampini]

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/353a6272/attachment-0001.html>