[petsc-users] Understanding index sets for PCGASM

Thu May 4 12:43:20 CDT 2023

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.

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 <
> leonardo.mutti01 at universitadipavia.it> wrote:
>
>> Thank you for the help.
>> Adding to my example:
>>
>>
>> *      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
>>>
>>>    Your now working test code is in src/ksp/ksp/tests/ex71f.F90  I had
>>> 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 <
>>> leonardo.mutti01 at universitadipavia.it> wrote:
>>>
>>> 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.
>>>
>>> Thanks in advance,
>>> Leonardo
>>>
>>>
>>>
>
> --
> 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
>
> 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/20230504/047a1187/attachment.html>