[petsc-users] Understanding index sets for PCGASM

[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.

Thanks in advance,
Leonardo