Actual source code: ex44f.F90

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

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

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

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

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

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

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

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

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

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

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