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.eq.0 .or. j.eq.0 .or. i.eq.mx-1 .or. j.eq.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:             endif
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 .ne. B) then
131:          PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr))
132:          PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr))
133:        endif
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*/