1:       program main   !   Solves the linear system  J x = f
  2: #include <finclude/petscdef.h>
  3:       use petscksp; use petscdm
  4:       Vec x,f
  5:       Mat J
  6:       DM da
  7:       KSP ksp
  8:       PetscErrorCode ierr
  9:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)

 11:       call DMDACreate1d(MPI_COMM_WORLD,DMDA_BOUNDARY_NONE,8,1,1,        &
 12:      &  PETSC_NULL_INTEGER,da,ierr)
 13:       call DMCreateGlobalVector(da,x,ierr)
 14:       call VecDuplicate(x,f,ierr)
 15:       call DMCreateMatrix(da,MATAIJ,J,ierr)

 17:       call ComputeRHS(da,f,ierr)
 18:       call ComputeMatrix(da,J,ierr)

 20:       call KSPCreate(MPI_COMM_WORLD,ksp,ierr)
 21:       call KSPSetOperators(ksp,J,J,SAME_NONZERO_PATTERN,ierr)
 22:       call KSPSetFromOptions(ksp,ierr)
 23:       call KSPSolve(ksp,f,x,ierr)

 25:       call MatDestroy(J,ierr)
 26:       call VecDestroy(x,ierr)
 27:       call VecDestroy(f,ierr)
 28:       call KSPDestroy(ksp,ierr)
 29:       call DMDestroy(da,ierr)
 30:       call PetscFinalize(ierr)
 31:       end
 32:       subroutine  ComputeRHS(da,x,ierr)
 33: #include <finclude/petscdef.h>
 34:       use petscdm
 35:       DM da
 36:       Vec x
 37:       PetscErrorCode ierr
 38:       PetscInt xs,xm,i,mx
 39:       PetscScalar hx
 40:       PetscScalar, pointer :: xx(:)
 41:       call DMDAGetInfo(da,PETSC_NULL_INTEGER,mx,PETSC_NULL_INTEGER,     &
 42:      &  PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,       &
 43:      &  PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,       &
 44:      &  PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,       &
 45:      &  PETSC_NULL_INTEGER,ierr)
 46:       call DMDAGetCorners(da,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,  &
 47:      &  xm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)
 48:       hx     = 1.d0/(mx-1)
 49:       call VecGetArrayF90(x,xx,ierr)
 50:       do i=xs,xs+xm-1
 51:           xx(i) = i*hx
 52:       enddo
 53:       call VecRestoreArrayF90(x,xx,ierr)
 54:       return
 55:       end
 56:       subroutine ComputeMatrix(da,J,ierr)
 57: #include <finclude/petscdef.h>
 58:       use petscdm
 59:       Mat J
 60:       DM da
 61:       PetscErrorCode ierr
 62:       PetscInt xs,xm,i,mx
 63:       PetscScalar hx
 64:       call DMDAGetInfo(da,PETSC_NULL_INTEGER,mx,PETSC_NULL_INTEGER,     &
 65:      &  PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,       &
 66:      &  PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,       &
 67:      &  PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,       &
 68:      &  PETSC_NULL_INTEGER,ierr)
 69:       call DMDAGetCorners(da,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,  &
 70:      &  xm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)
 71:       hx     = 1.d0/(mx-1)
 72:       do i=xs,xs+xm-1
 73:         if ((i .eq. 0) .or. (i .eq. mx-1)) then
 74:           call MatSetValue(J,i,i,1d0,INSERT_VALUES,ierr)
 75:         else
 76:           call MatSetValue(J,i,i-1,-hx,INSERT_VALUES,ierr)
 77:           call MatSetValue(J,i,i+1,-hx,INSERT_VALUES,ierr)
 78:           call MatSetValue(J,i,i,2*hx,INSERT_VALUES,ierr)
 79:         endif
 80:       enddo
 81:       call MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY,ierr)
 82:       call MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY,ierr)
 83:       return
 84:       end