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*/