1:       program main
  2:       implicit none

  4: #include <petsc/finclude/petscsys.h>
  5: #include <petsc/finclude/petscvec.h>
  6: #include <petsc/finclude/petscmat.h>
  7: #include <petsc/finclude/petscpc.h>
  8: #include <petsc/finclude/petscksp.h>
  9: #include <petsc/finclude/petscdm.h>
 10: #include <petsc/finclude/petscdmda.h>

 12:        PetscInt is,js,iw,jw
 13:        PetscInt one,three
 14:        PetscErrorCode ierr
 15:        KSP ksp
 16:        DM dm
 17:        external ComputeRHS,ComputeMatrix,ComputeInitialGuess

 19:        one = 1
 20:        three = 3

 22:        call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 23:        call KSPCreate(MPI_COMM_WORLD,ksp,ierr)
 24:        call DMDACreate2D(MPI_COMM_WORLD, DM_BOUNDARY_NONE,              &
 25:      &    DM_BOUNDARY_NONE, DMDA_STENCIL_STAR,-three,-three,            &
 26:      &    PETSC_DECIDE,PETSC_DECIDE,one,one, PETSC_NULL_INTEGER,        &
 27:      &    PETSC_NULL_INTEGER, dm, ierr)
 28:        call KSPSetDM(ksp,dm,ierr)
 29:        call KSPSetComputeInitialGuess(ksp,ComputeInitialGuess,             &
 30:      &                                PETSC_NULL_OBJECT,ierr)
 31:        call KSPSetComputeRHS(ksp,ComputeRHS,PETSC_NULL_OBJECT,ierr)
 32:        call KSPSetComputeOperators(ksp,ComputeMatrix,                   &
 33:      &      PETSC_NULL_OBJECT,ierr)
 34:        call DMDAGetCorners(dm,is,js,PETSC_NULL_INTEGER,iw,jw,             &
 35:      &                     PETSC_NULL_INTEGER,ierr)
 36:        call KSPSetFromOptions(ksp,ierr)
 37:        call KSPSetUp(ksp,ierr)
 38:        call KSPSolve(ksp,PETSC_NULL_OBJECT,PETSC_NULL_OBJECT,ierr)
 39:        call KSPDestroy(ksp,ierr)
 40:        call DMDestroy(dm,ierr)
 41:        call PetscFinalize(ierr)
 42:        end


 45:        subroutine ComputeInitialGuess(ksp,b,ctx,ierr)
 46:        implicit none
 47:        PetscErrorCode  ierr
 48:        KSP ksp
 49:        PetscInt ctx(*)
 50:        Vec b
 51:        PetscScalar  h

 53:        h=0.0
 54:        call VecSet(b,h,ierr)
 55:        end subroutine

 57:        subroutine ComputeRHS(ksp,b,dummy,ierr)
 58:        implicit none
 59: #include <petsc/finclude/petscsys.h>
 60: #include <petsc/finclude/petscvec.h>
 61: #include <petsc/finclude/petscmat.h>
 62:        PetscErrorCode  ierr
 63:        KSP ksp
 64:        Vec b
 65:        integer dummy(*)
 66:        PetscScalar  h,Hx,Hy
 67:        PetscInt  mx,my
 68:        DM dm

 70:        call KSPGetDM(ksp,dm,ierr)
 71:        call DMDAGetInfo(dm,PETSC_NULL_INTEGER,mx,my,PETSC_NULL_INTEGER,  &
 72:      &                     PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,             &
 73:      &                     PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,             &
 74:      &                     PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,             &
 75:      &                     PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,             &
 76:      &                     PETSC_NULL_INTEGER,ierr)

 78:        Hx = 1.0 / real(mx-1)
 79:        Hy = 1.0 / real(my-1)
 80:        h = Hx*Hy
 81:        call VecSet(b,h,ierr)
 82:        end subroutine

 84:        subroutine ComputeMatrix(ksp,A,B,dummy,ierr)
 85:        implicit none
 86: #include <petsc/finclude/petscsys.h>
 87: #include <petsc/finclude/petscvec.h>
 88: #include <petsc/finclude/petscmat.h>
 89:        PetscErrorCode  ierr
 90:        KSP ksp
 91:        Mat A,B
 92:        integer dummy(*)
 93:        DM dm

 95:       PetscInt    i,j,mx,my,xm
 96:       PetscInt    ym,xs,ys,i1,i5
 97:       PetscScalar  v(5),Hx,Hy
 98:       PetscScalar  HxdHy,HydHx
 99:       MatStencil   row(4),col(4,5)

101:       i1 = 1
102:       i5 = 5
103:       call KSPGetDM(ksp,dm,ierr)
104:       call DMDAGetInfo(dm,PETSC_NULL_INTEGER,mx,my,PETSC_NULL_INTEGER,  &
105:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,             &
106:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,             &
107:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,             &
108:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,             &
109:      &               PETSC_NULL_INTEGER,ierr)

111:       Hx = 1.0 / real(mx-1)
112:       Hy = 1.0 / real(my-1)
113:       HxdHy = Hx/Hy
114:       HydHx = Hy/Hx
115:       call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,            &
116:      &     PETSC_NULL_INTEGER,ierr)
117:       do 10,j=ys,ys+ym-1
118:         do 20,i=xs,xs+xm-1
119:           row(MatStencil_i) = i
120:           row(MatStencil_j) = j
121:           if (i.eq.0 .or. j.eq.0 .or. i.eq.mx-1 .or. j.eq.my-1 ) then
122:             v(1) = 2.0*(HxdHy + HydHx)
123:             call MatSetValuesStencil(B,i1,row,i1,row,v,                 &
124:      &           INSERT_VALUES,ierr)
125:           else
126:             v(1) = -HxdHy
127:             col(MatStencil_i,1) = i
128:             col(MatStencil_j,1) = j-1
129:             v(2) = -HydHx
130:             col(MatStencil_i,2) = i-1
131:             col(MatStencil_j,2) = j
132:             v(3) = 2.0*(HxdHy + HydHx)
133:             col(MatStencil_i,3) = i
134:             col(MatStencil_j,3) = j
135:             v(4) = -HydHx
136:             col(MatStencil_i,4) = i+1
137:             col(MatStencil_j,4) = j
138:             v(5) = -HxdHy
139:             col(MatStencil_i,5) = i
140:             col(MatStencil_j,5) = j+1
141:             call MatSetValuesStencil(B,i1,row,i5,col,v,                 &
142:      &           INSERT_VALUES,ierr)
143:             endif
144:  20      continue
145:  10   continue
146:        call MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr)
147:        call MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr)
148:        if ( A .ne. B) then
149:          call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
150:          call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
151:        endif
152:        end subroutine