[petsc-users] Understanding index sets for PCGASM

Barry Smith bsmith at petsc.dev
Fri May 5 10:00:53 CDT 2023 

   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 <leonardo.mutti01 at universitadipavia.it <mailto:leonardo.mutti01 at universitadipavia.it>> wrote:
>> 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 <mailto:knepley at gmail.com>> ha scritto:
>>> On Thu, May 4, 2023 at 11:24 AM LEONARDO MUTTI <leonardo.mutti01 at universitadipavia.it <mailto: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 <mailto: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 <mailto: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/>
> 
> 
> -- 
> 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/20230505/abf26712/attachment.html>

More information about the petsc-users mailing list