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: module mymoduleex21f
10: #include <petsc/finclude/petscksp.h>
11: use petscksp
12: PC jacobi,sor
13: Vec work
14: end module
16: program main
17: use mymoduleex21f
18: implicit none
20: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
21: ! Variable declarations
22: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
23: !
24: ! Variables:
25: ! ksp - linear solver context
26: ! ksp - Krylov subspace method context
27: ! pc - preconditioner context
28: ! x, b, u - approx solution, right-hand-side, exact solution vectors
29: ! A - matrix that defines linear system
30: ! its - iterations for convergence
31: ! norm - norm of solution error
33: Vec x,b,u
34: Mat A
35: PC pc
36: KSP ksp
37: PetscScalar v,one,neg_one
38: PetscReal norm,tol
39: PetscInt i,j,II,JJ,Istart
40: PetscInt Iend,m,n,its,ione
41: PetscMPIInt rank
42: PetscBool flg
43: PetscErrorCode ierr
45: ! Note: Any user-defined Fortran routines MUST be declared as external.
47: external SampleShellPCSetUp,SampleShellPCApply
49: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
50: ! Beginning of program
51: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
53: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
54: if (ierr .ne. 0) then
55: print*,'Unable to initialize PETSc'
56: stop
57: endif
58: one = 1.0
59: neg_one = -1.0
60: m = 8
61: n = 7
62: ione = 1
63: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-m',m,flg,ierr)
64: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n',n,flg,ierr)
65: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
67: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
68: ! Compute the matrix and right-hand-side vector that define
69: ! the linear system, Ax = b.
70: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
72: ! Create parallel matrix, specifying only its global dimensions.
73: ! When using MatCreate(), the matrix format can be specified at
74: ! runtime. Also, the parallel partitioning of the matrix is
75: ! determined by PETSc at runtime.
77: call MatCreate(PETSC_COMM_WORLD,A,ierr)
78: call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n,ierr)
79: call MatSetFromOptions(A,ierr)
80: call MatSetUp(A,ierr)
82: ! Currently, all PETSc parallel matrix formats are partitioned by
83: ! contiguous chunks of rows across the processors. Determine which
84: ! rows of the matrix are locally owned.
86: call MatGetOwnershipRange(A,Istart,Iend,ierr)
88: ! Set matrix elements for the 2-D, five-point stencil in parallel.
89: ! - Each processor needs to insert only elements that it owns
90: ! locally (but any non-local elements will be sent to the
91: ! appropriate processor during matrix assembly).
92: ! - Always specify global row and columns of matrix entries.
93: ! - Note that MatSetValues() uses 0-based row and column numbers
94: ! in Fortran as well as in C.
96: do 10, II=Istart,Iend-1
97: v = -1.0
98: i = II/n
99: j = II - i*n
100: if (i.gt.0) then
101: JJ = II - n
102: call MatSetValues(A,ione,II,ione,JJ,v,ADD_VALUES,ierr)
103: endif
104: if (i.lt.m-1) then
105: JJ = II + n
106: call MatSetValues(A,ione,II,ione,JJ,v,ADD_VALUES,ierr)
107: endif
108: if (j.gt.0) then
109: JJ = II - 1
110: call MatSetValues(A,ione,II,ione,JJ,v,ADD_VALUES,ierr)
111: endif
112: if (j.lt.n-1) then
113: JJ = II + 1
114: call MatSetValues(A,ione,II,ione,JJ,v,ADD_VALUES,ierr)
115: endif
116: v = 4.0
117: call MatSetValues(A,ione,II,ione,II,v,ADD_VALUES,ierr)
118: 10 continue
120: ! Assemble matrix, using the 2-step process:
121: ! MatAssemblyBegin(), MatAssemblyEnd()
122: ! Computations can be done while messages are in transition,
123: ! by placing code between these two statements.
125: call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
126: call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
128: ! Create parallel vectors.
129: ! - Here, the parallel partitioning of the vector is determined by
130: ! PETSc at runtime. We could also specify the local dimensions
131: ! if desired -- or use the more general routine VecCreate().
132: ! - When solving a linear system, the vectors and matrices MUST
133: ! be partitioned accordingly. PETSc automatically generates
134: ! appropriately partitioned matrices and vectors when MatCreate()
135: ! and VecCreate() are used with the same communicator.
136: ! - Note: We form 1 vector from scratch and then duplicate as needed.
138: call VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,m*n,u,ierr)
139: call VecDuplicate(u,b,ierr)
140: call VecDuplicate(b,x,ierr)
142: ! Set exact solution; then compute right-hand-side vector.
144: call VecSet(u,one,ierr)
145: call MatMult(A,u,b,ierr)
147: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
148: ! Create the linear solver and set various options
149: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151: ! Create linear solver context
153: call KSPCreate(PETSC_COMM_WORLD,ksp,ierr)
155: ! Set operators. Here the matrix that defines the linear system
156: ! also serves as the preconditioning matrix.
158: call KSPSetOperators(ksp,A,A,ierr)
160: ! Set linear solver defaults for this problem (optional).
161: ! - By extracting the KSP and PC contexts from the KSP context,
162: ! we can then directly directly call any KSP and PC routines
163: ! to set various options.
165: call KSPGetPC(ksp,pc,ierr)
166: tol = 1.e-7
167: call KSPSetTolerances(ksp,tol,PETSC_DEFAULT_REAL,PETSC_DEFAULT_REAL,PETSC_DEFAULT_INTEGER,ierr)
169: !
170: ! Set a user-defined shell preconditioner
171: !
173: ! (Required) Indicate to PETSc that we are using a shell preconditioner
174: call PCSetType(pc,PCSHELL,ierr)
176: ! (Required) Set the user-defined routine for applying the preconditioner
177: call PCShellSetApply(pc,SampleShellPCApply,ierr)
179: ! (Optional) Do any setup required for the preconditioner
180: ! Note: if you use PCShellSetSetUp, this will be done for your
181: call SampleShellPCSetUp(pc,x,ierr)
183: ! Set runtime options, e.g.,
184: ! -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
185: ! These options will override those specified above as long as
186: ! KSPSetFromOptions() is called _after_ any other customization
187: ! routines.
189: call KSPSetFromOptions(ksp,ierr)
191: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
192: ! Solve the linear system
193: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
195: call KSPSolve(ksp,b,x,ierr)
197: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
198: ! Check solution and clean up
199: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
201: ! Check the error
203: call VecAXPY(x,neg_one,u,ierr)
204: call VecNorm(x,NORM_2,norm,ierr)
205: call KSPGetIterationNumber(ksp,its,ierr)
207: if (rank .eq. 0) then
208: if (norm .gt. 1.e-12) then
209: write(6,100) norm,its
210: else
211: write(6,110) its
212: endif
213: endif
214: 100 format('Norm of error ',1pe11.4,' iterations ',i5)
215: 110 format('Norm of error < 1.e-12,iterations ',i5)
217: ! Free work space. All PETSc objects should be destroyed when they
218: ! are no longer needed.
220: call KSPDestroy(ksp,ierr)
221: call VecDestroy(u,ierr)
222: call VecDestroy(x,ierr)
223: call VecDestroy(b,ierr)
224: call MatDestroy(A,ierr)
226: ! Free up PCShell data
227: call PCDestroy(sor,ierr)
228: call PCDestroy(jacobi,ierr)
229: call VecDestroy(work,ierr)
231: ! Always call PetscFinalize() before exiting a program.
233: call PetscFinalize(ierr)
234: end
236: !/***********************************************************************/
237: !/* Routines for a user-defined shell preconditioner */
238: !/***********************************************************************/
240: !
241: ! SampleShellPCSetUp - This routine sets up a user-defined
242: ! preconditioner context.
243: !
244: ! Input Parameters:
245: ! pc - preconditioner object
246: ! x - vector
247: !
248: ! Output Parameter:
249: ! ierr - error code (nonzero if error has been detected)
250: !
251: ! Notes:
252: ! In this example, we define the shell preconditioner to be Jacobi
253: ! method. Thus, here we create a work vector for storing the reciprocal
254: ! of the diagonal of the preconditioner matrix; this vector is then
255: ! used within the routine SampleShellPCApply().
256: !
257: subroutine SampleShellPCSetUp(pc,x,ierr)
258: use mymoduleex21f
259: implicit none
261: PC pc
262: Vec x
263: Mat pmat
264: PetscErrorCode ierr
266: call PCGetOperators(pc,PETSC_NULL_MAT,pmat,ierr)
267: call PCCreate(PETSC_COMM_WORLD,jacobi,ierr)
268: call PCSetType(jacobi,PCJACOBI,ierr)
269: call PCSetOperators(jacobi,pmat,pmat,ierr)
270: call PCSetUp(jacobi,ierr)
272: call PCCreate(PETSC_COMM_WORLD,sor,ierr)
273: call PCSetType(sor,PCSOR,ierr)
274: call PCSetOperators(sor,pmat,pmat,ierr)
275: ! call PCSORSetSymmetric(sor,SOR_LOCAL_SYMMETRIC_SWEEP,ierr)
276: call PCSetUp(sor,ierr)
278: call VecDuplicate(x,work,ierr)
280: end
282: ! -------------------------------------------------------------------
283: !
284: ! SampleShellPCApply - This routine demonstrates the use of a
285: ! user-provided preconditioner.
286: !
287: ! Input Parameters:
288: ! pc - preconditioner object
289: ! x - input vector
290: !
291: ! Output Parameters:
292: ! y - preconditioned vector
293: ! ierr - error code (nonzero if error has been detected)
294: !
295: ! Notes:
296: ! This code implements the Jacobi preconditioner plus the
297: ! SOR preconditioner
298: !
299: ! YOU CAN GET THE EXACT SAME EFFECT WITH THE PCCOMPOSITE preconditioner using
300: ! mpiexec -n 1 ex21f -ksp_monitor -pc_type composite -pc_composite_pcs jacobi,sor -pc_composite_type additive
301: !
302: subroutine SampleShellPCApply(pc,x,y,ierr)
303: use mymoduleex21f
304: implicit none
306: PC pc
307: Vec x,y
308: PetscErrorCode ierr
309: PetscScalar one
311: one = 1.0
312: call PCApply(jacobi,x,y,ierr)
313: call PCApply(sor,x,work,ierr)
314: call VecAXPY(y,one,work,ierr)
316: end
318: !/*TEST
319: !
320: ! test:
321: ! requires: !single
322: !
323: !TEST*/