1: #include <petsc/finclude/petscdmda.h>
  2: #include <petsc/finclude/petscksp.h>
  3: module ex45fmodule
  4:   use petscmat
  5:   use petscksp
  6:   use petscdmda
  7:   use petscvec
  8:   implicit none

 10: contains

 12:   subroutine ComputeInitialGuess(ksp, b, ctx, ierr)
 13:     PetscErrorCode ierr
 14:     KSP ksp
 15:     PetscInt ctx(*)
 16:     Vec b
 17:     PetscScalar h

 19:     h = 0.0
 20:     PetscCall(VecSet(b, h, ierr))
 21:   end subroutine

 23:   subroutine ComputeRHS(ksp, b, dummy, ierr)

 25:     PetscErrorCode ierr
 26:     KSP ksp
 27:     Vec b
 28:     integer dummy(*)
 29:     PetscScalar h, Hx, Hy
 30:     PetscInt mx, my
 31:     DM dm

 33:     PetscCall(KSPGetDM(ksp, dm, ierr))
 34:     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))

 36:     Hx = 1.0/real(mx - 1)
 37:     Hy = 1.0/real(my - 1)
 38:     h = Hx*Hy
 39:     PetscCall(VecSet(b, h, ierr))
 40:   end subroutine

 42:   subroutine ComputeMatrix(ksp, A, B, dummy, ierr)
 43:     PetscErrorCode ierr
 44:     KSP ksp
 45:     Mat A, B
 46:     integer dummy(*)
 47:     DM dm

 49:     PetscInt i, j, mx, my, xm
 50:     PetscInt ym, xs, ys, i1, i5
 51:     PetscScalar v(5), Hx, Hy
 52:     PetscScalar HxdHy, HydHx
 53:     MatStencil row(1), col(5)

 55:     i1 = 1
 56:     i5 = 5
 57:     PetscCall(KSPGetDM(ksp, dm, ierr))
 58:     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))

 60:     Hx = 1.0/real(mx - 1)
 61:     Hy = 1.0/real(my - 1)
 62:     HxdHy = Hx/Hy
 63:     HydHx = Hy/Hx
 64:     PetscCall(DMDAGetCorners(dm, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr))
 65:     do j = ys, ys + ym - 1
 66:       do i = xs, xs + xm - 1
 67:         row(1)%i = i
 68:         row(1)%j = j
 69:         if (i == 0 .or. j == 0 .or. i == mx - 1 .or. j == my - 1) then
 70:           v(1) = 2.0*(HxdHy + HydHx)
 71:           PetscCall(MatSetValuesStencil(B, i1, row, i1, row, v, INSERT_VALUES, ierr))
 72:         else
 73:           v(1) = -HxdHy
 74:           col(1)%i = i
 75:           col(1)%j = j - 1
 76:           v(2) = -HydHx
 77:           col(2)%i = i - 1
 78:           col(2)%j = j
 79:           v(3) = 2.0*(HxdHy + HydHx)
 80:           col(3)%i = i
 81:           col(3)%j = j
 82:           v(4) = -HydHx
 83:           col(4)%i = i + 1
 84:           col(4)%j = j
 85:           v(5) = -HxdHy
 86:           col(5)%i = i
 87:           col(5)%j = j + 1
 88:           PetscCall(MatSetValuesStencil(B, i1, row, i5, col, v, INSERT_VALUES, ierr))
 89:         end if
 90:       end do
 91:     end do
 92:     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY, ierr))
 93:     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY, ierr))
 94:     if (A /= B) then
 95:       PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
 96:       PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))
 97:     end if
 98:   end subroutine
 99: end module ex45fmodule

101: program main
102:   use petscdmda
103:   use petscksp
104:   use ex45fmodule
105:   implicit none

107:   PetscInt is, js, iw, jw
108:   PetscInt one, three
109:   PetscErrorCode ierr
110:   KSP ksp
111:   DM dm

113:   one = 1
114:   three = 3

116:   PetscCallA(PetscInitialize(ierr))
117:   PetscCallA(KSPCreate(PETSC_COMM_WORLD, ksp, ierr))
118:   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))
119:   PetscCallA(DMSetFromOptions(dm, ierr))
120:   PetscCallA(DMSetUp(dm, ierr))
121:   PetscCallA(KSPSetDM(ksp, dm, ierr))
122:   PetscCallA(KSPSetComputeInitialGuess(ksp, ComputeInitialGuess, 0, ierr))
123:   PetscCallA(KSPSetComputeRHS(ksp, ComputeRHS, 0, ierr))
124:   PetscCallA(KSPSetComputeOperators(ksp, ComputeMatrix, 0, ierr))
125:   PetscCallA(DMDAGetCorners(dm, is, js, PETSC_NULL_INTEGER, iw, jw, PETSC_NULL_INTEGER, ierr))
126:   PetscCallA(KSPSetFromOptions(ksp, ierr))
127:   PetscCallA(KSPSetUp(ksp, ierr))
128:   PetscCallA(KSPSolve(ksp, PETSC_NULL_VEC, PETSC_NULL_VEC, ierr))
129:   PetscCallA(KSPDestroy(ksp, ierr))
130:   PetscCallA(DMDestroy(dm, ierr))
131:   PetscCallA(PetscFinalize(ierr))
132: end

134: !/*TEST
135: !
136: !   test:
137: !      nsize: 4
138: !      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
139: !
140: !TEST*/