Actual source code: ex1f.F90

      1: !
      2: !   Description: Solves a tridiagonal linear system with KSP.
      3: !
      4: ! -----------------------------------------------------------------------
      5: !
      6: !  Demonstrate a custom KSP convergence test that calls the default convergence test
      7: !
      8: #include <petsc/finclude/petscksp.h>
      9: module ex1fmodule
     10:   use petscksp
     11:   implicit none

     13: contains
     14:   subroutine MyKSPConverged(ksp, n, rnorm, flag, defaultctx, ierr)
     15:     KSP ksp
     16:     PetscErrorCode, intent(out) :: ierr
     17:     PetscInt n
     18:     integer(8) defaultctx
     19:     KSPConvergedReason flag
     20:     PetscReal rnorm

     22:     ! Must call default convergence test on the 0th iteration
     23:     PetscCall(KSPConvergedDefault(ksp, n, rnorm, flag, defaultctx, ierr))
     24:   end subroutine MyKSPConverged
     25: end module ex1fmodule

     27: program main
     28:   use petscksp
     29:   use ex1fmodule
     30:   implicit none

     32: !
     33: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     34: !                   Variable declarations
     35: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     36: !
     37: !  Variables:
     38: !     ksp     - linear solver context
     39: !     ksp      - Krylov subspace method context
     40: !     pc       - preconditioner context
     41: !     x, b, u  - approx solution, right-hand side, exact solution vectors
     42: !     A        - matrix that defines linear system
     43: !     its      - iterations for convergence
     44: !     norm     - norm of error in solution
     45: !
     46:   Vec x, b, u
     47:   Mat A
     48:   KSP ksp
     49:   PC pc
     50:   PetscReal norm, tol
     51:   PetscErrorCode ierr
     52:   PetscInt i, n, col(3), its
     53:   PetscBool flg
     54:   PetscMPIInt size
     55:   PetscScalar, parameter :: none = -1.0, one = 1.0
     56:   PetscScalar value(3)
     57:   PetscLogStage stages(2)
     58:   PetscFortranAddr defaultctx

     60: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     61: !                 Beginning of program
     62: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     64:   PetscCallA(PetscInitialize(ierr))
     65:   PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, size, ierr))
     66:   PetscCheckA(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, 'This is a uniprocessor example only')
     67:   n = 10
     68:   PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', n, flg, ierr))

     70:   PetscCallA(PetscLogStageRegister('MatVec Assembly', stages(1), ierr))
     71:   PetscCallA(PetscLogStageRegister('KSP Solve', stages(2), ierr))
     72:   PetscCallA(PetscLogStagePush(stages(1), ierr))
     73: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     74: !         Compute the matrix and right-hand-side vector that define
     75: !         the linear system, Ax = b.
     76: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     78: !  Create matrix.  When using MatCreate(), the matrix format can
     79: !  be specified at runtime.

     81:   PetscCallA(MatCreate(PETSC_COMM_WORLD, A, ierr))
     82:   PetscCallA(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n, ierr))
     83:   PetscCallA(MatSetFromOptions(A, ierr))
     84:   PetscCallA(MatSetUp(A, ierr))

     86: !  Assemble matrix.
     87: !   - Note that MatSetValues() uses 0-based row and column numbers
     88: !     in Fortran as well as in C (as set here in the array "col").

     90:   value = [-1.0, 2.0, -1.0]
     91:   do i = 1, n - 2
     92:     col(1) = i - 1
     93:     col(2) = i
     94:     col(3) = i + 1
     95:     PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, [i], 3_PETSC_INT_KIND, col, value, INSERT_VALUES, ierr))
     96:   end do
     97:   i = n - 1
     98:   col(1) = n - 2
     99:   col(2) = n - 1
    100:   PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, [i], 2_PETSC_INT_KIND, col, value, INSERT_VALUES, ierr))
    101:   i = 0
    102:   col(1) = 0
    103:   col(2) = 1
    104:   value(1) = 2.0
    105:   value(2) = -1.0
    106:   PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, [i], 2_PETSC_INT_KIND, col, value, INSERT_VALUES, ierr))
    107:   PetscCallA(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
    108:   PetscCallA(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))

    110: !  Create vectors.  Note that we form 1 vector from scratch and
    111: !  then duplicate as needed.

    113:   PetscCallA(VecCreate(PETSC_COMM_WORLD, x, ierr))
    114:   PetscCallA(VecSetSizes(x, PETSC_DECIDE, n, ierr))
    115:   PetscCallA(VecSetFromOptions(x, ierr))
    116:   PetscCallA(VecDuplicate(x, b, ierr))
    117:   PetscCallA(VecDuplicate(x, u, ierr))

    119: !  Set exact solution; then compute right-hand-side vector.

    121:   PetscCallA(VecSet(u, one, ierr))
    122:   PetscCallA(MatMult(A, u, b, ierr))
    123:   PetscCallA(PetscLogStagePop(ierr))
    124:   PetscCallA(PetscLogStagePush(stages(2), ierr))

    126: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    127: !          Create the linear solver and set various options
    128: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    130: !  Create linear solver context

    132:   PetscCallA(KSPCreate(PETSC_COMM_WORLD, ksp, ierr))

    134: !  Set operators. Here the matrix that defines the linear system
    135: !  also serves as the matrix from which the preconditioner is constructed.

    137:   PetscCallA(KSPConvergedDefaultCreate(defaultctx, ierr))
    138:   ! use KSPConvergedDefaultDestroyCptr() since it has that is the Fortran interface definition
    139:   PetscCallA(KSPSetConvergenceTest(ksp, MyKSPConverged, defaultctx, KSPConvergedDefaultDestroyCptr, ierr))
    140:   PetscCallA(KSPSetOperators(ksp, A, A, ierr))

    142: !  Set linear solver defaults for this problem (optional).
    143: !   - By extracting the KSP and PC contexts from the KSP context,
    144: !     we can then directly call any KSP and PC routines
    145: !     to set various options.
    146: !   - The following four statements are optional; all of these
    147: !     parameters could alternatively be specified at runtime via
    148: !     KSPSetFromOptions()

    150:   PetscCallA(KSPGetPC(ksp, pc, ierr))
    151:   PetscCallA(PCSetType(pc, PCJACOBI, ierr))
    152:   tol = .0000001_PETSC_REAL_KIND
    153:   PetscCallA(KSPSetTolerances(ksp, tol, PETSC_CURRENT_REAL, PETSC_CURRENT_REAL, PETSC_CURRENT_INTEGER, ierr))
    154:   PetscCallA(KSPGetTolerances(ksp, PETSC_NULL_REAL, tol, PETSC_NULL_REAL, PETSC_NULL_INTEGER, ierr))

    156: !  Set runtime options, e.g.,
    157: !      -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
    158: !  These options will override those specified above as long as
    159: !  KSPSetFromOptions() is called _after_ any other customization
    160: !  routines.

    162:   PetscCallA(KSPSetFromOptions(ksp, ierr))

    164: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    165: !                      Solve the linear system
    166: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    168:   PetscCallA(KSPSolve(ksp, b, x, ierr))
    169:   PetscCallA(PetscLogStagePop(ierr))

    171: !  View solver converged reason; we could instead use the option -ksp_converged_reason
    172:   PetscCallA(KSPConvergedReasonView(ksp, PETSC_VIEWER_STDOUT_WORLD, ierr))

    174: !  View solver info; we could instead use the option -ksp_view

    176:   PetscCallA(KSPView(ksp, PETSC_VIEWER_STDOUT_WORLD, ierr))

    178: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    179: !                      Check solution and clean up
    180: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    182: !  Check the error

    184:   PetscCallA(VecAXPY(x, none, u, ierr))
    185:   PetscCallA(VecNorm(x, NORM_2, norm, ierr))
    186:   PetscCallA(KSPGetIterationNumber(ksp, its, ierr))
    187:   if (norm > 1.e-12) then
    188:     write (6, 100) norm, its
    189:   else
    190:     write (6, 200) its
    191:   end if
    192: 100 format('Norm of error ', e11.4, ',  Iterations = ', i5)
    193: 200 format('Norm of error < 1.e-12, Iterations = ', i5)

    195: !  Free work space.  All PETSc objects should be destroyed when they
    196: !  are no longer needed.

    198:   PetscCallA(VecDestroy(x, ierr))
    199:   PetscCallA(VecDestroy(u, ierr))
    200:   PetscCallA(VecDestroy(b, ierr))
    201:   PetscCallA(MatDestroy(A, ierr))
    202:   PetscCallA(KSPDestroy(ksp, ierr))
    203:   PetscCallA(PetscFinalize(ierr))

    205: end

    207: !/*TEST
    208: !
    209: !     test:
    210: !       requires: !single
    211: !       args: -ksp_monitor_short
    212: !
    213: !TEST*/