1:       program main
  2:       implicit none

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

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

 18:        call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 19:        call KSPCreate(MPI_COMM_WORLD,ksp,ierr)
 20:        call DMDACreate2D(MPI_COMM_WORLD, DMDA_BOUNDARY_NONE,              &
 21:      &    DMDA_BOUNDARY_NONE, DMDA_STENCIL_STAR, -3,-3, PETSC_DECIDE,       &
 22:      &    PETSC_DECIDE,1,1, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER,       &
 23:      &    dm, ierr)
 24:        call DMSetInitialGuess(dm,ComputeInitialGuess,ierr)
 25:        call KSPSetComputeRHS(ksp,ComputeRHS,PETSC_NULL_OBJECT,ierr)
 26:        call KSPSetComputeOperators(ksp,ComputeMatrix,                   &
 27:      &      PETSC_NULL_OBJECT,ierr)
 28:        call KSPSetDM(ksp,dm,ierr)
 29:        call DMDAGetCorners(dm,is,js,PETSC_NULL_INTEGER,iw,jw,             &
 30:      &                     PETSC_NULL_INTEGER,ierr)
 31:        call KSPSetFromOptions(ksp,ierr)
 32:        call KSPSetUp(ksp,ierr)
 33:        call KSPSolve(ksp,PETSC_NULL_OBJECT,PETSC_NULL_OBJECT,ierr)
 34:        call KSPDestroy(ksp,ierr)
 35:        call DMDestroy(dm,ierr)
 36:        call PetscFinalize(ierr)
 37:        end


 40:        subroutine ComputeInitialGuess(dm,b,ierr)
 41:        implicit none
 42:        PetscErrorCode  ierr
 43:        DM dm
 44:        Vec b
 45:        PetscScalar  h
 46:        h=0.0
 47:        call VecSet(b,h,ierr)
 48:        end subroutine

 50:        subroutine ComputeRHS(ksp,b,dummy,ierr)
 51:        implicit none
 52:        PetscErrorCode  ierr
 53:        KSP ksp
 54:        Vec b
 55:        integer dummy(*)
 56:        PetscScalar  h
 57:        h=1.0
 58:        call VecSet(b,h,ierr)
 59:        end subroutine

 61:        subroutine ComputeMatrix(ksp,A,B,str,dummy,ierr)
 62:        implicit none
 63: #include <finclude/petscsys.h>
 64: #include <finclude/petscvec.h>
 65: #include <finclude/petscmat.h>
 66:        PetscErrorCode  ierr
 67:        KSP ksp
 68:        Mat A,B
 69:        MatStructure str
 70:        integer dummy(*)
 71:        DM dm

 73:       PetscInt    i,j,mx,my,xm
 74:       PetscInt    ym,xs,ys,i1,i5
 75:       PetscScalar  v(5),Hx,Hy
 76:       PetscScalar  HxdHy,HydHx
 77:       MatStencil   row(4),col(4,5)

 79:       i1 = 1
 80:       i5 = 5
 81:       call KSPGetDM(ksp,dm,ierr)
 82:       call DMDAGetInfo(dm,PETSC_NULL_INTEGER,mx,my,PETSC_NULL_INTEGER,  &
 83:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,             &
 84:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,             &
 85:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,             &
 86:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,             &
 87:      &               PETSC_NULL_INTEGER,ierr)

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