1:       program main   !   Solves the linear system  J x = f
  2: #include <petsc/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,DM_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 DMSetMatType(da,MATAIJ,ierr)
 16:       call DMCreateMatrix(da,J,ierr)

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

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

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