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)
16: KSP ksp
17: PetscErrorCode ierr
18: PetscInt n
19: integer(8) defaultctx
20: KSPConvergedReason flag
21: PetscReal rnorm
23: ! Must call default convergence test on the 0th iteration
24: PetscCall(KSPConvergedDefault(ksp, n, rnorm, flag, defaultctx, ierr))
25: end subroutine MyKSPConverged
26: end module ex1fmodule
28: program main
29: use petscksp
30: use ex1fmodule
31: implicit none
33: !
34: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
35: ! Variable declarations
36: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
37: !
38: ! Variables:
39: ! ksp - linear solver context
40: ! ksp - Krylov subspace method context
41: ! pc - preconditioner context
42: ! x, b, u - approx solution, right-hand side, exact solution vectors
43: ! A - matrix that defines linear system
44: ! its - iterations for convergence
45: ! norm - norm of error in solution
46: !
47: Vec x, b, u
48: Mat A
49: KSP ksp
50: PC pc
51: PetscReal norm, tol
52: PetscErrorCode ierr
53: PetscInt i, n, col(3), its, i1, i2, i3
54: PetscBool flg
55: PetscMPIInt size
56: PetscScalar none, one, value(3)
57: PetscLogStage stages(2)
58: integer(8) defaultctx
59: external :: KSPConvergedDefaultDestroy ! has no interface (yet)
61: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
62: ! Beginning of program
63: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
65: PetscCallA(PetscInitialize(ierr))
66: PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, size, ierr))
67: PetscCheckA(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, 'This is a uniprocessor example only')
68: none = -1.0
69: one = 1.0
70: n = 10
71: i1 = 1
72: i2 = 2
73: i3 = 3
74: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', n, flg, ierr))
76: PetscCallA(PetscLogStageRegister('MatVec Assembly', stages(1), ierr))
77: PetscCallA(PetscLogStageRegister('KSP Solve', stages(2), ierr))
78: PetscCallA(PetscLogStagePush(stages(1), ierr))
79: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
80: ! Compute the matrix and right-hand-side vector that define
81: ! the linear system, Ax = b.
82: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
84: ! Create matrix. When using MatCreate(), the matrix format can
85: ! be specified at runtime.
87: PetscCallA(MatCreate(PETSC_COMM_WORLD, A, ierr))
88: PetscCallA(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n, ierr))
89: PetscCallA(MatSetFromOptions(A, ierr))
90: PetscCallA(MatSetUp(A, ierr))
92: ! Assemble matrix.
93: ! - Note that MatSetValues() uses 0-based row and column numbers
94: ! in Fortran as well as in C (as set here in the array "col").
96: value(1) = -1.0
97: value(2) = 2.0
98: value(3) = -1.0
99: do i = 1, n - 2
100: col(1) = i - 1
101: col(2) = i
102: col(3) = i + 1
103: PetscCallA(MatSetValues(A, i1, [i], i3, col, value, INSERT_VALUES, ierr))
104: end do
105: i = n - 1
106: col(1) = n - 2
107: col(2) = n - 1
108: PetscCallA(MatSetValues(A, i1, [i], i2, col, value, INSERT_VALUES, ierr))
109: i = 0
110: col(1) = 0
111: col(2) = 1
112: value(1) = 2.0
113: value(2) = -1.0
114: PetscCallA(MatSetValues(A, i1, [i], i2, col, value, INSERT_VALUES, ierr))
115: PetscCallA(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
116: PetscCallA(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))
118: ! Create vectors. Note that we form 1 vector from scratch and
119: ! then duplicate as needed.
121: PetscCallA(VecCreate(PETSC_COMM_WORLD, x, ierr))
122: PetscCallA(VecSetSizes(x, PETSC_DECIDE, n, ierr))
123: PetscCallA(VecSetFromOptions(x, ierr))
124: PetscCallA(VecDuplicate(x, b, ierr))
125: PetscCallA(VecDuplicate(x, u, ierr))
127: ! Set exact solution; then compute right-hand-side vector.
129: PetscCallA(VecSet(u, one, ierr))
130: PetscCallA(MatMult(A, u, b, ierr))
131: PetscCallA(PetscLogStagePop(ierr))
132: PetscCallA(PetscLogStagePush(stages(2), ierr))
134: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
135: ! Create the linear solver and set various options
136: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138: ! Create linear solver context
140: PetscCallA(KSPCreate(PETSC_COMM_WORLD, ksp, ierr))
142: ! Set operators. Here the matrix that defines the linear system
143: ! also serves as the matrix from which the preconditioner is constructed.
145: PetscCallA(KSPConvergedDefaultCreate(defaultctx, ierr))
146: PetscCallA(KSPSetConvergenceTest(ksp, MyKSPConverged, defaultctx, KSPConvergedDefaultDestroy, ierr))
148: PetscCallA(KSPSetOperators(ksp, A, A, ierr))
150: ! Set linear solver defaults for this problem (optional).
151: ! - By extracting the KSP and PC contexts from the KSP context,
152: ! we can then directly call any KSP and PC routines
153: ! to set various options.
154: ! - The following four statements are optional; all of these
155: ! parameters could alternatively be specified at runtime via
156: ! KSPSetFromOptions()
158: PetscCallA(KSPGetPC(ksp, pc, ierr))
159: PetscCallA(PCSetType(pc, PCJACOBI, ierr))
160: tol = .0000001
161: PetscCallA(KSPSetTolerances(ksp, tol, PETSC_CURRENT_REAL, PETSC_CURRENT_REAL, PETSC_CURRENT_INTEGER, ierr))
162: PetscCallA(KSPGetTolerances(ksp, PETSC_NULL_REAL, tol, PETSC_NULL_REAL, PETSC_NULL_INTEGER, ierr))
164: ! Set runtime options, e.g.,
165: ! -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
166: ! These options will override those specified above as long as
167: ! KSPSetFromOptions() is called _after_ any other customization
168: ! routines.
170: PetscCallA(KSPSetFromOptions(ksp, ierr))
172: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
173: ! Solve the linear system
174: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
176: PetscCallA(KSPSolve(ksp, b, x, ierr))
177: PetscCallA(PetscLogStagePop(ierr))
179: ! View solver converged reason; we could instead use the option -ksp_converged_reason
180: PetscCallA(KSPConvergedReasonView(ksp, PETSC_VIEWER_STDOUT_WORLD, ierr))
182: ! View solver info; we could instead use the option -ksp_view
184: PetscCallA(KSPView(ksp, PETSC_VIEWER_STDOUT_WORLD, ierr))
186: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
187: ! Check solution and clean up
188: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
190: ! Check the error
192: PetscCallA(VecAXPY(x, none, u, ierr))
193: PetscCallA(VecNorm(x, NORM_2, norm, ierr))
194: PetscCallA(KSPGetIterationNumber(ksp, its, ierr))
195: if (norm > 1.e-12) then
196: write (6, 100) norm, its
197: else
198: write (6, 200) its
199: end if
200: 100 format('Norm of error ', e11.4, ', Iterations = ', i5)
201: 200 format('Norm of error < 1.e-12, Iterations = ', i5)
203: ! Free work space. All PETSc objects should be destroyed when they
204: ! are no longer needed.
206: PetscCallA(VecDestroy(x, ierr))
207: PetscCallA(VecDestroy(u, ierr))
208: PetscCallA(VecDestroy(b, ierr))
209: PetscCallA(MatDestroy(A, ierr))
210: PetscCallA(KSPDestroy(ksp, ierr))
211: PetscCallA(PetscFinalize(ierr))
213: end
215: !/*TEST
216: !
217: ! test:
218: ! requires: !single
219: ! args: -ksp_monitor_short
220: !
221: !TEST*/