Actual source code: ex62f.F90
1: !
2: ! Solves a linear system in parallel with KSP. Also indicates
3: ! use of a user-provided preconditioner. Input parameters include:
4: !
5: !
7: !
8: ! -------------------------------------------------------------------------
9: #include <petsc/finclude/petscksp.h>
10: module ex62fmodule
11: use petscksp
12: implicit none
13: PC jacobi, sor
14: Vec work
15: Mat matwork
16: PetscInt matapply_calls, matapplyrichardson_calls
18: contains
19: !/***********************************************************************/
20: !/* Routines for a user-defined shell preconditioner */
21: !/***********************************************************************/
23: !
24: ! SampleShellPCSetUp - This routine sets up a user-defined
25: ! preconditioner context.
26: !
27: ! Input Parameters:
28: ! pc - preconditioner object
29: ! x - vector
30: !
31: ! Output Parameter:
32: ! ierr - error code (nonzero if error has been detected)
33: !
34: ! Notes:
35: ! In this example, we define the shell preconditioner to be Jacobi
36: ! method. Thus, here we create a work vector for storing the reciprocal
37: ! of the diagonal of the matrix used to compute the preconditioner; this vector is then
38: ! used within the routine SampleShellPCApply().
39: !
40: subroutine SampleShellPCSetUp(pc, x, ierr)
42: PC pc
43: Vec x
44: Mat pmat
45: PetscErrorCode ierr
47: PetscCallA(PCGetOperators(pc, PETSC_NULL_MAT, pmat, ierr))
48: PetscCallA(PCCreate(PETSC_COMM_WORLD, jacobi, ierr))
49: PetscCallA(PCSetType(jacobi, PCJACOBI, ierr))
50: PetscCallA(PCSetOperators(jacobi, pmat, pmat, ierr))
51: PetscCallA(PCSetUp(jacobi, ierr))
53: PetscCallA(PCCreate(PETSC_COMM_WORLD, sor, ierr))
54: PetscCallA(PCSetType(sor, PCSOR, ierr))
55: PetscCallA(PCSetOperators(sor, pmat, pmat, ierr))
56: ! PetscCallA(PCSORSetSymmetric(sor,SOR_LOCAL_SYMMETRIC_SWEEP,ierr))
57: PetscCallA(PCSetUp(sor, ierr))
59: PetscCallA(VecDuplicate(x, work, ierr))
60: end
62: ! -------------------------------------------------------------------
63: !
64: ! SampleShellPCApply - This routine demonstrates the use of a
65: ! user-provided preconditioner.
66: !
67: ! Input Parameters:
68: ! pc - preconditioner object
69: ! x - input vector
70: !
71: ! Output Parameters:
72: ! y - preconditioned vector
73: ! ierr - error code (nonzero if error has been detected)
74: !
75: ! Notes:
76: ! This code implements the Jacobi preconditioner plus the
77: ! SOR preconditioner
78: !
79: ! YOU CAN GET THE EXACT SAME EFFECT WITH THE PCCOMPOSITE preconditioner using
80: ! mpiexec -n 1 ex21f -ksp_monitor -pc_type composite -pc_composite_pcs jacobi,sor -pc_composite_type additive
81: !
82: subroutine SampleShellPCApply(pc, x, y, ierr)
84: PC pc
85: Vec x, y
86: PetscErrorCode ierr
87: PetscScalar, parameter :: one = 1.0
89: PetscCallA(PCApply(jacobi, x, y, ierr))
90: PetscCallA(PCApply(sor, x, work, ierr))
91: PetscCallA(VecAXPY(y, one, work, ierr))
92: end
94: ! -------------------------------------------------------------------
95: !
96: ! SampleShellPCMatApply - Block analog of SampleShellPCApply(), the
97: ! preconditioner is applied to a whole block of vectors at once
98: !
99: ! Input Parameters:
100: ! pc - preconditioner object
101: ! x - input block of vectors
102: !
103: ! Output Parameters:
104: ! y - preconditioned block of vectors
105: ! ierr - error code (nonzero if error has been detected)
106: !
107: subroutine SampleShellPCMatApply(pc, x, y, ierr)
109: PC pc
110: Mat x, y
111: PetscErrorCode ierr
112: PetscScalar, parameter :: one = 1.0
114: PetscCallA(PCMatApply(jacobi, x, y, ierr))
115: PetscCallA(PCMatApply(sor, x, matwork, ierr))
116: PetscCallA(MatAXPY(y, one, matwork, SAME_NONZERO_PATTERN, ierr))
117: matapply_calls = matapply_calls + 1
118: end
120: ! -------------------------------------------------------------------
121: !
122: ! SampleShellPCMatApplyRichardson - Richardson iteration on a whole block
123: ! of right-hand sides, x <- x + M (b - A x) with M the preconditioner of
124: ! SampleShellPCMatApply(), so that it computes the same solutions as the
125: ! block iteration of KSPRICHARDSON
126: !
127: ! Input Parameters:
128: ! pc - preconditioner object
129: ! b - block of right-hand sides
130: ! x - block of initial guesses
131: ! w - block of work vectors, KSPRICHARDSON passes a null Mat so the work blocks are created here
132: ! rtol - relative tolerance, unused since the iteration is not tested for convergence
133: ! abstol - absolute tolerance, unused
134: ! dtol - divergence tolerance, unused
135: ! maxits - number of iterations to perform
136: ! guesszero - PETSC_TRUE if x is zero on entry
137: !
138: ! Output Parameters:
139: ! x - block of solutions
140: ! outits - number of iterations performed
141: ! reason - why the iteration stopped
142: ! ierr - error code (nonzero if error has been detected)
143: !
144: subroutine SampleShellPCMatApplyRichardson(pc, b, x, w, rtol, abstol, dtol, maxits, guesszero, outits, reason, ierr)
146: PC pc
147: Mat b, x, w
148: PetscReal rtol, abstol, dtol
149: PetscInt maxits, outits
150: PetscBool guesszero
151: PCRichardsonConvergedReason reason
152: PetscErrorCode ierr
153: Mat amat, r, z
154: PetscInt i
155: PetscScalar, parameter :: one = 1.0, neg_one = -1.0
157: PetscCallA(PCGetOperators(pc, amat, PETSC_NULL_MAT, ierr))
158: PetscCallA(MatDuplicate(b, MAT_DO_NOT_COPY_VALUES, z, ierr))
159: PetscCallA(MatDuplicate(b, MAT_DO_NOT_COPY_VALUES, r, ierr))
160: ! set the product A x up once so that only its numeric phase is run in the loop below
161: PetscCallA(MatProductCreateWithMat(amat, x, PETSC_NULL_MAT, r, ierr))
162: PetscCallA(MatProductSetType(r, MATPRODUCT_AB, ierr))
163: PetscCallA(MatProductSetFromOptions(r, ierr))
164: PetscCallA(MatProductSymbolic(r, ierr))
165: do i = 1, maxits
166: if (i == 1 .and. guesszero) then
167: PetscCallA(MatCopy(b, r, SAME_NONZERO_PATTERN, ierr))
168: else
169: PetscCallA(MatProductNumeric(r, ierr))
170: PetscCallA(MatAYPX(r, neg_one, b, SAME_NONZERO_PATTERN, ierr))
171: end if
172: PetscCallA(PCMatApply(jacobi, r, z, ierr))
173: PetscCallA(PCMatApply(sor, r, matwork, ierr))
174: PetscCallA(MatAXPY(z, one, matwork, SAME_NONZERO_PATTERN, ierr))
175: PetscCallA(MatAXPY(x, one, z, SAME_NONZERO_PATTERN, ierr))
176: end do
177: PetscCallA(MatProductClear(r, ierr))
178: PetscCallA(MatDestroy(r, ierr))
179: PetscCallA(MatDestroy(z, ierr))
180: outits = maxits
181: reason = PCRICHARDSON_CONVERGED_ITS
182: matapplyrichardson_calls = matapplyrichardson_calls + 1
183: end
185: end module
187: program main
188: use ex62fmodule
189: implicit none
191: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
192: ! Variable declarations
193: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
194: !
195: ! Variables:
196: ! ksp - linear solver context
197: ! ksp - Krylov subspace method context
198: ! pc - preconditioner context
199: ! x, b, u - approx solution, right-hand side, exact solution vectors
200: ! A - matrix that defines linear system
201: ! its - iterations for convergence
202: ! norm - norm of solution error
204: Vec x, b, u
205: Mat A, blockb, blockx, blocky
206: PC pc, blockpc
207: KSP ksp, blockksp
208: PetscScalar v
209: PetscScalar, parameter :: one = 1.0, neg_one = -1.0
210: PetscReal, parameter :: tol = 1e-7
211: PetscReal norm, blocknorm
212: PetscInt i, j, II, JJ, Istart, Iend, its
213: PetscInt m, n, mloc, nrhs, nits
214: PetscMPIInt rank
215: PetscBool flg, blocksolve
216: PetscErrorCode ierr
218: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
219: ! Beginning of program
220: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
222: PetscCallA(PetscInitialize(ierr))
223: m = 8
224: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-m', m, flg, ierr))
225: n = 7
226: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', n, flg, ierr))
227: blocksolve = PETSC_FALSE
228: PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-matsolve', blocksolve, flg, ierr))
229: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
231: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
232: ! Compute the matrix and right-hand-side vector that define
233: ! the linear system, Ax = b.
234: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
236: ! Create parallel matrix, specifying only its global dimensions.
237: ! When using MatCreate(), the matrix format can be specified at
238: ! runtime. Also, the parallel partitioning of the matrix is
239: ! determined by PETSc at runtime.
241: PetscCallA(MatCreate(PETSC_COMM_WORLD, A, ierr))
242: PetscCallA(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m*n, m*n, ierr))
243: PetscCallA(MatSetFromOptions(A, ierr))
244: PetscCallA(MatSetUp(A, ierr))
246: ! Currently, all PETSc parallel matrix formats are partitioned by
247: ! contiguous chunks of rows across the processors. Determine which
248: ! rows of the matrix are locally owned.
250: PetscCallA(MatGetOwnershipRange(A, Istart, Iend, ierr))
252: ! Set matrix elements for the 2-D, five-point stencil in parallel.
253: ! - Each processor needs to insert only elements that it owns
254: ! locally (but any non-local elements will be sent to the
255: ! appropriate processor during matrix assembly).
256: ! - Always specify global row and columns of matrix entries.
257: ! - Note that MatSetValues() uses 0-based row and column numbers
258: ! in Fortran as well as in C.
260: do II = Istart, Iend - 1
261: v = -1.0
262: i = II/n
263: j = II - i*n
264: if (i > 0) then
265: JJ = II - n
266: PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, [II], 1_PETSC_INT_KIND, [JJ], [v], ADD_VALUES, ierr))
267: end if
268: if (i < m - 1) then
269: JJ = II + n
270: PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, [II], 1_PETSC_INT_KIND, [JJ], [v], ADD_VALUES, ierr))
271: end if
272: if (j > 0) then
273: JJ = II - 1
274: PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, [II], 1_PETSC_INT_KIND, [JJ], [v], ADD_VALUES, ierr))
275: end if
276: if (j < n - 1) then
277: JJ = II + 1
278: PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, [II], 1_PETSC_INT_KIND, [JJ], [v], ADD_VALUES, ierr))
279: end if
280: v = 4.0
281: PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, [II], 1_PETSC_INT_KIND, [II], [v], ADD_VALUES, ierr))
282: end do
284: ! Assemble matrix, using the 2-step process:
285: ! MatAssemblyBegin(), MatAssemblyEnd()
286: ! Computations can be done while messages are in transition,
287: ! by placing code between these two statements.
289: PetscCallA(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
290: PetscCallA(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))
292: ! Create parallel vectors.
293: ! - Here, the parallel partitioning of the vector is determined by
294: ! PETSc at runtime. We could also specify the local dimensions
295: ! if desired -- or use the more general routine VecCreate().
296: ! - When solving a linear system, the vectors and matrices MUST
297: ! be partitioned accordingly. PETSc automatically generates
298: ! appropriately partitioned matrices and vectors when MatCreate()
299: ! and VecCreate() are used with the same communicator.
300: ! - Note: We form 1 vector from scratch and then duplicate as needed.
302: PetscCallA(VecCreateFromOptions(PETSC_COMM_WORLD, PETSC_NULL_CHARACTER, 1_PETSC_INT_KIND, PETSC_DECIDE, m*n, u, ierr))
303: PetscCallA(VecDuplicate(u, b, ierr))
304: PetscCallA(VecDuplicate(b, x, ierr))
306: ! Set exact solution; then compute right-hand-side vector.
308: PetscCallA(VecSet(u, one, ierr))
309: PetscCallA(MatMult(A, u, b, ierr))
311: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
312: ! Create the linear solver and set various options
313: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
315: ! Create linear solver context
317: PetscCallA(KSPCreate(PETSC_COMM_WORLD, ksp, ierr))
319: ! Set operators. Here the matrix that defines the linear system
320: ! also serves as the matrix from which the preconditioner is constructed.
322: PetscCallA(KSPSetOperators(ksp, A, A, ierr))
324: ! Set linear solver defaults for this problem (optional).
325: ! - By extracting the KSP and PC contexts from the KSP context,
326: ! we can then directly call any KSP and PC routines
327: ! to set various options.
329: PetscCallA(KSPGetPC(ksp, pc, ierr))
330: PetscCallA(KSPSetTolerances(ksp, tol, PETSC_CURRENT_REAL, PETSC_CURRENT_REAL, PETSC_CURRENT_INTEGER, ierr))
332: !
333: ! Set a user-defined shell preconditioner
334: !
336: ! (Required) Indicate to PETSc that we are using a shell preconditioner
337: PetscCallA(PCSetType(pc, PCSHELL, ierr))
339: ! (Required) Set the user-defined routine for applying the preconditioner
340: PetscCallA(PCShellSetApply(pc, SampleShellPCApply, ierr))
342: ! (Optional) Do any setup required for the preconditioner
343: ! Note: if you use PCShellSetSetUp, this will be done for your
344: PetscCallA(SampleShellPCSetUp(pc, x, ierr))
346: ! Set runtime options, e.g.,
347: ! -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
348: ! These options will override those specified above as long as
349: ! KSPSetFromOptions() is called _after_ any other customization
350: ! routines.
352: PetscCallA(KSPSetFromOptions(ksp, ierr))
354: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
355: ! Solve the linear system
356: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
358: PetscCallA(KSPSolve(ksp, b, x, ierr))
360: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
361: ! Check solution and clean up
362: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
364: ! Check the error
365: PetscCallA(VecAXPY(x, neg_one, u, ierr))
366: PetscCallA(VecNorm(x, NORM_2, norm, ierr))
367: PetscCallA(KSPGetIterationNumber(ksp, its, ierr))
369: if (rank == 0) then
370: write (6, 100) norm, its
371: end if
372: 100 format('Norm of error ', 1pe11.4, ' iterations ', i5)
374: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
375: ! Solve with a block of right-hand sides using the block
376: ! callbacks of the shell preconditioner
377: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
379: if (blocksolve) then
380: nrhs = 3
381: nits = 5
382: matapply_calls = 0
383: matapplyrichardson_calls = 0
385: PetscCallA(MatGetLocalSize(A, mloc, PETSC_NULL_INTEGER, ierr))
386: PetscCallA(MatCreateDense(PETSC_COMM_WORLD, mloc, PETSC_DECIDE, PETSC_DECIDE, nrhs, PETSC_NULL_SCALAR_ARRAY, blockb, ierr))
387: PetscCallA(MatSetRandom(blockb, PETSC_NULL_RANDOM, ierr))
388: PetscCallA(MatDuplicate(blockb, MAT_DO_NOT_COPY_VALUES, blockx, ierr))
389: PetscCallA(MatDuplicate(blockb, MAT_DO_NOT_COPY_VALUES, blocky, ierr))
390: PetscCallA(MatDuplicate(blockb, MAT_DO_NOT_COPY_VALUES, matwork, ierr))
392: ! A fixed number of Richardson iterations, so that the two block solves below perform the same operations
393: PetscCallA(KSPCreate(PETSC_COMM_WORLD, blockksp, ierr))
394: PetscCallA(KSPSetOperators(blockksp, A, A, ierr))
395: PetscCallA(KSPSetType(blockksp, KSPRICHARDSON, ierr))
396: PetscCallA(KSPSetNormType(blockksp, KSP_NORM_NONE, ierr))
397: PetscCallA(KSPSetTolerances(blockksp, PETSC_CURRENT_REAL, PETSC_CURRENT_REAL, PETSC_CURRENT_REAL, nits, ierr))
398: PetscCallA(KSPGetPC(blockksp, blockpc, ierr))
399: PetscCallA(PCSetType(blockpc, PCSHELL, ierr))
400: PetscCallA(PCShellSetApply(blockpc, SampleShellPCApply, ierr))
401: PetscCallA(PCShellSetMatApply(blockpc, SampleShellPCMatApply, ierr))
403: ! Without PCMatApplyRichardson(), the block iteration of KSPRICHARDSON calls PCMatApply() once per iteration
404: PetscCallA(KSPMatSolve(blockksp, blockb, blockx, ierr))
406: ! With PCMatApplyRichardson(), KSPRICHARDSON delegates the whole block iteration to the preconditioner
407: PetscCallA(PCShellSetMatApplyRichardson(blockpc, SampleShellPCMatApplyRichardson, ierr))
408: PetscCallA(KSPMatSolve(blockksp, blockb, blocky, ierr))
410: PetscCallA(MatAXPY(blocky, neg_one, blockx, SAME_NONZERO_PATTERN, ierr))
411: PetscCallA(MatNorm(blocky, NORM_FROBENIUS, blocknorm, ierr))
412: if (rank == 0) then
413: write (6, 120) matapply_calls, matapplyrichardson_calls
414: if (blocknorm < 1.e-12) then
415: write (6, 130)
416: else
417: write (6, 140)
418: end if
419: end if
421: PetscCallA(KSPDestroy(blockksp, ierr))
422: PetscCallA(MatDestroy(matwork, ierr))
423: PetscCallA(MatDestroy(blocky, ierr))
424: PetscCallA(MatDestroy(blockx, ierr))
425: PetscCallA(MatDestroy(blockb, ierr))
426: end if
427: 120 format('SampleShellPCMatApply() calls ', i5, ' SampleShellPCMatApplyRichardson() calls ', i5)
428: 130 format('KSPMatSolve() with and without PCMatApplyRichardson() agree')
429: 140 format('KSPMatSolve() with and without PCMatApplyRichardson() disagree')
431: ! Free work space. All PETSc objects should be destroyed when they
432: ! are no longer needed.
433: PetscCallA(KSPDestroy(ksp, ierr))
434: PetscCallA(VecDestroy(u, ierr))
435: PetscCallA(VecDestroy(x, ierr))
436: PetscCallA(VecDestroy(b, ierr))
437: PetscCallA(MatDestroy(A, ierr))
439: ! Free up PCShell data
440: PetscCallA(PCDestroy(sor, ierr))
441: PetscCallA(PCDestroy(jacobi, ierr))
442: PetscCallA(VecDestroy(work, ierr))
444: ! Always call PetscFinalize() before exiting a program.
445: PetscCallA(PetscFinalize(ierr))
446: end
448: !/*TEST
449: !
450: ! test:
451: ! requires: !single
452: !
453: ! test:
454: ! suffix: matsolve
455: ! requires: !single
456: ! args: -matsolve
457: !
458: !TEST*/