# [petsc-users] Understanding index sets for PCGASM

Fri May 5 10:17:17 CDT 2023

```Thanks a lot.
If this can help, we should need (not much more than) the functionalities
from https://petsc.org/release/src/ksp/ksp/tutorials/ex62.c.html:

- PCGASMSetSubdomains
- PCGASMDestroySubdomains
- PCGASMGetSubKSP

Best,
Leonardo

Il giorno ven 5 mag 2023 alle ore 17:00 Barry Smith <bsmith at petsc.dev> ha
scritto:

>
>    I will add the two interfaces you requested today. (Likely you may need
> more also).
>
>   Barry
>
>
> On May 4, 2023, at 6:01 PM, Matthew Knepley <knepley at gmail.com> wrote:
>
> On Thu, May 4, 2023 at 1:43 PM LEONARDO MUTTI <
>
>> Of course, I'll try to explain.
>>
>> I am solving a parabolic equation with space-time FEM and I want an
>> efficient solver/preconditioner for the resulting system.
>> The corresponding matrix, call it X, has an e.g. block bi-diagonal
>> structure, if the cG(1)-dG(0) method is used (i.e. implicit Euler solved in
>> batch).
>> Every block-row of X corresponds to a time instant.
>>
>> I want to introduce parallelism in time by subdividing X into overlapping
>> submatrices of e.g 2x2 or 3x3 blocks, along the block diagonal.
>> For instance, call X_i the individual blocks. The submatrices would be,
>> for various i, (X_{i-1,i-1},X_{i-1,i};X_{i,i-1},X_{i,i}).
>> I'd like each submatrix to be solved in parallel, to combine the various
>> results together in an ASM like fashion.
>> Every submatrix has thus a predecessor and a successor, and it overlaps
>> with both, so that as far as I could understand, GASM has to be used in
>> place of ASM.
>>
>
> Yes, ordered that way you need GASM. I wonder if inverting the ordering
> would be useful, namely putting the time index on the inside.
> Then the blocks would be over all time, but limited space, which is more
> the spirit of ASM I think.
>
> Have you considered waveform relaxation for this problem?
>
>    Thanks,
>
>      Matt
>
>
>> Hope this helps.
>> Best,
>> Leonardo
>>
>> Il giorno gio 4 mag 2023 alle ore 18:05 Matthew Knepley <
>> knepley at gmail.com> ha scritto:
>>
>>> On Thu, May 4, 2023 at 11:24 AM LEONARDO MUTTI <
>>>
>>>> Thank you for the help.
>>>>
>>>>
>>>> *      call PCGASMSetSubdomains(pc,NSub, subdomains_IS,
>>>> inflated_IS,ierr)      call
>>>> PCGASMDestroySubdomains(NSub,subdomains_IS,inflated_IS,ierr)*
>>>> results in:
>>>>
>>>> *      Error LNK2019 unresolved external symbol PCGASMDESTROYSUBDOMAINS
>>>> referenced in function ...  *
>>>>
>>>> *      Error LNK2019 unresolved external symbol PCGASMSETSUBDOMAINS
>>>> referenced in function ...  *
>>>> I'm not sure if the interfaces are missing or if I have a compilation
>>>> problem.
>>>>
>>>
>>> I just want to make sure you really want GASM. It sounded like you might
>>> able to do what you want just with ASM.
>>> Can you tell me again what you want to do overall?
>>>
>>>   Thanks,
>>>
>>>     Matt
>>>
>>>
>>>> Thank you again.
>>>> Best,
>>>> Leonardo
>>>>
>>>> Il giorno sab 29 apr 2023 alle ore 20:30 Barry Smith <bsmith at petsc.dev>
>>>> ha scritto:
>>>>
>>>>>
>>>>>    Thank you for the test code. I have a fix in the branch
>>>>> barry/2023-04-29/fix-pcasmcreatesubdomains2d
>>>>> <https://gitlab.com/petsc/petsc/-/commits/barry/2023-04-29/fix-pcasmcreatesubdomains2d> with
>>>>> merge request https://gitlab.com/petsc/petsc/-/merge_requests/6394
>>>>>
>>>>>    The functions did not have proper Fortran stubs and interfaces so I
>>>>> had to provide them manually in the new branch.
>>>>>
>>>>>    Use
>>>>>
>>>>>    git fetch
>>>>>    git checkout barry/2023-04-29/fix-pcasmcreatesubdomains2d
>>>>> <https://gitlab.com/petsc/petsc/-/commits/barry/2023-04-29/fix-pcasmcreatesubdomains2d>
>>>>>    ./configure etc
>>>>>
>>>>> to change things slightly and I updated the error handling for the latest
>>>>> version.
>>>>>
>>>>>    Please let us know if you have any later questions.
>>>>>
>>>>>   Barry
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> On Apr 28, 2023, at 12:07 PM, LEONARDO MUTTI <
>>>>>
>>>>> Hello. I am having a hard time understanding the index sets to feed
>>>>> PCGASMSetSubdomains, and I am working in Fortran (as a PETSc novice). To
>>>>> get more intuition on how the IS objects behave I tried the following
>>>>> minimal (non) working example, which should tile a 16x16 matrix into 16
>>>>> square, non-overlapping submatrices:
>>>>>
>>>>> #include <petsc/finclude/petscmat.h>
>>>>> #include <petsc/finclude/petscksp.h>
>>>>> #include <petsc/finclude/petscpc.h>
>>>>>       USE petscmat
>>>>>       USE petscksp
>>>>>       USE petscpc
>>>>>
>>>>>       Mat   :: A
>>>>>       PetscInt :: M, NSubx, dof, overlap, NSub
>>>>>       INTEGER :: I,J
>>>>>       PetscErrorCode      :: ierr
>>>>>       PetscScalar :: v
>>>>>       KSP            :: ksp
>>>>>       PC             :: pc
>>>>>       IS :: subdomains_IS, inflated_IS
>>>>>
>>>>>       call PetscInitialize(PETSC_NULL_CHARACTER , ierr)
>>>>>
>>>>> !-----Create a dummy matrix
>>>>>       M = 16
>>>>>       call MatCreateAIJ(MPI_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE,
>>>>>      & M, M,
>>>>>      & PETSC_DEFAULT_INTEGER, PETSC_NULL_INTEGER,
>>>>>      & PETSC_DEFAULT_INTEGER, PETSC_NULL_INTEGER,
>>>>>      & A, ierr)
>>>>>
>>>>>       DO I=1,M
>>>>>          DO J=1,M
>>>>>             v = I*J
>>>>>             CALL MatSetValue (A,I-1,J-1,v,
>>>>>      &                     INSERT_VALUES , ierr)
>>>>>          END DO
>>>>>       END DO
>>>>>
>>>>>       call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY , ierr)
>>>>>       call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY , ierr)
>>>>>
>>>>> !-----Create KSP and PC
>>>>>       call KSPCreate(PETSC_COMM_WORLD,ksp, ierr)
>>>>>       call KSPSetOperators(ksp,A,A, ierr)
>>>>>       call KSPSetType(ksp,"bcgs",ierr)
>>>>>       call KSPGetPC(ksp,pc,ierr)
>>>>>       call KSPSetUp(ksp, ierr)
>>>>>       call PCSetType(pc,PCGASM, ierr)
>>>>>       call PCSetUp(pc , ierr)
>>>>>
>>>>> !-----GASM setup
>>>>>       NSubx = 4
>>>>>       dof = 1
>>>>>       overlap = 0
>>>>>
>>>>>       call PCGASMCreateSubdomains2D(pc,
>>>>>      &      M, M,
>>>>>      &      NSubx, NSubx,
>>>>>      &      dof, overlap,
>>>>>      &      NSub, subdomains_IS, inflated_IS, ierr)
>>>>>
>>>>>       call ISView(subdomains_IS, PETSC_VIEWER_STDOUT_WORLD, ierr)
>>>>>
>>>>>       call KSPDestroy(ksp, ierr)
>>>>>       call PetscFinalize(ierr)
>>>>>
>>>>> Running this on one processor, I get NSub = 4.
>>>>> If PCASM and PCASMCreateSubdomains2D are used instead, I get NSub = 16
>>>>> as expected.
>>>>> Moreover, I get in the end "forrtl: severe (157): Program Exception -
>>>>> access violation". So:
>>>>> 1) why do I get two different results with ASM, and GASM?
>>>>> 2) why do I get access violation and how can I solve this?
>>>>> In fact, in C, subdomains_IS, inflated_IS should pointers to IS
>>>>> objects. As I see on the Fortran interface, the arguments to
>>>>> PCGASMCreateSubdomains2D are IS objects:
>>>>>
>>>>>        subroutine PCGASMCreateSubdomains2D(a,b,c,d,e,f,g,h,i,j,z)
>>>>>        import tPC,tIS
>>>>>        PC a ! PC
>>>>>        PetscInt b ! PetscInt
>>>>>        PetscInt c ! PetscInt
>>>>>        PetscInt d ! PetscInt
>>>>>        PetscInt e ! PetscInt
>>>>>        PetscInt f ! PetscInt
>>>>>        PetscInt g ! PetscInt
>>>>>        PetscInt h ! PetscInt
>>>>>        IS i ! IS
>>>>>        IS j ! IS
>>>>>        PetscErrorCode z
>>>>>        end subroutine PCGASMCreateSubdomains2D
>>>>> Thus:
>>>>> 3) what should be inside e.g., subdomains_IS? I expect it to contain,
>>>>> for every created subdomain, the list of rows and columns defining the
>>>>> subblock in the matrix, am I right?
>>>>>
>>>>> Context: I have a block-tridiagonal system arising from space-time
>>>>> finite elements, and I want to solve it with GMRES+PCGASM preconditioner,
>>>>> where each overlapping submatrix is on the diagonal and of size 3x3 blocks
>>>>> (and spanning multiple processes). This is PETSc 3.17.1 on Windows.
>>>>>
>>>>> Leonardo
>>>>>
>>>>>
>>>>>
>>>
>>> --
>>> What most experimenters take for granted before they begin their
>>> experiments is infinitely more interesting than any results to which their
>>> -- Norbert Wiener
>>>
>>> https://www.cse.buffalo.edu/~knepley/
>>> <http://www.cse.buffalo.edu/~knepley/>
>>>
>>
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which their
> -- Norbert Wiener
>
> https://www.cse.buffalo.edu/~knepley/
> <http://www.cse.buffalo.edu/~knepley/>
>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20230505/bc8b4fae/attachment.html>
```