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: 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: none = -1.0
68: one = 1.0
69: n = 10
70: i1 = 1
71: i2 = 2
72: i3 = 3
73: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', n, flg, ierr))
75: PetscCallA(PetscLogStageRegister('MatVec Assembly', stages(1), ierr))
76: PetscCallA(PetscLogStageRegister('KSP Solve', stages(2), ierr))
77: PetscCallA(PetscLogStagePush(stages(1), ierr))
78: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
79: ! Compute the matrix and right-hand-side vector that define
80: ! the linear system, Ax = b.
81: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
83: ! Create matrix. When using MatCreate(), the matrix format can
84: ! be specified at runtime.
86: PetscCallA(MatCreate(PETSC_COMM_WORLD, A, ierr))
87: PetscCallA(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n, ierr))
88: PetscCallA(MatSetFromOptions(A, ierr))
89: PetscCallA(MatSetUp(A, ierr))
91: ! Assemble matrix.
92: ! - Note that MatSetValues() uses 0-based row and column numbers
93: ! in Fortran as well as in C (as set here in the array "col").
95: value(1) = -1.0
96: value(2) = 2.0
97: value(3) = -1.0
98: do i = 1, n - 2
99: col(1) = i - 1
100: col(2) = i
101: col(3) = i + 1
102: PetscCallA(MatSetValues(A, i1, [i], i3, col, value, INSERT_VALUES, ierr))
103: end do
104: i = n - 1
105: col(1) = n - 2
106: col(2) = n - 1
107: PetscCallA(MatSetValues(A, i1, [i], i2, col, value, INSERT_VALUES, ierr))
108: i = 0
109: col(1) = 0
110: col(2) = 1
111: value(1) = 2.0
112: value(2) = -1.0
113: PetscCallA(MatSetValues(A, i1, [i], i2, col, value, INSERT_VALUES, ierr))
114: PetscCallA(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
115: PetscCallA(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))
117: ! Create vectors. Note that we form 1 vector from scratch and
118: ! then duplicate as needed.
120: PetscCallA(VecCreate(PETSC_COMM_WORLD, x, ierr))
121: PetscCallA(VecSetSizes(x, PETSC_DECIDE, n, ierr))
122: PetscCallA(VecSetFromOptions(x, ierr))
123: PetscCallA(VecDuplicate(x, b, ierr))
124: PetscCallA(VecDuplicate(x, u, ierr))
126: ! Set exact solution; then compute right-hand-side vector.
128: PetscCallA(VecSet(u, one, ierr))
129: PetscCallA(MatMult(A, u, b, ierr))
130: PetscCallA(PetscLogStagePop(ierr))
131: PetscCallA(PetscLogStagePush(stages(2), ierr))
133: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
134: ! Create the linear solver and set various options
135: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
137: ! Create linear solver context
139: PetscCallA(KSPCreate(PETSC_COMM_WORLD, ksp, ierr))
141: ! Set operators. Here the matrix that defines the linear system
142: ! also serves as the matrix from which the preconditioner is constructed.
144: PetscCallA(KSPConvergedDefaultCreate(defaultctx, ierr))
145: ! use KSPConvergedDefaultDestroyCptr() since it has that is the Fortran interface definition
146: PetscCallA(KSPSetConvergenceTest(ksp, MyKSPConverged, defaultctx, KSPConvergedDefaultDestroyCptr, ierr))
147: PetscCallA(KSPSetOperators(ksp, A, A, ierr))
149: ! Set linear solver defaults for this problem (optional).
150: ! - By extracting the KSP and PC contexts from the KSP context,
151: ! we can then directly call any KSP and PC routines
152: ! to set various options.
153: ! - The following four statements are optional; all of these
154: ! parameters could alternatively be specified at runtime via
155: ! KSPSetFromOptions()
157: PetscCallA(KSPGetPC(ksp, pc, ierr))
158: PetscCallA(PCSetType(pc, PCJACOBI, ierr))
159: tol = .0000001
160: PetscCallA(KSPSetTolerances(ksp, tol, PETSC_CURRENT_REAL, PETSC_CURRENT_REAL, PETSC_CURRENT_INTEGER, ierr))
161: PetscCallA(KSPGetTolerances(ksp, PETSC_NULL_REAL, tol, PETSC_NULL_REAL, PETSC_NULL_INTEGER, ierr))
163: ! Set runtime options, e.g.,
164: ! -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
165: ! These options will override those specified above as long as
166: ! KSPSetFromOptions() is called _after_ any other customization
167: ! routines.
169: PetscCallA(KSPSetFromOptions(ksp, ierr))
171: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
172: ! Solve the linear system
173: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
175: PetscCallA(KSPSolve(ksp, b, x, ierr))
176: PetscCallA(PetscLogStagePop(ierr))
178: ! View solver converged reason; we could instead use the option -ksp_converged_reason
179: PetscCallA(KSPConvergedReasonView(ksp, PETSC_VIEWER_STDOUT_WORLD, ierr))
181: ! View solver info; we could instead use the option -ksp_view
183: PetscCallA(KSPView(ksp, PETSC_VIEWER_STDOUT_WORLD, ierr))
185: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
186: ! Check solution and clean up
187: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
189: ! Check the error
191: PetscCallA(VecAXPY(x, none, u, ierr))
192: PetscCallA(VecNorm(x, NORM_2, norm, ierr))
193: PetscCallA(KSPGetIterationNumber(ksp, its, ierr))
194: if (norm > 1.e-12) then
195: write (6, 100) norm, its
196: else
197: write (6, 200) its
198: end if
199: 100 format('Norm of error ', e11.4, ', Iterations = ', i5)
200: 200 format('Norm of error < 1.e-12, Iterations = ', i5)
202: ! Free work space. All PETSc objects should be destroyed when they
203: ! are no longer needed.
205: PetscCallA(VecDestroy(x, ierr))
206: PetscCallA(VecDestroy(u, ierr))
207: PetscCallA(VecDestroy(b, ierr))
208: PetscCallA(MatDestroy(A, ierr))
209: PetscCallA(KSPDestroy(ksp, ierr))
210: PetscCallA(PetscFinalize(ierr))
212: end
214: !/*TEST
215: !
216: ! test:
217: ! requires: !single
218: ! args: -ksp_monitor_short
219: !
220: !TEST*/