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 :: M, M2, NSubx, dof, overlap, NSub, i1
  9:   PetscInt :: 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:   M = 16
 20:   M2 = M*M
 21:   i1 = 1
 22:   PetscCallA(MatCreateFromOptions(PETSC_COMM_WORLD, PETSC_NULL_CHARACTER, i1, PETSC_DECIDE, PETSC_DECIDE, M2, M2, A, ierr))
 23:   do I = 1, M2
 24:     do J = 1, M2
 25:       v = I*J
 26:       PetscCallA(MatSetValue(A, I - 1, J - 1, v, INSERT_VALUES, ierr))
 27:     end do
 28:   end do

 30:   PetscCallA(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
 31:   PetscCallA(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))
 32:   PetscCallA(PCCreate(PETSC_COMM_WORLD, pc, ierr))
 33:   PetscCallA(PCSetOperators(pc, A, A, ierr))
 34:   PetscCallA(PCSetType(pc, PCGASM, ierr))

 36:   NSubx = 4
 37:   dof = 1
 38:   overlap = 0

 40:   PetscCallA(PCGASMCreateSubdomains2D(pc, M, M, NSubx, NSubx, dof, overlap, NSub, subdomains_IS, inflated_IS, ierr))
 41:   PetscCallA(PCGASMSetSubdomains(pc, NSub, subdomains_IS, inflated_IS, ierr))
 42:   PetscCallA(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, singleton, ierr))
 43:   PetscCallA(PetscViewerASCIIPrintf(singleton, 'GASM index sets from this MPI process\n', ierr))
 44:   do i = 1, Nsub
 45:     PetscCallA(ISView(subdomains_IS(i), singleton, ierr))
 46:   end do
 47:   PetscCallA(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, singleton, ierr))
 48:   PetscCallA(PCGASMDestroySubdomains(NSub, subdomains_IS, inflated_IS, ierr))

 50:   if (size == 1) then
 51:     ! this routine only works on one rank
 52:     PetscCallA(PCASMCreateSubdomains2D(M, M, NSubx, NSubx, dof, overlap, NSub, subdomains_IS, inflated_IS, ierr))
 53:     do i = 1, Nsub
 54:       PetscCallA(ISView(subdomains_IS(i), PETSC_VIEWER_STDOUT_SELF, ierr))
 55:     end do
 56:     PetscCallA(PCASMDestroySubdomains(NSub, subdomains_IS, inflated_IS, ierr))
 57:   end if

 59:   PetscCallA(MatDestroy(A, ierr))
 60:   PetscCallA(PCDestroy(pc, ierr))
 61:   PetscCallA(PetscFinalize(ierr))
 62: end

 64: !/*TEST
 65: !
 66: !   test:
 67: !     suffix: 1
 68: !
 69: !   test:
 70: !     suffix: 2
 71: !     nsize: 2
 72: !
 73: !TEST*/