Actual source code: ex15f.F90
1: !
2: ! Solves a linear system in parallel with KSP. Also indicates
3: ! use of a user-provided preconditioner. Input parameters include:
4: ! -user_defined_pc : Activate a user-defined preconditioner
5: !
6: ! -------------------------------------------------------------------------
7: !
8: ! Module contains diag needed by shell preconditioner
9: !
10: #include <petsc/finclude/petscksp.h>
11: module ex15fmodule
12: use petscksp
13: implicit none
14: Vec diag
16: contains
18: !/***********************************************************************/
19: !/* Routines for a user-defined shell preconditioner */
20: !/***********************************************************************/
22: !
23: ! SampleShellPCSetUp - This routine sets up a user-defined
24: ! preconditioner context.
25: !
26: ! Input Parameters:
27: ! pc - preconditioner object
28: !
29: ! Output Parameter:
30: ! ierr - error code (nonzero if error has been detected)
31: !
32: ! Notes:
33: ! In this example, we define the shell preconditioner to be Jacobi
34: ! method. Thus, here we create a work vector for storing the reciprocal
35: ! of the diagonal of the matrix; this vector is then
36: ! used within the routine SampleShellPCApply().
37: !
38: subroutine SampleShellPCSetUp(pc, ierr)
40: PC pc
41: Mat pmat
42: PetscErrorCode ierr
44: PetscCallA(PCGetOperators(pc, PETSC_NULL_MAT, pmat, ierr))
45: PetscCallA(MatCreateVecs(pmat, diag, PETSC_NULL_VEC, ierr))
46: PetscCallA(MatGetDiagonal(pmat, diag, ierr))
47: PetscCallA(VecReciprocal(diag, ierr))
49: end
51: ! -------------------------------------------------------------------
52: !
53: ! SampleShellPCApply - This routine demonstrates the use of a
54: ! user-provided preconditioner.
55: !
56: ! Input Parameters:
57: ! pc - preconditioner object
58: ! x - input vector
59: !
60: ! Output Parameters:
61: ! y - preconditioned vector
62: ! ierr - error code (nonzero if error has been detected)
63: !
64: ! Notes:
65: ! This code implements the Jacobi preconditioner, merely as an
66: ! example of working with a PCSHELL. Note that the Jacobi method
67: ! is already provided within PETSc.
68: !
69: subroutine SampleShellPCApply(pc, x, y, ierr)
71: PC pc
72: Vec x, y
73: PetscErrorCode ierr
75: PetscCallA(VecPointwiseMult(y, x, diag, ierr))
77: end
79: !/***********************************************************************/
80: !/* Routines for a user-defined shell preconditioner */
81: !/***********************************************************************/
83: !
84: ! SampleShellPCDestroy - This routine destroys (frees the memory of) any
85: ! objects we made for the preconditioner
86: !
87: ! Input Parameters:
88: ! pc - for this example we use the actual PC as our shell context
89: !
90: ! Output Parameter:
91: ! ierr - error code (nonzero if error has been detected)
92: !
94: subroutine SampleShellPCDestroy(pc, ierr)
96: PC pc
97: PetscErrorCode ierr
99: ! Normally we would recommend storing all the work data (like diag) in
100: ! the context set with PCShellSetContext()
102: PetscCallA(VecDestroy(diag, ierr))
104: end
106: end module
108: program main
109: use ex15fmodule
110: implicit none
112: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
113: ! Variable declarations
114: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
115: !
116: ! Variables:
117: ! ksp - linear solver context
118: ! ksp - Krylov subspace method context
119: ! pc - preconditioner context
120: ! x, b, u - approx solution, right-hand side, exact solution vectors
121: ! A - matrix that defines linear system
122: ! its - iterations for convergence
123: ! norm - norm of solution error
125: Vec x, b, u
126: Mat A
127: PC pc
128: KSP ksp
129: PetscScalar v, one, neg_one
130: PetscReal norm, tol
131: PetscErrorCode ierr
132: PetscInt i, j, II, JJ, Istart
133: PetscInt Iend, m, n, i1, its, five
134: PetscMPIInt rank
135: PetscBool user_defined_pc, flg
137: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138: ! Beginning of program
139: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
141: PetscCallA(PetscInitialize(ierr))
142: one = 1.0
143: neg_one = -1.0
144: i1 = 1
145: m = 8
146: n = 7
147: five = 5
148: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-m', m, flg, ierr))
149: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', n, flg, ierr))
150: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
152: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
153: ! Compute the matrix and right-hand-side vector that define
154: ! the linear system, Ax = b.
155: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
157: ! Create parallel matrix, specifying only its global dimensions.
158: ! When using MatCreate(), the matrix format can be specified at
159: ! runtime. Also, the parallel partitioning of the matrix is
160: ! determined by PETSc at runtime.
162: PetscCallA(MatCreate(PETSC_COMM_WORLD, A, ierr))
163: PetscCallA(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m*n, m*n, ierr))
164: PetscCallA(MatSetType(A, MATAIJ, ierr))
165: PetscCallA(MatSetFromOptions(A, ierr))
166: PetscCallA(MatMPIAIJSetPreallocation(A, five, PETSC_NULL_INTEGER_ARRAY, five, PETSC_NULL_INTEGER_ARRAY, ierr))
167: PetscCallA(MatSeqAIJSetPreallocation(A, five, PETSC_NULL_INTEGER_ARRAY, ierr))
169: ! Currently, all PETSc parallel matrix formats are partitioned by
170: ! contiguous chunks of rows across the processors. Determine which
171: ! rows of the matrix are locally owned.
173: PetscCallA(MatGetOwnershipRange(A, Istart, Iend, ierr))
175: ! Set matrix elements for the 2-D, five-point stencil in parallel.
176: ! - Each processor needs to insert only elements that it owns
177: ! locally (but any non-local elements will be sent to the
178: ! appropriate processor during matrix assembly).
179: ! - Always specify global row and columns of matrix entries.
180: ! - Note that MatSetValues() uses 0-based row and column numbers
181: ! in Fortran as well as in C.
183: do II = Istart, Iend - 1
184: v = -1.0
185: i = II/n
186: j = II - i*n
187: if (i > 0) then
188: JJ = II - n
189: PetscCallA(MatSetValues(A, i1, [II], i1, [JJ], [v], ADD_VALUES, ierr))
190: end if
191: if (i < m - 1) then
192: JJ = II + n
193: PetscCallA(MatSetValues(A, i1, [II], i1, [JJ], [v], ADD_VALUES, ierr))
194: end if
195: if (j > 0) then
196: JJ = II - 1
197: PetscCallA(MatSetValues(A, i1, [II], i1, [JJ], [v], ADD_VALUES, ierr))
198: end if
199: if (j < n - 1) then
200: JJ = II + 1
201: PetscCallA(MatSetValues(A, i1, [II], i1, [JJ], [v], ADD_VALUES, ierr))
202: end if
203: v = 4.0
204: PetscCallA(MatSetValues(A, i1, [II], i1, [II], [v], ADD_VALUES, ierr))
205: end do
207: ! Assemble matrix, using the 2-step process:
208: ! MatAssemblyBegin(), MatAssemblyEnd()
209: ! Computations can be done while messages are in transition,
210: ! by placing code between these two statements.
212: PetscCallA(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
213: PetscCallA(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))
215: ! Create parallel vectors.
216: ! - Here, the parallel partitioning of the vector is determined by
217: ! PETSc at runtime. We could also specify the local dimensions
218: ! if desired -- or use the more general routine VecCreate().
219: ! - When solving a linear system, the vectors and matrices MUST
220: ! be partitioned accordingly. PETSc automatically generates
221: ! appropriately partitioned matrices and vectors when MatCreate()
222: ! and VecCreate() are used with the same communicator.
223: ! - Note: We form 1 vector from scratch and then duplicate as needed.
225: PetscCallA(VecCreateFromOptions(PETSC_COMM_WORLD, PETSC_NULL_CHARACTER, i1, PETSC_DECIDE, m*n, u, ierr))
226: PetscCallA(VecDuplicate(u, b, ierr))
227: PetscCallA(VecDuplicate(b, x, ierr))
229: ! Set exact solution; then compute right-hand-side vector.
231: PetscCallA(VecSet(u, one, ierr))
232: PetscCallA(MatMult(A, u, b, ierr))
234: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
235: ! Create the linear solver and set various options
236: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
238: ! Create linear solver context
240: PetscCallA(KSPCreate(PETSC_COMM_WORLD, ksp, ierr))
242: ! Set operators. Here the matrix that defines the linear system
243: ! also serves as the matrix from which the preconditioner is constructed.
245: PetscCallA(KSPSetOperators(ksp, A, A, ierr))
247: ! Set linear solver defaults for this problem (optional).
248: ! - By extracting the KSP and PC contexts from the KSP context,
249: ! we can then directly call any KSP and PC routines
250: ! to set various options.
252: PetscCallA(KSPGetPC(ksp, pc, ierr))
253: tol = 1.e-7
254: PetscCallA(KSPSetTolerances(ksp, tol, PETSC_CURRENT_REAL, PETSC_CURRENT_REAL, PETSC_CURRENT_INTEGER, ierr))
256: !
257: ! Set a user-defined shell preconditioner if desired
258: !
259: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-user_defined_pc', user_defined_pc, ierr))
261: if (user_defined_pc) then
263: ! (Required) Indicate to PETSc that we are using a shell preconditioner
264: PetscCallA(PCSetType(pc, PCSHELL, ierr))
266: ! (Required) Set the user-defined routine for applying the preconditioner
267: PetscCallA(PCShellSetApply(pc, SampleShellPCApply, ierr))
269: ! (Optional) Do any setup required for the preconditioner
270: PetscCallA(PCShellSetSetUp(pc, SampleShellPCSetUp, ierr))
272: ! (Optional) Frees any objects we created for the preconditioner
273: PetscCallA(PCShellSetDestroy(pc, SampleShellPCDestroy, ierr))
275: else
276: PetscCallA(PCSetType(pc, PCJACOBI, ierr))
277: end if
279: ! Set runtime options, e.g.,
280: ! -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
281: ! These options will override those specified above as long as
282: ! KSPSetFromOptions() is called _after_ any other customization
283: ! routines.
285: PetscCallA(KSPSetFromOptions(ksp, ierr))
287: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
288: ! Solve the linear system
289: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
291: PetscCallA(KSPSolve(ksp, b, x, ierr))
293: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
294: ! Check solution and clean up
295: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
297: ! Check the error
299: PetscCallA(VecAXPY(x, neg_one, u, ierr))
300: PetscCallA(VecNorm(x, NORM_2, norm, ierr))
301: PetscCallA(KSPGetIterationNumber(ksp, its, ierr))
303: if (rank == 0) then
304: if (norm > 1.e-12) then
305: write (6, 100) norm, its
306: else
307: write (6, 110) its
308: end if
309: end if
310: 100 format('Norm of error ', 1pe11.4, ' iterations ', i5)
311: 110 format('Norm of error < 1.e-12,iterations ', i5)
313: ! Free work space. All PETSc objects should be destroyed when they
314: ! are no longer needed.
316: PetscCallA(KSPDestroy(ksp, ierr))
317: PetscCallA(VecDestroy(u, ierr))
318: PetscCallA(VecDestroy(x, ierr))
319: PetscCallA(VecDestroy(b, ierr))
320: PetscCallA(MatDestroy(A, ierr))
322: ! Always call PetscFinalize() before exiting a program.
324: PetscCallA(PetscFinalize(ierr))
325: end program
326: !
327: !/*TEST
328: !
329: ! test:
330: ! nsize: 2
331: ! args: -ksp_view -user_defined_pc -ksp_gmres_cgs_refinement_type refine_always
332: !
333: !TEST*/