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

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

 14:        one = 1
 15:        three = 3

 17:        PetscCallA(PetscInitialize(ierr))
 18:        PetscCallA(KSPCreate(PETSC_COMM_WORLD,ksp,ierr))
 19:        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))
 20:        PetscCallA(DMSetFromOptions(dm,ierr))
 21:        PetscCallA(DMSetUp(dm,ierr))
 22:        PetscCallA(KSPSetDM(ksp,dm,ierr))
 23:        PetscCallA(KSPSetComputeInitialGuess(ksp,ComputeInitialGuess,0,ierr))
 24:        PetscCallA(KSPSetComputeRHS(ksp,ComputeRHS,0,ierr))
 25:        PetscCallA(KSPSetComputeOperators(ksp,ComputeMatrix,0,ierr))
 26:        PetscCallA(DMDAGetCorners(dm,is,js,PETSC_NULL_INTEGER,iw,jw,PETSC_NULL_INTEGER,ierr))
 27:        PetscCallA(KSPSetFromOptions(ksp,ierr))
 28:        PetscCallA(KSPSetUp(ksp,ierr))
 29:        PetscCallA(KSPSolve(ksp,PETSC_NULL_VEC,PETSC_NULL_VEC,ierr))
 30:        PetscCallA(KSPDestroy(ksp,ierr))
 31:        PetscCallA(DMDestroy(dm,ierr))
 32:        PetscCallA(PetscFinalize(ierr))
 33:        end

 35:        subroutine ComputeInitialGuess(ksp,b,ctx,ierr)
 36:        use petscksp
 37:        implicit none
 38:        PetscErrorCode  ierr
 39:        KSP ksp
 40:        PetscInt ctx(*)
 41:        Vec b
 42:        PetscScalar  h

 44:        h=0.0
 45:        PetscCall(VecSet(b,h,ierr))
 46:        end subroutine

 48:        subroutine ComputeRHS(ksp,b,dummy,ierr)
 49:        use petscksp
 50:        implicit none

 52:        PetscErrorCode  ierr
 53:        KSP ksp
 54:        Vec b
 55:        integer dummy(*)
 56:        PetscScalar  h,Hx,Hy
 57:        PetscInt  mx,my
 58:        DM dm

 60:        PetscCall(KSPGetDM(ksp,dm,ierr))
 61:        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_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))

 63:        Hx = 1.0 / real(mx-1)
 64:        Hy = 1.0 / real(my-1)
 65:        h = Hx*Hy
 66:        PetscCall(VecSet(b,h,ierr))
 67:        end subroutine

 69:       subroutine ComputeMatrix(ksp,A,B,dummy,ierr)
 70:       use petscksp
 71:        implicit none
 72:        PetscErrorCode  ierr
 73:        KSP ksp
 74:        Mat A,B
 75:        integer dummy(*)
 76:        DM dm

 78:       PetscInt    i,j,mx,my,xm
 79:       PetscInt    ym,xs,ys,i1,i5
 80:       PetscScalar  v(5),Hx,Hy
 81:       PetscScalar  HxdHy,HydHx
 82:       MatStencil   row(4),col(4,5)

 84:       i1 = 1
 85:       i5 = 5
 86:       PetscCall(KSPGetDM(ksp,dm,ierr))
 87:       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_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))

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

129: !/*TEST
130: !
131: !   test:
132: !      nsize: 4
133: !      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
134: !
135: !TEST*/