Actual source code: ex57f.F90
1: !
2: ! Description: Modified from ex2f.F and ex52.c to illustrate how use external packages MUMPS
3: ! Solves a linear system in parallel with KSP (Fortran code).
4: ! Also shows how to set a user-defined monitoring routine.
5: !
6: ! -----------------------------------------------------------------------
7: #include <petsc/finclude/petscksp.h>
8: module ex57fmodule
9: use petscksp
10: implicit none
12: contains
13: ! --------------------------------------------------------------
14: !
15: ! MyKSPMonitor - This is a user-defined routine for monitoring
16: ! the KSP iterative solvers.
17: !
18: ! Input Parameters:
19: ! ksp - iterative context
20: ! n - iteration number
21: ! rnorm - 2-norm (preconditioned) residual value (may be estimated)
22: ! unused - optional user-defined monitor context (unused here)
23: !
24: subroutine MyKSPMonitor(ksp, n, rnorm, unused, ierr)
26: KSP ksp
27: Vec x
28: PetscErrorCode ierr
29: PetscInt n, unused
30: PetscMPIInt rank
31: PetscReal rnorm
33: ! Build the solution vector
35: PetscCallA(KSPBuildSolution(ksp, PETSC_NULL_VEC, x, ierr))
37: ! Write the solution vector and residual norm to stdout
38: ! - Note that the parallel viewer PETSC_VIEWER_STDOUT_WORLD
39: ! handles data from multiple processors so that the
40: ! output is not jumbled.
42: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
43: if (rank == 0) write (6, 100) n
44: PetscCallA(VecView(x, PETSC_VIEWER_STDOUT_WORLD, ierr))
45: if (rank == 0) write (6, 200) n, rnorm
47: 100 format('iteration ', i5, ' solution vector:')
48: 200 format('iteration ', i5, ' residual norm ', e11.4)
49: ierr = 0
50: end
52: ! --------------------------------------------------------------
53: !
54: ! MyKSPConverged - This is a user-defined routine for testing
55: ! convergence of the KSP iterative solvers.
56: !
57: ! Input Parameters:
58: ! ksp - iterative context
59: ! n - iteration number
60: ! rnorm - 2-norm (preconditioned) residual value (may be estimated)
61: ! unused - optional user-defined monitor context (unused here)
62: !
63: subroutine MyKSPConverged(ksp, n, rnorm, flag, unused, ierr)
65: KSP ksp
66: PetscErrorCode, intent(out) :: ierr
67: PetscInt n, unused
68: KSPConvergedReason, intent(out) :: flag
69: PetscReal, intent(in) :: rnorm
71: if (rnorm <= .05) then
72: flag = KSP_CONVERGED_RTOL
73: else
74: flag = KSP_CONVERGED_ITERATING
75: end if
76: ierr = 0
78: end
79: end module ex57fmodule
81: program main
82: use petscksp
83: use ex57fmodule
84: implicit none
86: !
87: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
88: ! Variable declarations
89: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
90: !
91: ! Variables:
92: ! ksp - linear solver context
93: ! ksp - Krylov subspace method context
94: ! pc - preconditioner context
95: ! x, b, u - approx solution, right-hand side, exact solution vectors
96: ! A - matrix that defines linear system
97: ! its - iterations for convergence
98: ! norm - norm of error in solution
99: ! rctx - random number generator context
100: !
101: ! Note that vectors are declared as PETSc "Vec" objects. These vectors
102: ! are mathematical objects that contain more than just an array of
103: ! double precision numbers. I.e., vectors in PETSc are not just
104: ! double precision x(*).
105: ! However, local vector data can be easily accessed via VecGetArray().
106: ! See the Fortran section of the PETSc users manual for details.
107: !
108: #ifdef PETSC_HAVE_MUMPS
109: PetscInt icntl, ival
110: Mat F
111: #endif
112: PC pc
113: PetscReal norm
114: PetscInt i, j, II, JJ, m, n, its
115: PetscInt Istart, Iend
116: PetscErrorCode ierr
117: PetscMPIInt rank, size
118: PetscBool flg
119: PetscScalar v
120: PetscScalar, parameter :: one = 1.0, neg_one = -1.0
121: Vec x, b, u
122: Mat A
123: KSP ksp
124: PetscRandom rctx
125: character(len=80) ksptype
127: ! These variables are not currently used.
128: ! PC pc
129: ! PCType ptype
130: ! double precision tol
132: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
133: ! Beginning of program
134: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136: PetscCallA(PetscInitialize(ierr))
137: m = 3
138: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-m', m, flg, ierr))
139: n = 3
140: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', n, flg, ierr))
141: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
142: PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, size, ierr))
144: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
145: ! Compute the matrix and right-hand-side vector that define
146: ! the linear system, Ax = b.
147: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
149: ! Create parallel matrix, specifying only its global dimensions.
150: ! When using MatCreate(), the matrix format can be specified at
151: ! runtime. Also, the parallel partitioning of the matrix is
152: ! determined by PETSc at runtime.
154: PetscCallA(MatCreate(PETSC_COMM_WORLD, A, ierr))
155: PetscCallA(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m*n, m*n, ierr))
156: PetscCallA(MatSetFromOptions(A, ierr))
157: PetscCallA(MatSetUp(A, ierr))
159: ! Currently, all PETSc parallel matrix formats are partitioned by
160: ! contiguous chunks of rows across the processors. Determine which
161: ! rows of the matrix are locally owned.
163: PetscCallA(MatGetOwnershipRange(A, Istart, Iend, ierr))
165: ! Set matrix elements for the 2-D, five-point stencil in parallel.
166: ! - Each processor needs to insert only elements that it owns
167: ! locally (but any non-local elements will be sent to the
168: ! appropriate processor during matrix assembly).
169: ! - Always specify global row and columns of matrix entries.
170: ! - Note that MatSetValues() uses 0-based row and column numbers
171: ! in Fortran as well as in C.
173: ! Note: this uses the less common natural ordering that orders first
174: ! all the unknowns for x = h then for x = 2h etc; Hence you see JH = II +- n
175: ! instead of JJ = II +- m as you might expect. The more standard ordering
176: ! would first do all variables for y = h, then y = 2h etc.
178: do II = Istart, Iend - 1
179: v = -1.0
180: i = II/n
181: j = II - i*n
182: if (i > 0) then
183: JJ = II - n
184: PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, [II], 1_PETSC_INT_KIND, [JJ], [v], INSERT_VALUES, ierr))
185: end if
186: if (i < m - 1) then
187: JJ = II + n
188: PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, [II], 1_PETSC_INT_KIND, [JJ], [v], INSERT_VALUES, ierr))
189: end if
190: if (j > 0) then
191: JJ = II - 1
192: PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, [II], 1_PETSC_INT_KIND, [JJ], [v], INSERT_VALUES, ierr))
193: end if
194: if (j < n - 1) then
195: JJ = II + 1
196: PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, [II], 1_PETSC_INT_KIND, [JJ], [v], INSERT_VALUES, ierr))
197: end if
198: v = 4.0
199: PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, [II], 1_PETSC_INT_KIND, [II], [v], INSERT_VALUES, ierr))
200: end do
201: PetscCallA(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
202: PetscCallA(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))
204: ! Check if A is symmetric
205: if (size == 1) then
206: PetscCallA(MatIsSymmetric(A, 0.0_PETSC_REAL_KIND, flg, ierr))
207: if (flg .eqv. PETSC_FALSE) then
208: write (6, 120)
209: end if
210: end if
212: ! Create parallel vectors.
213: ! - Here, the parallel partitioning of the vector is determined by
214: ! PETSc at runtime. We could also specify the local dimensions
215: ! if desired -- or use the more general routine VecCreate().
216: ! - When solving a linear system, the vectors and matrices MUST
217: ! be partitioned accordingly. PETSc automatically generates
218: ! appropriately partitioned matrices and vectors when MatCreate()
219: ! and VecCreate() are used with the same communicator.
220: ! - Note: We form 1 vector from scratch and then duplicate as needed.
222: PetscCallA(VecCreateFromOptions(PETSC_COMM_WORLD, PETSC_NULL_CHARACTER, 1_PETSC_INT_KIND, PETSC_DECIDE, m*n, u, ierr))
223: PetscCallA(VecSetFromOptions(u, ierr))
224: PetscCallA(VecDuplicate(u, b, ierr))
225: PetscCallA(VecDuplicate(b, x, ierr))
227: ! Set exact solution; then compute right-hand-side vector.
228: ! By default we use an exact solution of a vector with all
229: ! elements of 1.0; Alternatively, using the runtime option
230: ! -random_sol forms a solution vector with random components.
232: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-random_exact_sol', flg, ierr))
233: if (flg) then
234: PetscCallA(PetscRandomCreate(PETSC_COMM_WORLD, rctx, ierr))
235: PetscCallA(PetscRandomSetFromOptions(rctx, ierr))
236: PetscCallA(VecSetRandom(u, rctx, ierr))
237: PetscCallA(PetscRandomDestroy(rctx, ierr))
238: else
239: PetscCallA(VecSet(u, one, ierr))
240: end if
241: PetscCallA(MatMult(A, u, b, ierr))
243: ! View the exact solution vector if desired
245: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-view_exact_sol', flg, ierr))
246: if (flg) then
247: PetscCallA(VecView(u, PETSC_VIEWER_STDOUT_WORLD, ierr))
248: end if
250: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
251: ! Create the linear solver and set various options
252: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
254: ! Create linear solver context
256: PetscCallA(KSPCreate(PETSC_COMM_WORLD, ksp, ierr))
258: ! Set operators. Here the matrix that defines the linear system
259: ! also serves as the matrix from which the preconditioner is constructed.
261: PetscCallA(KSPSetOperators(ksp, A, A, ierr))
263: PetscCallA(KSPSetType(ksp, KSPPREONLY, ierr))
264: PetscCallA(KSPGetType(ksp, ksptype, ierr))
265: PetscCheckA(ksptype == KSPPREONLY, PETSC_COMM_WORLD, PETSC_ERR_PLIB, 'Error')
266: PetscCallA(KSPGetPC(ksp, pc, ierr))
267: PetscCallA(PCSetType(pc, PCCHOLESKY, ierr))
268: #ifdef PETSC_HAVE_MUMPS
269: PetscCallA(PCFactorSetMatSolverType(pc, MATSOLVERMUMPS, ierr))
270: PetscCallA(PCFactorSetUpMatSolverType(pc, ierr))
271: PetscCallA(PCFactorGetMatrix(pc, F, ierr))
272: PetscCallA(KSPSetFromOptions(ksp, ierr))
273: icntl = 7; ival = 2
274: PetscCallA(MatMumpsSetIcntl(F, icntl, ival, ierr))
275: #endif
277: ! Set runtime options, e.g.,
278: ! -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
279: ! These options will override those specified above as long as
280: ! KSPSetFromOptions() is called _after_ any other customization
281: ! routines.
283: PetscCallA(KSPSetFromOptions(ksp, ierr))
285: ! Set convergence test routine if desired
287: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-my_ksp_convergence', flg, ierr))
288: if (flg) then
289: PetscCallA(KSPSetConvergenceTest(ksp, MyKSPConverged, 0, PETSC_NULL_FUNCTION, ierr))
290: end if
291: !
292: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
293: ! Solve the linear system
294: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
296: PetscCallA(KSPSolve(ksp, b, x, ierr))
298: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
299: ! Check solution and clean up
300: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
302: ! Check the error
303: PetscCallA(VecAXPY(x, neg_one, u, ierr))
304: PetscCallA(VecNorm(x, NORM_2, norm, ierr))
305: PetscCallA(KSPGetIterationNumber(ksp, its, ierr))
306: if (rank == 0) then
307: write (6, 100) norm, its
308: end if
309: 100 format('Norm of error ', e11.4, ' iterations ', i5)
310: 120 format('Matrix A is non-symmetric ')
312: ! Free work space. All PETSc objects should be destroyed when they
313: ! are no longer needed.
315: PetscCallA(KSPDestroy(ksp, ierr))
316: PetscCallA(VecDestroy(u, ierr))
317: PetscCallA(VecDestroy(x, ierr))
318: PetscCallA(VecDestroy(b, ierr))
319: PetscCallA(MatDestroy(A, ierr))
321: ! Always call PetscFinalize() before exiting a program. This routine
322: ! - finalizes the PETSc libraries as well as MPI
323: ! - provides summary and diagnostic information if certain runtime
324: ! options are chosen (e.g., -log_view). See PetscFinalize()
325: ! manpage for more information.
327: PetscCallA(PetscFinalize(ierr))
328: end
330: !/*TEST
331: !
332: ! test:
333: !
334: !TEST*/