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

 14:   PetscCallA(PetscInitialize(ierr))
 15:   PetscCallA(DMDACreate1d(MPI_COMM_WORLD, DM_BOUNDARY_NONE, 8_PETSC_INT_KIND, 1_PETSC_INT_KIND, 1_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, da, ierr))
 16:   PetscCallA(DMSetFromOptions(da, ierr))
 17:   PetscCallA(DMSetUp(da, ierr))
 18:   PetscCallA(DMCreateGlobalVector(da, x, ierr))
 19:   PetscCallA(VecDuplicate(x, f, ierr))
 20:   PetscCallA(DMSetMatType(da, MATAIJ, ierr))
 21:   PetscCallA(DMCreateMatrix(da, J, ierr))

 23:   PetscCallA(ComputeRHS(da, f, ierr))
 24:   PetscCallA(ComputeMatrix(da, J, ierr))

 26:   PetscCallA(KSPCreate(MPI_COMM_WORLD, ksp, ierr))
 27:   PetscCallA(KSPSetOperators(ksp, J, J, ierr))
 28:   PetscCallA(KSPSetFromOptions(ksp, ierr))
 29:   PetscCallA(KSPSolve(ksp, f, x, ierr))

 31:   PetscCallA(MatDestroy(J, ierr))
 32:   PetscCallA(VecDestroy(x, ierr))
 33:   PetscCallA(VecDestroy(f, ierr))
 34:   PetscCallA(KSPDestroy(ksp, ierr))
 35:   PetscCallA(DMDestroy(da, ierr))
 36:   PetscCallA(PetscFinalize(ierr))

 38: contains
 39:   subroutine ComputeRHS(da, x, ierr)
 40:     DM da
 41:     Vec x
 42:     PetscErrorCode, intent(out) :: ierr
 43:     PetscInt xs, xm, i, mx
 44:     PetscScalar hx
 45:     PetscScalar, pointer :: xx(:)
 46:     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))
 47:     PetscCall(DMDAGetCorners(da, xs, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, xm, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
 48:     hx = 1.0_PETSC_REAL_KIND/(mx - 1)
 49:     PetscCall(DMDAVecGetArray(da, x, xx, ierr))
 50:     do i = xs, xs + xm - 1
 51:       xx(i) = i*hx
 52:     end do
 53:     PetscCall(DMDAVecRestoreArray(da, x, xx, ierr))
 54:   end

 56:   subroutine ComputeMatrix(da, J, ierr)
 57:     Mat J
 58:     DM da
 59:     PetscErrorCode, intent(out) :: ierr
 60:     PetscInt xs, xm, i, mx
 61:     PetscScalar hx
 62:     PetscScalar, parameter :: one = 1.0

 64:     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))
 65:     PetscCall(DMDAGetCorners(da, xs, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, xm, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
 66:     hx = 1.0_PETSC_REAL_KIND/(mx - 1)
 67:     do i = xs, xs + xm - 1
 68:       if ((i == 0) .or. (i == mx - 1)) then
 69:         PetscCall(MatSetValue(J, i, i, one, INSERT_VALUES, ierr))
 70:       else
 71:         PetscCall(MatSetValue(J, i, i - 1, -hx, INSERT_VALUES, ierr))
 72:         PetscCall(MatSetValue(J, i, i + 1, -hx, INSERT_VALUES, ierr))
 73:         PetscCall(MatSetValue(J, i, i, 2*hx, INSERT_VALUES, ierr))
 74:       end if
 75:     end do
 76:     PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY, ierr))
 77:     PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY, ierr))
 78:   end
 79: end program
 80: !/*TEST
 81: !
 82: !   test:
 83: !      args: -ksp_converged_reason
 84: !
 85: !TEST*/