Actual source code: ex5f.F90
1: !Solves two linear systems in parallel with KSP. The code
2: !illustrates repeated solution of linear systems with the same preconditioner
3: !method but different matrices (having the same nonzero structure). The code
4: !also uses multiple profiling stages. Input arguments are
5: ! -m <size> : problem size
6: ! -mat_nonsym : use nonsymmetric matrix (default is symmetric)
7: #include <petsc/finclude/petscksp.h>
8: program main
9: use petscksp
11: implicit none
12: KSP :: ksp ! linear solver context
13: Mat :: C, Ctmp ! matrix
14: Vec :: x, u, b ! approx solution, RHS, exact solution
15: PetscReal :: norm, bnorm ! norm of solution residual
16: PetscScalar :: v
17: PetscScalar, parameter :: myNone = -1.0
18: PetscInt :: Ii, JJ, ldim, low, high, iglobal, Istart, Iend
19: PetscErrorCode :: ierr
20: PetscInt :: i, j, its, n, m, orthog
21: PetscMPIInt :: size, rank
22: PetscBool :: &
23: testnewC = PETSC_FALSE, &
24: testscaledMat = PETSC_FALSE, &
25: mat_nonsymmetric = PETSC_FALSE
26: PetscBool :: flg
27: PetscRandom :: rctx
28: PetscLogStage, dimension(0:1) :: stages
29: character(len=PETSC_MAX_PATH_LEN) :: outputString
31: PetscCallA(PetscInitialize(ierr))
33: orthog = 0
34: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-orthog', orthog, flg, ierr))
35: m = 3
36: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-m', m, flg, ierr))
37: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
38: PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, size, ierr))
39: n = 2*size
41: ! Set flag if we are doing a nonsymmetric problem; the default is symmetric.
43: PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-mat_nonsym', mat_nonsymmetric, flg, ierr))
44: PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-test_scaledMat', testscaledMat, flg, ierr))
46: ! Register two stages for separate profiling of the two linear solves.
47: ! Use the runtime option -log_view for a printout of performance
48: ! statistics at the program's conlusion.
50: PetscCallA(PetscLogStageRegister('Original Solve', stages(0), ierr))
51: PetscCallA(PetscLogStageRegister('Second Solve', stages(1), ierr))
53: ! -------------- Stage 0: Solve Original System ----------------------
54: ! Indicate to PETSc profiling that we're beginning the first stage
56: PetscCallA(PetscLogStagePush(stages(0), ierr))
58: ! Create parallel matrix, specifying only its global dimensions.
59: ! When using MatCreate(), the matrix format can be specified at
60: ! runtime. Also, the parallel partitioning of the matrix is
61: ! determined by PETSc at runtime.
63: PetscCallA(MatCreate(PETSC_COMM_WORLD, C, ierr))
64: PetscCallA(MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, m*n, m*n, ierr))
65: PetscCallA(MatSetFromOptions(C, ierr))
66: PetscCallA(MatSetUp(C, ierr))
68: ! Currently, all PETSc parallel matrix formats are partitioned by
69: ! contiguous chunks of rows across the processors. Determine which
70: ! rows of the matrix are locally owned.
72: PetscCallA(MatGetOwnershipRange(C, Istart, Iend, ierr))
74: ! Set matrix entries matrix in parallel.
75: ! - Each processor needs to insert only elements that it owns
76: ! locally (but any non-local elements will be sent to the
77: ! appropriate processor during matrix assembly).
78: !- Always specify global row and columns of matrix entries.
80: intitializeC: do Ii = Istart, Iend - 1
81: v = -1.0; i = Ii/n; j = Ii - i*n
82: if (i > 0) then
83: JJ = Ii - n
84: PetscCallA(MatSetValues(C, 1_PETSC_INT_KIND, [Ii], 1_PETSC_INT_KIND, [JJ], [v], ADD_VALUES, ierr))
85: end if
87: if (i < m - 1) then
88: JJ = Ii + n
89: PetscCallA(MatSetValues(C, 1_PETSC_INT_KIND, [Ii], 1_PETSC_INT_KIND, [JJ], [v], ADD_VALUES, ierr))
90: end if
92: if (j > 0) then
93: JJ = Ii - 1
94: PetscCallA(MatSetValues(C, 1_PETSC_INT_KIND, [Ii], 1_PETSC_INT_KIND, [JJ], [v], ADD_VALUES, ierr))
95: end if
97: if (j < n - 1) then
98: JJ = Ii + 1
99: PetscCallA(MatSetValues(C, 1_PETSC_INT_KIND, [Ii], 1_PETSC_INT_KIND, [JJ], [v], ADD_VALUES, ierr))
100: end if
102: v = 4.0
103: PetscCallA(MatSetValues(C, 1_PETSC_INT_KIND, [Ii], 1_PETSC_INT_KIND, [Ii], [v], ADD_VALUES, ierr))
104: end do intitializeC
106: ! Make the matrix nonsymmetric if desired
107: if (mat_nonsymmetric) then
108: do Ii = Istart, Iend - 1
109: v = -1.5; i = Ii/n
110: if (i > 1) then
111: JJ = Ii - n - 1
112: PetscCallA(MatSetValues(C, 1_PETSC_INT_KIND, [Ii], 1_PETSC_INT_KIND, [JJ], [v], ADD_VALUES, ierr))
113: end if
114: end do
115: else
116: PetscCallA(MatSetOption(C, MAT_SYMMETRIC, PETSC_TRUE, ierr))
117: PetscCallA(MatSetOption(C, MAT_SYMMETRY_ETERNAL, PETSC_TRUE, ierr))
118: end if
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: PetscCallA(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY, ierr))
126: PetscCallA(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY, ierr))
128: ! Create parallel vectors.
129: ! - When using VecSetSizes(), we specify only the vector's global
130: ! dimension; the parallel partitioning is determined at runtime.
131: ! - Note: We form 1 vector from scratch and then duplicate as needed.
133: PetscCallA(VecCreate(PETSC_COMM_WORLD, u, ierr))
134: PetscCallA(VecSetSizes(u, PETSC_DECIDE, m*n, ierr))
135: PetscCallA(VecSetFromOptions(u, ierr))
136: PetscCallA(VecDuplicate(u, b, ierr))
137: PetscCallA(VecDuplicate(b, x, ierr))
139: ! Currently, all parallel PETSc vectors are partitioned by
140: ! contiguous chunks across the processors. Determine which
141: ! range of entries are locally owned.
143: PetscCallA(VecGetOwnershipRange(x, low, high, ierr))
145: !Set elements within the exact solution vector in parallel.
146: ! - Each processor needs to insert only elements that it owns
147: ! locally (but any non-local entries will be sent to the
148: ! appropriate processor during vector assembly).
149: ! - Always specify global locations of vector entries.
151: PetscCallA(VecGetLocalSize(x, ldim, ierr))
152: do i = 0, ldim - 1
153: iglobal = i + low
154: v = real(i + 100*rank)
155: PetscCallA(VecSetValues(u, 1_PETSC_INT_KIND, [iglobal], [v], INSERT_VALUES, ierr))
156: end do
158: ! Assemble vector, using the 2-step process:
159: ! VecAssemblyBegin(), VecAssemblyEnd()
160: ! Computations can be done while messages are in transition,
161: ! by placing code between these two statements.
162: PetscCallA(VecAssemblyBegin(u, ierr))
163: PetscCallA(VecAssemblyEnd(u, ierr))
165: ! Compute right-hand-side vector
167: PetscCallA(MatMult(C, u, b, ierr))
169: ! Create linear solver context
171: PetscCallA(KSPCreate(PETSC_COMM_WORLD, ksp, ierr))
173: ! Set operators. Here the matrix that defines the linear system
174: ! also serves as the matrix used to construct the preconditioner.
176: PetscCallA(KSPSetOperators(ksp, C, C, ierr))
178: ! Set runtime options (e.g., -ksp_type <type> -pc_type <type>)
180: PetscCallA(KSPSetFromOptions(ksp, ierr))
182: ! Solve linear system. Here we explicitly call KSPSetUp() for more
183: ! detailed performance monitoring of certain preconditioners, such
184: ! as ICC and ILU. This call is optional, as KSPSetUp() will
185: ! automatically be called within KSPSolve() if it hasn't been
186: ! called already.
188: PetscCallA(KSPSetUp(ksp, ierr))
190: ! Do not do this in application code, use -ksp_gmres_modifiedgramschmidt or -ksp_gmres_modifiedgramschmidt
191: if (orthog == 1) then
192: PetscCallA(KSPGMRESSetOrthogonalization(ksp, KSPGMRESModifiedGramSchmidtOrthogonalization, ierr))
193: else if (orthog == 2) then
194: PetscCallA(KSPGMRESSetOrthogonalization(ksp, KSPGMRESClassicalGramSchmidtOrthogonalization, ierr))
195: end if
197: PetscCallA(KSPSolve(ksp, b, x, ierr))
199: ! Check the residual
200: PetscCallA(VecAXPY(x, myNone, u, ierr))
201: PetscCallA(VecNorm(x, NORM_2, norm, ierr))
202: PetscCallA(VecNorm(b, NORM_2, bnorm, ierr))
204: PetscCallA(KSPGetIterationNumber(ksp, its, ierr))
205: if (.not. testscaledMat .or. norm/bnorm > PETSC_SMALL) then
206: write (outputString, '(a,f11.9,a,i2.2,a)') 'Relative norm of residual ', norm/bnorm, ', Iterations ', its, '\n'
207: PetscCallA(PetscPrintf(PETSC_COMM_WORLD, outputString, ierr))
208: end if
210: ! -------------- Stage 1: Solve Second System ----------------------
212: ! Solve another linear system with the same method. We reuse the KSP
213: ! context, matrix and vector data structures, and hence save the
214: ! overhead of creating new ones.
216: ! Indicate to PETSc profiling that we're concluding the first
217: ! stage with PetscLogStagePop(), and beginning the second stage with
218: ! PetscLogStagePush().
220: PetscCallA(PetscLogStagePop(ierr))
221: PetscCallA(PetscLogStagePush(stages(1), ierr))
223: ! Initialize all matrix entries to zero. MatZeroEntries() retains the
224: ! nonzero structure of the matrix for sparse formats.
226: PetscCallA(MatZeroEntries(C, ierr))
228: ! Assemble matrix again. Note that we retain the same matrix data
229: ! structure and the same nonzero pattern; we just change the values
230: ! of the matrix entries.
232: do i = 0, m - 1
233: do j = 2*rank, 2*rank + 1
234: v = -1.0; Ii = j + n*i
235: if (i > 0) then
236: JJ = Ii - n
237: PetscCallA(MatSetValues(C, 1_PETSC_INT_KIND, [Ii], 1_PETSC_INT_KIND, [JJ], [v], ADD_VALUES, ierr))
238: end if
240: if (i < m - 1) then
241: JJ = Ii + n
242: PetscCallA(MatSetValues(C, 1_PETSC_INT_KIND, [Ii], 1_PETSC_INT_KIND, [JJ], [v], ADD_VALUES, ierr))
243: end if
245: if (j > 0) then
246: JJ = Ii - 1
247: PetscCallA(MatSetValues(C, 1_PETSC_INT_KIND, [Ii], 1_PETSC_INT_KIND, [JJ], [v], ADD_VALUES, ierr))
248: end if
250: if (j < n - 1) then
251: JJ = Ii + 1
252: PetscCallA(MatSetValues(C, 1_PETSC_INT_KIND, [Ii], 1_PETSC_INT_KIND, [JJ], [v], ADD_VALUES, ierr))
253: end if
255: v = 6.0
256: PetscCallA(MatSetValues(C, 1_PETSC_INT_KIND, [Ii], 1_PETSC_INT_KIND, [Ii], [v], ADD_VALUES, ierr))
257: end do
258: end do
260: ! Make the matrix nonsymmetric if desired
262: if (mat_nonsymmetric) then
263: do Ii = Istart, Iend - 1
264: v = -1.5; i = Ii/n
265: if (i > 1) then
266: JJ = Ii - n - 1
267: PetscCallA(MatSetValues(C, 1_PETSC_INT_KIND, [Ii], 1_PETSC_INT_KIND, [JJ], [v], ADD_VALUES, ierr))
268: end if
269: end do
270: end if
272: ! Assemble matrix, using the 2-step process:
273: ! MatAssemblyBegin(), MatAssemblyEnd()
274: ! Computations can be done while messages are in transition
275: ! by placing code between these two statements.
277: PetscCallA(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY, ierr))
278: PetscCallA(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY, ierr))
280: if (testscaledMat) then
281: ! Scale a(0,0) and a(M-1,M-1)
283: if (rank /= 0) then
284: v = 6.0*0.00001; Ii = 0; JJ = 0
285: PetscCallA(MatSetValues(C, 1_PETSC_INT_KIND, [Ii], 1_PETSC_INT_KIND, [JJ], [v], INSERT_VALUES, ierr))
286: elseif (rank == size - 1) then
287: v = 6.0*0.00001; Ii = m*n - 1; JJ = m*n - 1
288: PetscCallA(MatSetValues(C, 1_PETSC_INT_KIND, [Ii], 1_PETSC_INT_KIND, [JJ], [v], INSERT_VALUES, ierr))
290: end if
292: PetscCallA(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY, ierr))
293: PetscCallA(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY, ierr))
295: ! Compute a new right-hand-side vector
297: PetscCallA(VecDestroy(u, ierr))
298: PetscCallA(VecCreate(PETSC_COMM_WORLD, u, ierr))
299: PetscCallA(VecSetSizes(u, PETSC_DECIDE, m*n, ierr))
300: PetscCallA(VecSetFromOptions(u, ierr))
302: PetscCallA(PetscRandomCreate(PETSC_COMM_WORLD, rctx, ierr))
303: PetscCallA(PetscRandomSetFromOptions(rctx, ierr))
304: PetscCallA(VecSetRandom(u, rctx, ierr))
305: PetscCallA(PetscRandomDestroy(rctx, ierr))
306: PetscCallA(VecAssemblyBegin(u, ierr))
307: PetscCallA(VecAssemblyEnd(u, ierr))
309: end if
311: PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-test_newMat', testnewC, flg, ierr))
313: if (testnewC) then
314: ! User may use a new matrix C with same nonzero pattern, e.g.
315: ! ex5 -ksp_monitor -mat_type sbaij -pc_type cholesky -pc_factor_mat_solver_type mumps -test_newMat
317: PetscCallA(MatDuplicate(C, MAT_COPY_VALUES, Ctmp, ierr))
318: PetscCallA(MatDestroy(C, ierr))
319: PetscCallA(MatDuplicate(Ctmp, MAT_COPY_VALUES, C, ierr))
320: PetscCallA(MatDestroy(Ctmp, ierr))
321: end if
323: PetscCallA(MatMult(C, u, b, ierr))
325: ! Set operators. Here the matrix that defines the linear system
326: ! also serves as the matrix used to construct the preconditioner.
328: PetscCallA(KSPSetOperators(ksp, C, C, ierr))
330: ! Solve linear system
331: PetscCallA(KSPSetUp(ksp, ierr))
332: PetscCallA(KSPSolve(ksp, b, x, ierr))
333: ! Check the residual
335: PetscCallA(VecAXPY(x, myNone, u, ierr))
336: PetscCallA(VecNorm(x, NORM_2, norm, ierr))
337: PetscCallA(VecNorm(b, NORM_2, bnorm, ierr))
338: PetscCallA(KSPGetIterationNumber(ksp, its, ierr))
339: if (.not. testscaledMat .or. norm/bnorm > PETSC_SMALL) then
340: write (outputString, '(a,f11.9,a,i2.2,a)') 'Relative norm of residual ', norm/bnorm, ', Iterations ', its, '\n'
341: PetscCallA(PetscPrintf(PETSC_COMM_WORLD, outputString, ierr))
342: end if
344: ! Free work space. All PETSc objects should be destroyed when they
345: ! are no longer needed.
347: PetscCallA(KSPDestroy(ksp, ierr))
348: PetscCallA(VecDestroy(u, ierr))
349: PetscCallA(VecDestroy(x, ierr))
350: PetscCallA(VecDestroy(b, ierr))
351: PetscCallA(MatDestroy(C, ierr))
353: ! Indicate to PETSc profiling that we're concluding the second stage
355: PetscCallA(PetscLogStagePop(ierr))
356: PetscCallA(PetscFinalize(ierr))
357: end program main
359: !/*TEST
360: !
361: ! test:
362: ! args: -pc_type jacobi -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always
363: !
364: ! test:
365: ! suffix: 2
366: ! nsize: 2
367: ! args: -pc_type jacobi -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always -ksp_rtol .000001
368: !
369: ! test:
370: ! suffix: 5
371: ! nsize: 2
372: ! args: -ksp_gmres_cgs_refinement_type refine_always -ksp_monitor draw::draw_lg -ksp_monitor_true_residual draw::draw_lg
373: ! output_file: output/ex5f_5.out
374: !
375: ! test:
376: ! suffix: asm
377: ! nsize: 4
378: ! args: -pc_type asm
379: ! output_file: output/ex5f_asm.out
380: !
381: ! test:
382: ! suffix: asm_baij
383: ! nsize: 4
384: ! args: -pc_type asm -mat_type baij
385: ! output_file: output/ex5f_asm.out
386: !
387: ! test:
388: ! suffix: redundant_0
389: ! args: -m 1000 -pc_type redundant -pc_redundant_number 1 -redundant_ksp_type gmres -redundant_pc_type jacobi
390: !
391: ! test:
392: ! suffix: redundant_1
393: ! nsize: 5
394: ! args: -pc_type redundant -pc_redundant_number 1 -redundant_ksp_type gmres -redundant_pc_type jacobi
395: !
396: ! test:
397: ! suffix: redundant_2
398: ! nsize: 5
399: ! args: -pc_type redundant -pc_redundant_number 3 -redundant_ksp_type gmres -redundant_pc_type jacobi
400: !
401: ! test:
402: ! suffix: redundant_3
403: ! nsize: 5
404: ! args: -pc_type redundant -pc_redundant_number 5 -redundant_ksp_type gmres -redundant_pc_type jacobi
405: !
406: ! test:
407: ! suffix: redundant_4
408: ! nsize: 5
409: ! args: -pc_type redundant -pc_redundant_number 3 -redundant_ksp_type gmres -redundant_pc_type jacobi -psubcomm_type interlaced
410: !
411: ! test:
412: ! suffix: superlu_dist
413: ! nsize: 15
414: ! requires: superlu_dist
415: ! args: -pc_type lu -pc_factor_mat_solver_type superlu_dist -mat_superlu_dist_equil false -m 150 -mat_superlu_dist_r 3 -mat_superlu_dist_c 5 -test_scaledMat
416: ! output_file: output/empty.out
417: !
418: ! test:
419: ! suffix: superlu_dist_2
420: ! nsize: 15
421: ! requires: superlu_dist
422: ! args: -pc_type lu -pc_factor_mat_solver_type superlu_dist -mat_superlu_dist_equil false -m 150 -mat_superlu_dist_r 3 -mat_superlu_dist_c 5 -test_scaledMat -mat_superlu_dist_fact SamePattern_SameRowPerm
423: ! output_file: output/empty.out
424: !
425: ! test:
426: ! suffix: superlu_dist_3
427: ! nsize: 15
428: ! requires: superlu_dist
429: ! args: -pc_type lu -pc_factor_mat_solver_type superlu_dist -mat_superlu_dist_equil false -m 500 -mat_superlu_dist_r 3 -mat_superlu_dist_c 5 -test_scaledMat -mat_superlu_dist_fact DOFACT
430: ! output_file: output/empty.out
431: !
432: ! test:
433: ! suffix: superlu_dist_0
434: ! nsize: 1
435: ! requires: superlu_dist
436: ! args: -pc_type lu -pc_factor_mat_solver_type superlu_dist -test_scaledMat
437: ! output_file: output/empty.out
438: !
439: ! test:
440: ! suffix: orthog1
441: ! args: -orthog 1 -ksp_view
442: !
443: ! test:
444: ! suffix: orthog2
445: ! args: -orthog 2 -ksp_view
446: !
447: !TEST*/