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