Actual source code: ex15f.F
petsc-3.3-p7 2013-05-11
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: ! Program usage: mpiexec ex15f [-help] [all PETSc options]
7: !
8: !/*T
9: ! Concepts: KSP^basic parallel example
10: ! Concepts: PC^setting a user-defined shell preconditioner
11: ! Processors: n
12: !T*/
13: !
14: ! -------------------------------------------------------------------------
16: program main
17: implicit none
19: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
20: ! Include files
21: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
22: !
23: ! petscsys.h - base PETSc routines petscvec.h - vectors
24: ! petscmat.h - matrices
25: ! petscksp.h - Krylov subspace methods petscpc.h - preconditioners
27: #include <finclude/petscsys.h>
28: #include <finclude/petscvec.h>
29: #include <finclude/petscmat.h>
30: #include <finclude/petscpc.h>
31: #include <finclude/petscksp.h>
33: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
34: ! Variable declarations
35: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
36: !
37: ! Variables:
38: ! ksp - linear solver context
39: ! ksp - Krylov subspace method context
40: ! pc - preconditioner context
41: ! x, b, u - approx solution, right-hand-side, exact solution vectors
42: ! A - matrix that defines linear system
43: ! its - iterations for convergence
44: ! norm - norm of solution error
46: Vec x,b,u
47: Mat A
48: PC pc
49: KSP ksp
50: PetscScalar v,one,neg_one
51: double precision norm,tol
52: PetscErrorCode ierr
53: PetscInt i,j,II,JJ,Istart
54: PetscInt Iend,m,n,i1,its,five
55: PetscMPIInt rank
56: PetscBool user_defined_pc,flg
58: ! Note: Any user-defined Fortran routines MUST be declared as external.
60: external SampleShellPCSetUp, SampleShellPCApply
61: external SampleShellPCDestroy
63: ! Common block to store data for user-provided preconditioner
64: common /myshellpc/ diag
65: Vec diag
67: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
68: ! Beginning of program
69: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
71: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
72: one = 1.0
73: neg_one = -1.0
74: i1 = 1
75: m = 8
76: n = 7
77: five = 5
78: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-m',m,flg,ierr)
79: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-n',n,flg,ierr)
80: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
82: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
83: ! Compute the matrix and right-hand-side vector that define
84: ! the linear system, Ax = b.
85: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
87: ! Create parallel matrix, specifying only its global dimensions.
88: ! When using MatCreate(), the matrix format can be specified at
89: ! runtime. Also, the parallel partitioning of the matrix is
90: ! determined by PETSc at runtime.
92: call MatCreate(PETSC_COMM_WORLD,A,ierr)
93: call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n,ierr)
94: call MatSetType(A, MATAIJ,ierr)
95: call MatSetFromOptions(A,ierr)
96: call MatMPIAIJSetPreallocation(A,five,PETSC_NULL_INTEGER,five, &
97: & PETSC_NULL_INTEGER,ierr)
98: call MatSeqAIJSetPreallocation(A,five,PETSC_NULL_INTEGER,ierr)
100: ! Currently, all PETSc parallel matrix formats are partitioned by
101: ! contiguous chunks of rows across the processors. Determine which
102: ! rows of the matrix are locally owned.
104: call MatGetOwnershipRange(A,Istart,Iend,ierr)
106: ! Set matrix elements for the 2-D, five-point stencil in parallel.
107: ! - Each processor needs to insert only elements that it owns
108: ! locally (but any non-local elements will be sent to the
109: ! appropriate processor during matrix assembly).
110: ! - Always specify global row and columns of matrix entries.
111: ! - Note that MatSetValues() uses 0-based row and column numbers
112: ! in Fortran as well as in C.
114: do 10, II=Istart,Iend-1
115: v = -1.0
116: i = II/n
117: j = II - i*n
118: if (i.gt.0) then
119: JJ = II - n
120: call MatSetValues(A,i1,II,i1,JJ,v,ADD_VALUES,ierr)
121: endif
122: if (i.lt.m-1) then
123: JJ = II + n
124: call MatSetValues(A,i1,II,i1,JJ,v,ADD_VALUES,ierr)
125: endif
126: if (j.gt.0) then
127: JJ = II - 1
128: call MatSetValues(A,i1,II,i1,JJ,v,ADD_VALUES,ierr)
129: endif
130: if (j.lt.n-1) then
131: JJ = II + 1
132: call MatSetValues(A,i1,II,i1,JJ,v,ADD_VALUES,ierr)
133: endif
134: v = 4.0
135: call MatSetValues(A,i1,II,i1,II,v,ADD_VALUES,ierr)
136: 10 continue
138: ! Assemble matrix, using the 2-step process:
139: ! MatAssemblyBegin(), MatAssemblyEnd()
140: ! Computations can be done while messages are in transition,
141: ! by placing code between these two statements.
143: call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
144: call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
146: ! Create parallel vectors.
147: ! - Here, the parallel partitioning of the vector is determined by
148: ! PETSc at runtime. We could also specify the local dimensions
149: ! if desired -- or use the more general routine VecCreate().
150: ! - When solving a linear system, the vectors and matrices MUST
151: ! be partitioned accordingly. PETSc automatically generates
152: ! appropriately partitioned matrices and vectors when MatCreate()
153: ! and VecCreate() are used with the same communicator.
154: ! - Note: We form 1 vector from scratch and then duplicate as needed.
156: call VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,m*n,u,ierr)
157: call VecDuplicate(u,b,ierr)
158: call VecDuplicate(b,x,ierr)
160: ! Set exact solution; then compute right-hand-side vector.
162: call VecSet(u,one,ierr)
163: call MatMult(A,u,b,ierr)
165: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
166: ! Create the linear solver and set various options
167: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
169: ! Create linear solver context
171: call KSPCreate(PETSC_COMM_WORLD,ksp,ierr)
173: ! Set operators. Here the matrix that defines the linear system
174: ! also serves as the preconditioning matrix.
176: call KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN,ierr)
178: ! Set linear solver defaults for this problem (optional).
179: ! - By extracting the KSP and PC contexts from the KSP context,
180: ! we can then directly directly call any KSP and PC routines
181: ! to set various options.
183: call KSPGetPC(ksp,pc,ierr)
184: tol = 1.e-7
185: call KSPSetTolerances(ksp,tol,PETSC_DEFAULT_DOUBLE_PRECISION, &
186: & PETSC_DEFAULT_DOUBLE_PRECISION,PETSC_DEFAULT_INTEGER,ierr)
188: !
189: ! Set a user-defined shell preconditioner if desired
190: !
191: call PetscOptionsHasName(PETSC_NULL_CHARACTER,'-user_defined_pc', &
192: & user_defined_pc,ierr)
194: if (user_defined_pc) then
196: ! (Required) Indicate to PETSc that we are using a shell preconditioner
197: call PCSetType(pc,PCSHELL,ierr)
199: ! (Required) Set the user-defined routine for applying the preconditioner
200: call PCShellSetApply(pc,SampleShellPCApply,ierr)
202: ! (Optional) Do any setup required for the preconditioner
203: call PCShellSetSetUp(pc,SampleShellPCSetUp,ierr)
205: ! (Optional) Frees any objects we created for the preconditioner
206: call PCShellSetDestroy(pc,SampleShellPCDestroy,ierr)
208: else
209: call PCSetType(pc,PCJACOBI,ierr)
210: endif
212: ! Set runtime options, e.g.,
213: ! -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
214: ! These options will override those specified above as long as
215: ! KSPSetFromOptions() is called _after_ any other customization
216: ! routines.
218: call KSPSetFromOptions(ksp,ierr)
220: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
221: ! Solve the linear system
222: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
224: call KSPSolve(ksp,b,x,ierr)
226: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
227: ! Check solution and clean up
228: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
230: ! Check the error
232: call VecAXPY(x,neg_one,u,ierr)
233: call VecNorm(x,NORM_2,norm,ierr)
234: call KSPGetIterationNumber(ksp,its,ierr)
236: if (rank .eq. 0) then
237: if (norm .gt. 1.e-12) then
238: write(6,100) norm,its
239: else
240: write(6,110) its
241: endif
242: endif
243: 100 format('Norm of error ',1pe11.4,' iterations ',i5)
244: 110 format('Norm of error < 1.e-12,iterations ',i5)
246: ! Free work space. All PETSc objects should be destroyed when they
247: ! are no longer needed.
249: call KSPDestroy(ksp,ierr)
250: call VecDestroy(u,ierr)
251: call VecDestroy(x,ierr)
252: call VecDestroy(b,ierr)
253: call MatDestroy(A,ierr)
255: ! Always call PetscFinalize() before exiting a program.
257: call PetscFinalize(ierr)
258: end
260: !/***********************************************************************/
261: !/* Routines for a user-defined shell preconditioner */
262: !/***********************************************************************/
264: !
265: ! SampleShellPCSetUp - This routine sets up a user-defined
266: ! preconditioner context.
267: !
268: ! Input Parameters:
269: ! pc - preconditioner object
270: !
271: ! Output Parameter:
272: ! ierr - error code (nonzero if error has been detected)
273: !
274: ! Notes:
275: ! In this example, we define the shell preconditioner to be Jacobi
276: ! method. Thus, here we create a work vector for storing the reciprocal
277: ! of the diagonal of the preconditioner matrix; this vector is then
278: ! used within the routine SampleShellPCApply().
279: !
280: subroutine SampleShellPCSetUp(pc,ierr)
282: implicit none
284: #include <finclude/petscsys.h>
285: #include <finclude/petscvec.h>
286: #include <finclude/petscmat.h>
287: PC pc
289: Mat pmat
290: integer ierr
292: ! Common block to store data for user-provided preconditioner
293: ! Normally we would recommend storing all the work data (like diag) in
294: ! the context set with PCShellSetContext()
296: common /myshellpc/ diag
297: Vec diag
299: call PCGetOperators(pc,PETSC_NULL_OBJECT,pmat,PETSC_NULL_INTEGER, &
300: & ierr)
301: call MatGetVecs(pmat,diag,PETSC_NULL_OBJECT,ierr)
302: call MatGetDiagonal(pmat,diag,ierr)
303: call VecReciprocal(diag,ierr)
305: end
307: ! -------------------------------------------------------------------
308: !
309: ! SampleShellPCApply - This routine demonstrates the use of a
310: ! user-provided preconditioner.
311: !
312: ! Input Parameters:
313: ! pc - preconditioner object
314: ! x - input vector
315: !
316: ! Output Parameters:
317: ! y - preconditioned vector
318: ! ierr - error code (nonzero if error has been detected)
319: !
320: ! Notes:
321: ! This code implements the Jacobi preconditioner, merely as an
322: ! example of working with a PCSHELL. Note that the Jacobi method
323: ! is already provided within PETSc.
324: !
325: subroutine SampleShellPCApply(pc,x,y,ierr)
327: implicit none
329: #include <finclude/petscsys.h>
330: #include <finclude/petscvec.h>
332: PC pc
333: Vec x,y
334: integer ierr
336: ! Common block to store data for user-provided preconditioner
337: common /myshellpc/ diag
338: Vec diag
340: call VecPointwiseMult(y,x,diag,ierr)
342: end
344: !/***********************************************************************/
345: !/* Routines for a user-defined shell preconditioner */
346: !/***********************************************************************/
348: !
349: ! SampleShellPCDestroy - This routine destroys (frees the memory of) any
350: ! objects we made for the preconditioner
351: !
352: ! Input Parameters:
353: ! pc - for this example we use the actual PC as our shell context
354: !
355: ! Output Parameter:
356: ! ierr - error code (nonzero if error has been detected)
357: !
359: subroutine SampleShellPCDestroy(pc,ierr)
361: implicit none
363: #include <finclude/petscsys.h>
364: #include <finclude/petscvec.h>
365: #include <finclude/petscmat.h>
366: PC pc
367: PetscErrorCode ierr
369: ! Common block to store data for user-provided preconditioner
370: ! Normally we would recommend storing all the work data (like diag) in
371: ! the context set with PCShellSetContext()
373: common /myshellpc/ diag
374: Vec diag
376: call VecDestroy(diag,ierr)
378: end