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
16: contains
17: !/***********************************************************************/
18: !/* Routines for a user-defined shell preconditioner */
19: !/***********************************************************************/
21: !
22: ! SampleShellPCSetUp - This routine sets up a user-defined
23: ! preconditioner context.
24: !
25: ! Input Parameters:
26: ! pc - preconditioner object
27: ! x - vector
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 used to compute the preconditioner; this vector is then
36: ! used within the routine SampleShellPCApply().
37: !
38: subroutine SampleShellPCSetUp(pc, x, ierr)
40: PC pc
41: Vec x
42: Mat pmat
43: PetscErrorCode ierr
45: PetscCallA(PCGetOperators(pc, PETSC_NULL_MAT, pmat, ierr))
46: PetscCallA(PCCreate(PETSC_COMM_WORLD, jacobi, ierr))
47: PetscCallA(PCSetType(jacobi, PCJACOBI, ierr))
48: PetscCallA(PCSetOperators(jacobi, pmat, pmat, ierr))
49: PetscCallA(PCSetUp(jacobi, ierr))
51: PetscCallA(PCCreate(PETSC_COMM_WORLD, sor, ierr))
52: PetscCallA(PCSetType(sor, PCSOR, ierr))
53: PetscCallA(PCSetOperators(sor, pmat, pmat, ierr))
54: ! PetscCallA(PCSORSetSymmetric(sor,SOR_LOCAL_SYMMETRIC_SWEEP,ierr))
55: PetscCallA(PCSetUp(sor, ierr))
57: PetscCallA(VecDuplicate(x, work, ierr))
58: end
60: ! -------------------------------------------------------------------
61: !
62: ! SampleShellPCApply - This routine demonstrates the use of a
63: ! user-provided preconditioner.
64: !
65: ! Input Parameters:
66: ! pc - preconditioner object
67: ! x - input vector
68: !
69: ! Output Parameters:
70: ! y - preconditioned vector
71: ! ierr - error code (nonzero if error has been detected)
72: !
73: ! Notes:
74: ! This code implements the Jacobi preconditioner plus the
75: ! SOR preconditioner
76: !
77: ! YOU CAN GET THE EXACT SAME EFFECT WITH THE PCCOMPOSITE preconditioner using
78: ! mpiexec -n 1 ex21f -ksp_monitor -pc_type composite -pc_composite_pcs jacobi,sor -pc_composite_type additive
79: !
80: subroutine SampleShellPCApply(pc, x, y, ierr)
82: PC pc
83: Vec x, y
84: PetscErrorCode ierr
85: PetscScalar, parameter :: one = 1.0
87: PetscCallA(PCApply(jacobi, x, y, ierr))
88: PetscCallA(PCApply(sor, x, work, ierr))
89: PetscCallA(VecAXPY(y, one, work, ierr))
90: end
92: end module
94: program main
95: use ex62fmodule
96: implicit none
98: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
99: ! Variable declarations
100: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101: !
102: ! Variables:
103: ! ksp - linear solver context
104: ! ksp - Krylov subspace method context
105: ! pc - preconditioner context
106: ! x, b, u - approx solution, right-hand side, exact solution vectors
107: ! A - matrix that defines linear system
108: ! its - iterations for convergence
109: ! norm - norm of solution error
111: Vec x, b, u
112: Mat A
113: PC pc
114: KSP ksp
115: PetscScalar v
116: PetscScalar, parameter :: one = 1.0, neg_one = -1.0
117: PetscReal, parameter :: tol = 1e-7
118: PetscReal norm
119: PetscInt i, j, II, JJ, Istart, Iend, its
120: PetscInt m, n
121: PetscMPIInt rank
122: PetscBool flg
123: PetscErrorCode ierr
125: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126: ! Beginning of program
127: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
129: PetscCallA(PetscInitialize(ierr))
130: m = 8
131: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-m', m, flg, ierr))
132: n = 7
133: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', n, flg, ierr))
134: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
136: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
137: ! Compute the matrix and right-hand-side vector that define
138: ! the linear system, Ax = b.
139: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
141: ! Create parallel matrix, specifying only its global dimensions.
142: ! When using MatCreate(), the matrix format can be specified at
143: ! runtime. Also, the parallel partitioning of the matrix is
144: ! determined by PETSc at runtime.
146: PetscCallA(MatCreate(PETSC_COMM_WORLD, A, ierr))
147: PetscCallA(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m*n, m*n, ierr))
148: PetscCallA(MatSetFromOptions(A, ierr))
149: PetscCallA(MatSetUp(A, ierr))
151: ! Currently, all PETSc parallel matrix formats are partitioned by
152: ! contiguous chunks of rows across the processors. Determine which
153: ! rows of the matrix are locally owned.
155: PetscCallA(MatGetOwnershipRange(A, Istart, Iend, ierr))
157: ! Set matrix elements for the 2-D, five-point stencil in parallel.
158: ! - Each processor needs to insert only elements that it owns
159: ! locally (but any non-local elements will be sent to the
160: ! appropriate processor during matrix assembly).
161: ! - Always specify global row and columns of matrix entries.
162: ! - Note that MatSetValues() uses 0-based row and column numbers
163: ! in Fortran as well as in C.
165: do II = Istart, Iend - 1
166: v = -1.0
167: i = II/n
168: j = II - i*n
169: if (i > 0) then
170: JJ = II - n
171: PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, [II], 1_PETSC_INT_KIND, [JJ], [v], ADD_VALUES, ierr))
172: end if
173: if (i < m - 1) then
174: JJ = II + n
175: PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, [II], 1_PETSC_INT_KIND, [JJ], [v], ADD_VALUES, ierr))
176: end if
177: if (j > 0) then
178: JJ = II - 1
179: PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, [II], 1_PETSC_INT_KIND, [JJ], [v], ADD_VALUES, ierr))
180: end if
181: if (j < n - 1) then
182: JJ = II + 1
183: PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, [II], 1_PETSC_INT_KIND, [JJ], [v], ADD_VALUES, ierr))
184: end if
185: v = 4.0
186: PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, [II], 1_PETSC_INT_KIND, [II], [v], ADD_VALUES, ierr))
187: end do
189: ! Assemble matrix, using the 2-step process:
190: ! MatAssemblyBegin(), MatAssemblyEnd()
191: ! Computations can be done while messages are in transition,
192: ! by placing code between these two statements.
194: PetscCallA(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
195: PetscCallA(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))
197: ! Create parallel vectors.
198: ! - Here, the parallel partitioning of the vector is determined by
199: ! PETSc at runtime. We could also specify the local dimensions
200: ! if desired -- or use the more general routine VecCreate().
201: ! - When solving a linear system, the vectors and matrices MUST
202: ! be partitioned accordingly. PETSc automatically generates
203: ! appropriately partitioned matrices and vectors when MatCreate()
204: ! and VecCreate() are used with the same communicator.
205: ! - Note: We form 1 vector from scratch and then duplicate as needed.
207: PetscCallA(VecCreateFromOptions(PETSC_COMM_WORLD, PETSC_NULL_CHARACTER, 1_PETSC_INT_KIND, PETSC_DECIDE, m*n, u, ierr))
208: PetscCallA(VecDuplicate(u, b, ierr))
209: PetscCallA(VecDuplicate(b, x, ierr))
211: ! Set exact solution; then compute right-hand-side vector.
213: PetscCallA(VecSet(u, one, ierr))
214: PetscCallA(MatMult(A, u, b, ierr))
216: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
217: ! Create the linear solver and set various options
218: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
220: ! Create linear solver context
222: PetscCallA(KSPCreate(PETSC_COMM_WORLD, ksp, ierr))
224: ! Set operators. Here the matrix that defines the linear system
225: ! also serves as the matrix from which the preconditioner is constructed.
227: PetscCallA(KSPSetOperators(ksp, A, A, ierr))
229: ! Set linear solver defaults for this problem (optional).
230: ! - By extracting the KSP and PC contexts from the KSP context,
231: ! we can then directly call any KSP and PC routines
232: ! to set various options.
234: PetscCallA(KSPGetPC(ksp, pc, ierr))
235: PetscCallA(KSPSetTolerances(ksp, tol, PETSC_CURRENT_REAL, PETSC_CURRENT_REAL, PETSC_CURRENT_INTEGER, ierr))
237: !
238: ! Set a user-defined shell preconditioner
239: !
241: ! (Required) Indicate to PETSc that we are using a shell preconditioner
242: PetscCallA(PCSetType(pc, PCSHELL, ierr))
244: ! (Required) Set the user-defined routine for applying the preconditioner
245: PetscCallA(PCShellSetApply(pc, SampleShellPCApply, ierr))
247: ! (Optional) Do any setup required for the preconditioner
248: ! Note: if you use PCShellSetSetUp, this will be done for your
249: PetscCallA(SampleShellPCSetUp(pc, x, ierr))
251: ! Set runtime options, e.g.,
252: ! -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
253: ! These options will override those specified above as long as
254: ! KSPSetFromOptions() is called _after_ any other customization
255: ! routines.
257: PetscCallA(KSPSetFromOptions(ksp, ierr))
259: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
260: ! Solve the linear system
261: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
263: PetscCallA(KSPSolve(ksp, b, x, ierr))
265: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
266: ! Check solution and clean up
267: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
269: ! Check the error
270: PetscCallA(VecAXPY(x, neg_one, u, ierr))
271: PetscCallA(VecNorm(x, NORM_2, norm, ierr))
272: PetscCallA(KSPGetIterationNumber(ksp, its, ierr))
274: if (rank == 0) then
275: if (norm > 1.e-12) then
276: write (6, 100) norm, its
277: else
278: write (6, 110) its
279: end if
280: end if
281: 100 format('Norm of error ', 1pe11.4, ' iterations ', i5)
282: 110 format('Norm of error < 1.e-12,iterations ', i5)
284: ! Free work space. All PETSc objects should be destroyed when they
285: ! are no longer needed.
286: PetscCallA(KSPDestroy(ksp, ierr))
287: PetscCallA(VecDestroy(u, ierr))
288: PetscCallA(VecDestroy(x, ierr))
289: PetscCallA(VecDestroy(b, ierr))
290: PetscCallA(MatDestroy(A, ierr))
292: ! Free up PCShell data
293: PetscCallA(PCDestroy(sor, ierr))
294: PetscCallA(PCDestroy(jacobi, ierr))
295: PetscCallA(VecDestroy(work, ierr))
297: ! Always call PetscFinalize() before exiting a program.
298: PetscCallA(PetscFinalize(ierr))
299: end
301: !/*TEST
302: !
303: ! test:
304: ! requires: !single
305: !
306: !TEST*/