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