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