1: !     Contributed by leonardo.mutti01@universitadipavia.it
  2: #include <petsc/finclude/petscksp.h>
  3: program main
  4:   use petscksp
  5:   implicit none

  7:   Mat :: A
  8:   PetscInt, parameter :: M = 16, dof = 1, overlap = 0, NSubx = 4
  9:   PetscInt :: NSub, i, j
 10:   PetscMPIInt :: size
 11:   PetscErrorCode :: ierr
 12:   PetscScalar :: v
 13:   PC :: pc
 14:   IS, pointer :: subdomains_IS(:) => null(), inflated_IS(:) => null()
 15:   PetscViewer :: singleton

 17:   PetscCallA(PetscInitialize(PETSC_NULL_CHARACTER, ierr))
 18:   PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, size, ierr))
 19:   PetscCallA(MatCreateFromOptions(PETSC_COMM_WORLD, PETSC_NULL_CHARACTER, 1_PETSC_INT_KIND, PETSC_DECIDE, PETSC_DECIDE, M**2, M**2, A, ierr))
 20:   do i = 1, M**2
 21:     do j = 1, M**2
 22:       v = i*j
 23:       PetscCallA(MatSetValue(A, i - 1, j - 1, v, INSERT_VALUES, ierr))
 24:     end do
 25:   end do

 27:   PetscCallA(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
 28:   PetscCallA(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))
 29:   PetscCallA(PCCreate(PETSC_COMM_WORLD, pc, ierr))
 30:   PetscCallA(PCSetOperators(pc, A, A, ierr))
 31:   PetscCallA(PCSetType(pc, PCGASM, ierr))

 33:   PetscCallA(PCGASMCreateSubdomains2D(pc, M, M, NSubx, NSubx, dof, overlap, NSub, subdomains_IS, inflated_IS, ierr))
 34:   PetscCallA(PCGASMSetSubdomains(pc, NSub, subdomains_IS, inflated_IS, ierr))
 35:   PetscCallA(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, singleton, ierr))
 36:   PetscCallA(PetscViewerASCIIPrintf(singleton, 'GASM index sets from this MPI process\n', ierr))
 37:   do i = 1, Nsub
 38:     PetscCallA(ISView(subdomains_IS(i), singleton, ierr))
 39:   end do
 40:   PetscCallA(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, singleton, ierr))
 41:   PetscCallA(PCGASMDestroySubdomains(NSub, subdomains_IS, inflated_IS, ierr))

 43:   if (size == 1) then
 44:     ! this routine only works on one rank
 45:     PetscCallA(PCASMCreateSubdomains2D(M, M, NSubx, NSubx, dof, overlap, NSub, subdomains_IS, inflated_IS, ierr))
 46:     do i = 1, Nsub
 47:       PetscCallA(ISView(subdomains_IS(i), PETSC_VIEWER_STDOUT_SELF, ierr))
 48:     end do
 49:     PetscCallA(PCASMDestroySubdomains(NSub, subdomains_IS, inflated_IS, ierr))
 50:   end if

 52:   PetscCallA(MatDestroy(A, ierr))
 53:   PetscCallA(PCDestroy(pc, ierr))
 54:   PetscCallA(PetscFinalize(ierr))
 55: end

 57: !/*TEST
 58: !
 59: !   test:
 60: !     suffix: 1
 61: !
 62: !   test:
 63: !     suffix: 2
 64: !     nsize: 2
 65: !
 66: !TEST*/