Actual source code: ex71f.F90

  1: !     Contributed by leonardo.mutti01@universitadipavia.it
  2:       program main
  3: #include <petsc/finclude/petscsys.h>
  4: #include <petsc/finclude/petscmat.h>
  5: #include <petsc/finclude/petscpc.h>
  6:       USE petscmat
  7:       USE petscpc
  8:       implicit none

 10:       Mat :: A
 11:       PetscInt :: M, M2, NSubx, dof, overlap, NSub,i1
 12:       PetscInt :: I, J
 13:       PetscMPIInt :: size
 14:       PetscErrorCode :: ierr
 15:       PetscScalar :: v
 16:       PC :: pc
 17:       IS :: subdomains_IS(20), inflated_IS(20)
 18:       PetscViewer :: singleton

 20:       PetscCallA(PetscInitialize(PETSC_NULL_CHARACTER, ierr))
 21:       PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, size, ierr))
 22:       M = 16
 23:       M2 = M*M
 24:       i1 = 1
 25:       PetscCallA(MatCreateFromOptions(PETSC_COMM_WORLD, PETSC_NULL_CHARACTER,i1, PETSC_DECIDE, PETSC_DECIDE, M2, M2,A, ierr))
 26:       DO I=1,M2
 27:          DO J=1,M2
 28:             v = I*J
 29:             PetscCallA(MatSetValue(A, I-1, J-1, v, INSERT_VALUES, ierr))
 30:          END DO
 31:       END DO

 33:       PetscCallA(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
 34:       PetscCallA(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))
 35:       PetscCallA(PCCreate(PETSC_COMM_WORLD, pc, ierr))
 36:       PetscCallA(PCSetOperators(pc, A, A, ierr))
 37:       PetscCallA(PCSetType(pc, PCGASM, ierr))

 39:       NSubx = 4
 40:       dof = 1
 41:       overlap = 0

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

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

 62:       PetscCallA(MatDestroy(A, ierr))
 63:       PetscCallA(PCDestroy(pc, ierr))
 64:       PetscCallA(PetscFinalize(ierr))
 65:       end

 67: !/*TEST
 68: !
 69: !   test:
 70: !     suffix: 1
 71: !
 72: !   test:
 73: !     suffix: 2
 74: !     nsize: 2
 75: !
 76: !TEST*/