1: program main
  2: #include <petsc/finclude/petscdmda.h>
  3: #include <petsc/finclude/petscksp.h>
  4:   use petscdm
  5:   use petscdmda
  6:   use petscksp
  7:   implicit none

  9:   PetscInt is, js, iw, jw
 10:   PetscInt one, three
 11:   PetscErrorCode ierr
 12:   KSP ksp
 13:   DM dm
 14:   external ComputeRHS, ComputeMatrix, ComputeInitialGuess

 16:   one = 1
 17:   three = 3

 19:   PetscCallA(PetscInitialize(ierr))
 20:   PetscCallA(KSPCreate(PETSC_COMM_WORLD, ksp, ierr))
 21:   PetscCallA(DMDACreate2D(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, three, three, PETSC_DECIDE, PETSC_DECIDE, one, one, PETSC_NULL_INTEGER_ARRAY, PETSC_NULL_INTEGER_ARRAY, dm, ierr))
 22:   PetscCallA(DMSetFromOptions(dm, ierr))
 23:   PetscCallA(DMSetUp(dm, ierr))
 24:   PetscCallA(KSPSetDM(ksp, dm, ierr))
 25:   PetscCallA(KSPSetComputeInitialGuess(ksp, ComputeInitialGuess, 0, ierr))
 26:   PetscCallA(KSPSetComputeRHS(ksp, ComputeRHS, 0, ierr))
 27:   PetscCallA(KSPSetComputeOperators(ksp, ComputeMatrix, 0, ierr))
 28:   PetscCallA(DMDAGetCorners(dm, is, js, PETSC_NULL_INTEGER, iw, jw, PETSC_NULL_INTEGER, ierr))
 29:   PetscCallA(KSPSetFromOptions(ksp, ierr))
 30:   PetscCallA(KSPSetUp(ksp, ierr))
 31:   PetscCallA(KSPSolve(ksp, PETSC_NULL_VEC, PETSC_NULL_VEC, ierr))
 32:   PetscCallA(KSPDestroy(ksp, ierr))
 33:   PetscCallA(DMDestroy(dm, ierr))
 34:   PetscCallA(PetscFinalize(ierr))
 35: end

 37: subroutine ComputeInitialGuess(ksp, b, ctx, ierr)
 38:   use petsckspdef
 39:   use petscvec
 40:   implicit none
 41:   PetscErrorCode ierr
 42:   KSP ksp
 43:   PetscInt ctx(*)
 44:   Vec b
 45:   PetscScalar h

 47:   h = 0.0
 48:   PetscCall(VecSet(b, h, ierr))
 49: end subroutine

 51: subroutine ComputeRHS(ksp, b, dummy, ierr)
 52:   use petscksp
 53:   use petscdmda
 54:   use petscvec
 55:   implicit none

 57:   PetscErrorCode ierr
 58:   KSP ksp
 59:   Vec b
 60:   integer dummy(*)
 61:   PetscScalar h, Hx, Hy
 62:   PetscInt mx, my
 63:   DM dm

 65:   PetscCall(KSPGetDM(ksp, dm, ierr))
 66:   PetscCall(DMDAGetInfo(dm, PETSC_NULL_INTEGER, mx, my, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_DMBOUNDARYTYPE, PETSC_NULL_DMBOUNDARYTYPE, PETSC_NULL_DMBOUNDARYTYPE, PETSC_NULL_DMDASTENCILTYPE, ierr))

 68:   Hx = 1.0/real(mx - 1)
 69:   Hy = 1.0/real(my - 1)
 70:   h = Hx*Hy
 71:   PetscCall(VecSet(b, h, ierr))
 72: end subroutine

 74: subroutine ComputeMatrix(ksp, A, B, dummy, ierr)
 75:   use petscmat
 76:   use petscksp
 77:   use petscdmda
 78:   implicit none
 79:   PetscErrorCode ierr
 80:   KSP ksp
 81:   Mat A, B
 82:   integer dummy(*)
 83:   DM dm

 85:   PetscInt i, j, mx, my, xm
 86:   PetscInt ym, xs, ys, i1, i5
 87:   PetscScalar v(5), Hx, Hy
 88:   PetscScalar HxdHy, HydHx
 89:   MatStencil row(1), col(5)

 91:   i1 = 1
 92:   i5 = 5
 93:   PetscCall(KSPGetDM(ksp, dm, ierr))
 94:   PetscCall(DMDAGetInfo(dm, PETSC_NULL_INTEGER, mx, my, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_DMBOUNDARYTYPE, PETSC_NULL_DMBOUNDARYTYPE, PETSC_NULL_DMBOUNDARYTYPE, PETSC_NULL_DMDASTENCILTYPE, ierr))

 96:   Hx = 1.0/real(mx - 1)
 97:   Hy = 1.0/real(my - 1)
 98:   HxdHy = Hx/Hy
 99:   HydHx = Hy/Hx
100:   PetscCall(DMDAGetCorners(dm, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr))
101:   do 10, j = ys, ys + ym - 1
102:     do 20, i = xs, xs + xm - 1
103:       row(1)%i = i
104:       row(1)%j = j
105:       if (i == 0 .or. j == 0 .or. i == mx - 1 .or. j == my - 1) then
106:         v(1) = 2.0*(HxdHy + HydHx)
107:         PetscCall(MatSetValuesStencil(B, i1, row, i1, row, v, INSERT_VALUES, ierr))
108:       else
109:         v(1) = -HxdHy
110:         col(1)%i = i
111:         col(1)%j = j - 1
112:         v(2) = -HydHx
113:         col(2)%i = i - 1
114:         col(2)%j = j
115:         v(3) = 2.0*(HxdHy + HydHx)
116:         col(3)%i = i
117:         col(3)%j = j
118:         v(4) = -HydHx
119:         col(4)%i = i + 1
120:         col(4)%j = j
121:         v(5) = -HxdHy
122:         col(5)%i = i
123:         col(5)%j = j + 1
124:         PetscCall(MatSetValuesStencil(B, i1, row, i5, col, v, INSERT_VALUES, ierr))
125:       end if
126: 20    continue
127: 10    continue
128:       PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY, ierr))
129:       PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY, ierr))
130:       if (A /= B) then
131:         PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
132:         PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))
133:       end if
134:       end subroutine

136: !/*TEST
137: !
138: !   test:
139: !      nsize: 4
140: !      args: -ksp_monitor_short -da_refine 5 -pc_type mg -pc_mg_levels 5 -mg_levels_ksp_type chebyshev -mg_levels_ksp_max_it 2 -mg_levels_pc_type jacobi -ksp_pc_side right
141: !
142: !TEST*/