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