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: subroutine MyKSPConverged(ksp, n, rnorm, flag, defaultctx, ierr)
      9: #include <petsc/finclude/petscksp.h>
     10:   use petscksp

     12:   KSP ksp
     13:   PetscErrorCode ierr
     14:   PetscInt n
     15:   integer*8 defaultctx
     16:   KSPConvergedReason flag
     17:   PetscReal rnorm

     19:   ! Must call default convergence test on the 0th iteration
     20:   PetscCall(KSPConvergedDefault(ksp, n, rnorm, flag, defaultctx, ierr))
     21: end subroutine MyKSPConverged

     23: program main
     24: #include <petsc/finclude/petscksp.h>
     25:   use petscksp
     26:   implicit none

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

     56: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     57: !                 Beginning of program
     58: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     60:   PetscCallA(PetscInitialize(ierr))
     61:   PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, size, ierr))
     62:   PetscCheckA(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, 'This is a uniprocessor example only')
     63:   none = -1.0
     64:   one = 1.0
     65:   n = 10
     66:   i1 = 1
     67:   i2 = 2
     68:   i3 = 3
     69:   PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', n, flg, ierr))

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

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

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

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

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

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

    116:     PetscCallA(VecCreate(PETSC_COMM_WORLD, x, ierr))
    117:     PetscCallA(VecSetSizes(x, PETSC_DECIDE, n, ierr))
    118:     PetscCallA(VecSetFromOptions(x, ierr))
    119:     PetscCallA(VecDuplicate(x, b, ierr))
    120:     PetscCallA(VecDuplicate(x, u, ierr))

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

    124:     PetscCallA(VecSet(u, one, ierr))
    125:     PetscCallA(MatMult(A, u, b, ierr))
    126:     PetscCallA(PetscLogStagePop(ierr))
    127:     PetscCallA(PetscLogStagePush(stages(2), ierr))

    129: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    130: !          Create the linear solver and set various options
    131: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    133: !  Create linear solver context

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

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

    140:     PetscCallA(KSPConvergedDefaultCreate(defaultctx, ierr))
    141:     PetscCallA(KSPSetConvergenceTest(ksp, MyKSPConverged, defaultctx, KSPConvergedDefaultDestroy, ierr))

    143:     PetscCallA(KSPSetOperators(ksp, A, A, ierr))

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

    153:     PetscCallA(KSPGetPC(ksp, pc, ierr))
    154:     PetscCallA(PCSetType(pc, PCJACOBI, ierr))
    155:     tol = .0000001
    156:     PetscCallA(KSPSetTolerances(ksp, tol, PETSC_CURRENT_REAL, PETSC_CURRENT_REAL, PETSC_CURRENT_INTEGER, ierr))
    157:     PetscCallA(KSPGetTolerances(ksp, PETSC_NULL_REAL, tol, PETSC_NULL_REAL, PETSC_NULL_INTEGER, ierr))

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

    165:     PetscCallA(KSPSetFromOptions(ksp, ierr))

    167: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    168: !                      Solve the linear system
    169: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    171:     PetscCallA(KSPSolve(ksp, b, x, ierr))
    172:     PetscCallA(PetscLogStagePop(ierr))

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

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

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

    181: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    182: !                      Check solution and clean up
    183: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    185: !  Check the error

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

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

    201:     PetscCallA(VecDestroy(x, ierr))
    202:     PetscCallA(VecDestroy(u, ierr))
    203:     PetscCallA(VecDestroy(b, ierr))
    204:     PetscCallA(MatDestroy(A, ierr))
    205:     PetscCallA(KSPDestroy(ksp, ierr))
    206:     PetscCallA(PetscFinalize(ierr))

    208:   end

    210: !/*TEST
    211: !
    212: !     test:
    213: !       requires: !single
    214: !       args: -ksp_monitor_short
    215: !
    216: !TEST*/