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)
8: program main
9: #include <petsc/finclude/petscksp.h>
10: use petscksp
12: implicit none
13: KSP :: ksp ! linear solver context
14: Mat :: C,Ctmp ! matrix
15: Vec :: x,u,b ! approx solution, RHS, exact solution
16: PetscReal :: norm ! norm of solution error
17: PetscScalar :: v
18: PetscScalar, parameter :: myNone = -1.0
19: PetscInt :: Ii,JJ,ldim,low,high,iglobal,Istart,Iend
20: PetscErrorCode :: ierr
21: PetscInt :: i,j,its,n
22: PetscInt :: m = 3, orthog = 0
23: PetscMPIInt :: size,rank
24: PetscBool :: &
25: testnewC = PETSC_FALSE, &
26: testscaledMat = PETSC_FALSE, &
27: mat_nonsymmetric = PETSC_FALSE
28: PetscBool :: flg
29: PetscRandom :: rctx
30: PetscLogStage,dimension(0:1) :: stages
31: character(len=PETSC_MAX_PATH_LEN) :: outputString
32: PetscInt,parameter :: one = 1
34: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
35: if (ierr /= 0) then
36: write(6,*)'Unable to initialize PETSc'
37: stop
38: endif
40: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-orthog',orthog,flg,ierr)
41: CHKERRA(ierr)
42: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-m',m,flg,ierr)
43: CHKERRA(ierr)
44: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
45: CHKERRA(ierr)
46: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
47: CHKERRA(ierr)
48: n=2*size
50: ! Set flag if we are doing a nonsymmetric problem; the default is symmetric.
52: call PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,"-mat_nonsym",mat_nonsymmetric,flg,ierr)
53: CHKERRA(ierr)
54: call PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,"-test_scaledMat",testscaledMat,flg,ierr)
55: CHKERRA(ierr)
57: ! Register two stages for separate profiling of the two linear solves.
58: ! Use the runtime option -log_view for a printout of performance
59: ! statistics at the program's conlusion.
61: call PetscLogStageRegister("Original Solve",stages(0),ierr)
62: CHKERRA(ierr)
63: call PetscLogStageRegister("Second Solve",stages(1),ierr)
64: CHKERRA(ierr)
66: ! -------------- Stage 0: Solve Original System ----------------------
67: ! Indicate to PETSc profiling that we're beginning the first stage
69: call PetscLogStagePush(stages(0),ierr)
70: CHKERRA(ierr)
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,C,ierr)
78: CHKERRA(ierr)
79: call MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n,ierr)
80: CHKERRA(ierr)
81: call MatSetFromOptions(C,ierr)
82: CHKERRA(ierr)
83: call MatSetUp(C,ierr)
84: CHKERRA(ierr)
86: ! Currently, all PETSc parallel matrix formats are partitioned by
87: ! contiguous chunks of rows across the processors. Determine which
88: ! rows of the matrix are locally owned.
90: call MatGetOwnershipRange(C,Istart,Iend,ierr)
92: ! Set matrix entries matrix in parallel.
93: ! - Each processor needs to insert only elements that it owns
94: ! locally (but any non-local elements will be sent to the
95: ! appropriate processor during matrix assembly).
96: !- Always specify global row and columns of matrix entries.
98: intitializeC: do Ii=Istart,Iend-1
99: v =-1.0; i = Ii/n; j = Ii - i*n
100: if (i>0) then
101: JJ = Ii - n
102: call MatSetValues(C,one,Ii,one,JJ,v,ADD_VALUES,ierr)
103: CHKERRA(ierr)
104: endif
106: if (i<m-1) then
107: JJ = Ii + n
108: call MatSetValues(C,one,Ii,one,JJ,v,ADD_VALUES,ierr)
109: CHKERRA(ierr)
110: endif
112: if (j>0) then
113: JJ = Ii - 1
114: call MatSetValues(C,one,Ii,one,JJ,v,ADD_VALUES,ierr)
115: CHKERRA(ierr)
116: endif
118: if (j<n-1) then
119: JJ = Ii + 1
120: call MatSetValues(C,one,Ii,one,JJ,v,ADD_VALUES,ierr)
121: CHKERRA(ierr)
122: endif
124: v=4.0
125: call MatSetValues(C,one,Ii,one,Ii,v,ADD_VALUES,ierr)
126: CHKERRA(ierr)
128: enddo intitializeC
130: ! Make the matrix nonsymmetric if desired
131: if (mat_nonsymmetric) then
132: do Ii=Istart,Iend-1
133: v=-1.5; i=Ii/n
134: if (i>1) then
135: JJ=Ii-n-1
136: call MatSetValues(C,one,Ii,one,JJ,v,ADD_VALUES,ierr)
137: CHKERRA(ierr)
138: endif
139: enddo
140: else
141: call MatSetOption(C,MAT_SYMMETRIC,PETSC_TRUE,ierr)
142: CHKERRA(ierr)
143: call MatSetOption(C,MAT_SYMMETRY_ETERNAL,PETSC_TRUE,ierr)
144: CHKERRA(ierr)
145: endif
147: ! Assemble matrix, using the 2-step process:
148: ! MatAssemblyBegin(), MatAssemblyEnd()
149: ! Computations can be done while messages are in transition
150: ! by placing code between these two statements.
152: call MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY,ierr)
153: CHKERRA(ierr)
154: call MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY,ierr)
155: CHKERRA(ierr)
157: ! Create parallel vectors.
158: ! - When using VecSetSizes(), we specify only the vector's global
159: ! dimension; the parallel partitioning is determined at runtime.
160: ! - Note: We form 1 vector from scratch and then duplicate as needed.
162: call VecCreate(PETSC_COMM_WORLD,u,ierr)
163: call VecSetSizes(u,PETSC_DECIDE,m*n,ierr)
164: call VecSetFromOptions(u,ierr)
165: call VecDuplicate(u,b,ierr)
166: call VecDuplicate(b,x,ierr)
168: ! Currently, all parallel PETSc vectors are partitioned by
169: ! contiguous chunks across the processors. Determine which
170: ! range of entries are locally owned.
172: call VecGetOwnershipRange(x,low,high,ierr)
173: CHKERRA(ierr)
175: !Set elements within the exact solution vector in parallel.
176: ! - Each processor needs to insert only elements that it owns
177: ! locally (but any non-local entries will be sent to the
178: ! appropriate processor during vector assembly).
179: ! - Always specify global locations of vector entries.
181: call VecGetLocalSize(x,ldim,ierr)
182: CHKERRA(ierr)
183: do i=0,ldim-1
184: iglobal = i + low
185: v = real(i + 100*rank)
186: call VecSetValues(u,one,iglobal,v,INSERT_VALUES,ierr)
187: CHKERRA(ierr)
188: enddo
190: ! Assemble vector, using the 2-step process:
191: ! VecAssemblyBegin(), VecAssemblyEnd()
192: ! Computations can be done while messages are in transition,
193: ! by placing code between these two statements.
195: call VecAssemblyBegin(u,ierr)
196: CHKERRA(ierr)
197: call VecAssemblyEnd(u,ierr)
198: CHKERRA(ierr)
200: ! Compute right-hand-side vector
202: call MatMult(C,u,b,ierr)
204: CHKERRA(ierr)
206: ! Create linear solver context
208: call KSPCreate(PETSC_COMM_WORLD,ksp,ierr)
209: CHKERRA(ierr)
210: ! Set operators. Here the matrix that defines the linear system
211: ! also serves as the preconditioning matrix.
213: call KSPSetOperators(ksp,C,C,ierr)
214: CHKERRA(ierr)
215: ! Set runtime options (e.g., -ksp_type <type> -pc_type <type>)
217: call KSPSetFromOptions(ksp,ierr)
218: CHKERRA(ierr)
219: ! Solve linear system. Here we explicitly call KSPSetUp() for more
220: ! detailed performance monitoring of certain preconditioners, such
221: ! as ICC and ILU. This call is optional, as KSPSetUp() will
222: ! automatically be called within KSPSolve() if it hasn't been
223: ! called already.
225: call KSPSetUp(ksp,ierr)
226: CHKERRA(ierr)
228: ! Do not do this in application code, use -ksp_gmres_modifiedgramschmidt or -ksp_gmres_modifiedgramschmidt
229: if (orthog .eq. 1) then
230: call KSPGMRESSetOrthogonalization(ksp,KSPGMRESModifiedGramSchmidtOrthogonalization,ierr)
231: else if (orthog .eq. 2) then
232: call KSPGMRESSetOrthogonalization(ksp,KSPGMRESClassicalGramSchmidtOrthogonalization,ierr)
233: endif
234: CHKERRA(ierr)
236: call KSPSolve(ksp,b,x,ierr)
237: CHKERRA(ierr)
239: ! Check the error
241: call VecAXPY(x,myNone,u,ierr)
242: call VecNorm(x,NORM_2,norm,ierr)
244: call KSPGetIterationNumber(ksp,its,ierr)
245: if (.not. testscaledMat .or. norm > 1.e-7) then
246: write(outputString,'(a,f11.9,a,i2.2,a)') 'Norm of error ',norm,', Iterations ',its,'\n'
247: call PetscPrintf(PETSC_COMM_WORLD,outputString,ierr)
248: endif
250: ! -------------- Stage 1: Solve Second System ----------------------
252: ! Solve another linear system with the same method. We reuse the KSP
253: ! context, matrix and vector data structures, and hence save the
254: ! overhead of creating new ones.
256: ! Indicate to PETSc profiling that we're concluding the first
257: ! stage with PetscLogStagePop(), and beginning the second stage with
258: ! PetscLogStagePush().
260: call PetscLogStagePop(ierr)
261: CHKERRA(ierr)
262: call PetscLogStagePush(stages(1),ierr)
263: CHKERRA(ierr)
265: ! Initialize all matrix entries to zero. MatZeroEntries() retains the
266: ! nonzero structure of the matrix for sparse formats.
268: call MatZeroEntries(C,ierr)
269: CHKERRA(ierr)
271: ! Assemble matrix again. Note that we retain the same matrix data
272: ! structure and the same nonzero pattern; we just change the values
273: ! of the matrix entries.
275: do i=0,m-1
276: do j=2*rank,2*rank+1
277: v =-1.0; Ii=j + n*i
278: if (i>0) then
279: JJ = Ii - n
280: call MatSetValues(C,one,Ii,one,JJ,v,ADD_VALUES,ierr)
281: CHKERRA(ierr)
282: endif
284: if (i<m-1) then
285: JJ = Ii + n
286: call MatSetValues(C,one,Ii,one,JJ,v,ADD_VALUES,ierr)
287: CHKERRA(ierr)
288: endif
290: if (j>0) then
291: JJ = Ii - 1
292: call MatSetValues(C,one,Ii,one,JJ,v,ADD_VALUES,ierr)
293: CHKERRA(ierr)
294: endif
296: if (j<n-1) then
297: JJ = Ii + 1
298: call MatSetValues(C,one,Ii,one,JJ,v,ADD_VALUES,ierr)
299: CHKERRA(ierr)
300: endif
302: v=6.0
303: call MatSetValues(C,one,Ii,one,Ii,v,ADD_VALUES,ierr)
304: CHKERRA(ierr)
306: enddo
307: enddo
309: ! Make the matrix nonsymmetric if desired
311: if (mat_nonsymmetric) then
312: do Ii=Istart,Iend-1
313: v=-1.5; i=Ii/n
314: if (i>1) then
315: JJ=Ii-n-1
316: call MatSetValues(C,one,Ii,one,JJ,v,ADD_VALUES,ierr)
317: CHKERRA(ierr)
318: endif
319: enddo
320: endif
322: ! Assemble matrix, using the 2-step process:
323: ! MatAssemblyBegin(), MatAssemblyEnd()
324: ! Computations can be done while messages are in transition
325: ! by placing code between these two statements.
327: call MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY,ierr)
328: CHKERRA(ierr)
329: call MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY,ierr)
330: CHKERRA(ierr)
332: if (testscaledMat) then
333: ! Scale a(0,0) and a(M-1,M-1)
335: if (rank /= 0) then
336: v = 6.0*0.00001; Ii = 0; JJ = 0
337: call MatSetValues(C,one,Ii,one,JJ,v,INSERT_VALUES,ierr)
338: CHKERRA(ierr)
339: elseif (rank == size -1) then
340: v = 6.0*0.00001; Ii = m*n-1; JJ = m*n-1
341: call MatSetValues(C,one,Ii,one,JJ,v,INSERT_VALUES,ierr)
343: endif
345: call MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY,ierr)
346: call MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY,ierr)
348: ! Compute a new right-hand-side vector
350: call VecDestroy(u,ierr)
351: call VecCreate(PETSC_COMM_WORLD,u,ierr)
352: call VecSetSizes(u,PETSC_DECIDE,m*n,ierr)
353: call VecSetFromOptions(u,ierr)
355: call PetscRandomCreate(PETSC_COMM_WORLD,rctx,ierr)
356: call PetscRandomSetFromOptions(rctx,ierr)
357: call VecSetRandom(u,rctx,ierr)
358: call PetscRandomDestroy(rctx,ierr)
359: call VecAssemblyBegin(u,ierr)
360: call VecAssemblyEnd(u,ierr)
362: endif
364: call PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,"-test_newMat",testnewC,flg,ierr)
365: CHKERRA(ierr)
367: if (testnewC) then
368: ! User may use a new matrix C with same nonzero pattern, e.g.
369: ! ex5 -ksp_monitor -mat_type sbaij -pc_type cholesky -pc_factor_mat_solver_type mumps -test_newMat
371: call MatDuplicate(C,MAT_COPY_VALUES,Ctmp,ierr)
372: call MatDestroy(C,ierr)
373: call MatDuplicate(Ctmp,MAT_COPY_VALUES,C,ierr)
374: call MatDestroy(Ctmp,ierr)
375: endif
377: call MatMult(C,u,b,ierr);CHKERRA(ierr)
379: ! Set operators. Here the matrix that defines the linear system
380: ! also serves as the preconditioning matrix.
382: call KSPSetOperators(ksp,C,C,ierr);CHKERRA(ierr)
384: ! Solve linear system
386: call KSPSetUp(ksp,ierr); CHKERRA(ierr)
387: call KSPSolve(ksp,b,x,ierr); CHKERRA(ierr)
389: ! Check the error
391: call VecAXPY(x,myNone,u,ierr); CHKERRA(ierr)
392: call VecNorm(x,NORM_2,norm,ierr); CHKERRA(ierr)
393: call KSPGetIterationNumber(ksp,its,ierr); CHKERRA(ierr)
394: if (.not. testscaledMat .or. norm > 1.e-7) then
395: write(outputString,'(a,f11.9,a,i2.2,a)') 'Norm of error ',norm,', Iterations ',its,'\n'
396: call PetscPrintf(PETSC_COMM_WORLD,outputString,ierr)
397: endif
399: ! Free work space. All PETSc objects should be destroyed when they
400: ! are no longer needed.
402: call KSPDestroy(ksp,ierr); CHKERRA(ierr)
403: call VecDestroy(u,ierr); CHKERRA(ierr)
404: call VecDestroy(x,ierr); CHKERRA(ierr)
405: call VecDestroy(b,ierr); CHKERRA(ierr)
406: call MatDestroy(C,ierr); CHKERRA(ierr)
408: ! Indicate to PETSc profiling that we're concluding the second stage
410: call PetscLogStagePop(ierr)
411: CHKERRA(ierr)
413: call PetscFinalize(ierr)
415: !/*TEST
416: !
417: ! test:
418: ! args: -pc_type jacobi -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always
419: !
420: ! test:
421: ! suffix: 2
422: ! nsize: 2
423: ! args: -pc_type jacobi -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always -ksp_rtol .000001
424: !
425: ! test:
426: ! suffix: 5
427: ! nsize: 2
428: ! args: -ksp_gmres_cgs_refinement_type refine_always -ksp_monitor draw::draw_lg -ksp_monitor_true_residual draw::draw_lg
429: ! output_file: output/ex5f_5.out
430: !
431: ! test:
432: ! suffix: asm
433: ! nsize: 4
434: ! args: -pc_type asm
435: ! output_file: output/ex5f_asm.out
436: !
437: ! test:
438: ! suffix: asm_baij
439: ! nsize: 4
440: ! args: -pc_type asm -mat_type baij
441: ! output_file: output/ex5f_asm.out
442: !
443: ! test:
444: ! suffix: redundant_0
445: ! args: -m 1000 -pc_type redundant -pc_redundant_number 1 -redundant_ksp_type gmres -redundant_pc_type jacobi
446: !
447: ! test:
448: ! suffix: redundant_1
449: ! nsize: 5
450: ! args: -pc_type redundant -pc_redundant_number 1 -redundant_ksp_type gmres -redundant_pc_type jacobi
451: !
452: ! test:
453: ! suffix: redundant_2
454: ! nsize: 5
455: ! args: -pc_type redundant -pc_redundant_number 3 -redundant_ksp_type gmres -redundant_pc_type jacobi
456: !
457: ! test:
458: ! suffix: redundant_3
459: ! nsize: 5
460: ! args: -pc_type redundant -pc_redundant_number 5 -redundant_ksp_type gmres -redundant_pc_type jacobi
461: !
462: ! test:
463: ! suffix: redundant_4
464: ! nsize: 5
465: ! args: -pc_type redundant -pc_redundant_number 3 -redundant_ksp_type gmres -redundant_pc_type jacobi -psubcomm_type interlaced
466: !
467: ! test:
468: ! suffix: superlu_dist
469: ! nsize: 15
470: ! requires: superlu_dist
471: ! 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
472: !
473: ! test:
474: ! suffix: superlu_dist_2
475: ! nsize: 15
476: ! requires: superlu_dist
477: ! 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
478: ! output_file: output/ex5f_superlu_dist.out
479: !
480: ! test:
481: ! suffix: superlu_dist_3
482: ! nsize: 15
483: ! requires: superlu_dist
484: ! 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
485: ! output_file: output/ex5f_superlu_dist.out
486: !
487: ! test:
488: ! suffix: superlu_dist_0
489: ! nsize: 1
490: ! requires: superlu_dist
491: ! args: -pc_type lu -pc_factor_mat_solver_type superlu_dist -test_scaledMat
492: ! output_file: output/ex5f_superlu_dist.out
493: !
494: ! test:
495: ! suffix: orthog1
496: ! args: -orthog 1 -ksp_view
497: !
498: ! test:
499: ! suffix: orthog2
500: ! args: -orthog 2 -ksp_view
501: !
502: !TEST*/
504: end program main