1: #include <petsc/finclude/petscksp.h>
  2: #include <petsc/finclude/petscdmda.h>
  3: program main              !   Solves the linear system  J x = f
  4:   use petscksp
  5:   use petscdmda
  6:   implicit none

  8:   Vec x, f
  9:   Mat J
 10:   DM da
 11:   KSP ksp
 12:   PetscErrorCode ierr
 13:   PetscInt eight, one

 15:   eight = 8
 16:   one = 1
 17:   PetscCallA(PetscInitialize(ierr))
 18:   PetscCallA(DMDACreate1d(MPI_COMM_WORLD, DM_BOUNDARY_NONE, eight, one, one, PETSC_NULL_INTEGER_ARRAY, da, ierr))
 19:   PetscCallA(DMSetFromOptions(da, ierr))
 20:   PetscCallA(DMSetUp(da, ierr))
 21:   PetscCallA(DMCreateGlobalVector(da, x, ierr))
 22:   PetscCallA(VecDuplicate(x, f, ierr))
 23:   PetscCallA(DMSetMatType(da, MATAIJ, ierr))
 24:   PetscCallA(DMCreateMatrix(da, J, ierr))

 26:   PetscCallA(ComputeRHS(da, f, ierr))
 27:   PetscCallA(ComputeMatrix(da, J, ierr))

 29:   PetscCallA(KSPCreate(MPI_COMM_WORLD, ksp, ierr))
 30:   PetscCallA(KSPSetOperators(ksp, J, J, ierr))
 31:   PetscCallA(KSPSetFromOptions(ksp, ierr))
 32:   PetscCallA(KSPSolve(ksp, f, x, ierr))

 34:   PetscCallA(MatDestroy(J, ierr))
 35:   PetscCallA(VecDestroy(x, ierr))
 36:   PetscCallA(VecDestroy(f, ierr))
 37:   PetscCallA(KSPDestroy(ksp, ierr))
 38:   PetscCallA(DMDestroy(da, ierr))
 39:   PetscCallA(PetscFinalize(ierr))

 41: contains
 42:   subroutine ComputeRHS(da, x, ierr)
 43:     use petscdmda
 44:     implicit none

 46:     DM da
 47:     Vec x
 48:     PetscErrorCode ierr
 49:     PetscInt xs, xm, i, mx
 50:     PetscScalar hx
 51:     PetscScalar, pointer :: xx(:)
 52:     PetscCall(DMDAGetInfo(da, PETSC_NULL_INTEGER, mx, PETSC_NULL_INTEGER, 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))
 53:     PetscCall(DMDAGetCorners(da, xs, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, xm, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
 54:     hx = 1.0_PETSC_REAL_KIND/(mx - 1)
 55:     PetscCall(DMDAVecGetArray(da, x, xx, ierr))
 56:     do i = xs, xs + xm - 1
 57:       xx(i) = i*hx
 58:     end do
 59:     PetscCall(DMDAVecRestoreArray(da, x, xx, ierr))
 60:   end

 62:   subroutine ComputeMatrix(da, J, ierr)
 63:     use petscdmda
 64:     use petscmat
 65:     implicit none

 67:     Mat J
 68:     DM da
 69:     PetscErrorCode ierr
 70:     PetscInt xs, xm, i, mx
 71:     PetscScalar hx, one

 73:     one = 1.0
 74:     PetscCall(DMDAGetInfo(da, PETSC_NULL_INTEGER, mx, PETSC_NULL_INTEGER, 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))
 75:     PetscCall(DMDAGetCorners(da, xs, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, xm, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
 76:     hx = 1.0_PETSC_REAL_KIND/(mx - 1)
 77:     do i = xs, xs + xm - 1
 78:       if ((i == 0) .or. (i == mx - 1)) then
 79:         PetscCall(MatSetValue(J, i, i, one, INSERT_VALUES, ierr))
 80:       else
 81:         PetscCall(MatSetValue(J, i, i - 1, -hx, INSERT_VALUES, ierr))
 82:         PetscCall(MatSetValue(J, i, i + 1, -hx, INSERT_VALUES, ierr))
 83:         PetscCall(MatSetValue(J, i, i, 2*hx, INSERT_VALUES, ierr))
 84:       end if
 85:     end do
 86:     PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY, ierr))
 87:     PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY, ierr))
 88:   end
 89: end program
 90: !/*TEST
 91: !
 92: !   test:
 93: !      args: -ksp_converged_reason
 94: !
 95: !TEST*/