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 .eq. 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 preconditioning matrix.

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_DEFAULT_REAL,PETSC_DEFAULT_REAL,PETSC_DEFAULT_INTEGER,ierr))

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

164:       PetscCallA(KSPSetFromOptions(ksp,ierr))

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

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

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

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

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

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

184: !  Check the error

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

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

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

207:       end

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