Actual source code: ex5f.F90
petsc-3.12.0 2019-09-29
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
15:
16: implicit none
17: KSP :: myKsp ! 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
27: PetscMPIInt :: mySize,myRank
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=80) :: 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
43:
44: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-m',m,flg,ierr)
45: CHKERRA(ierr)
46: call MPI_Comm_rank(PETSC_COMM_WORLD,myRank,ierr)
47: CHKERRA(ierr)
48: call MPI_Comm_size(PETSC_COMM_WORLD,mySize,ierr)
49: CHKERRA(ierr)
50: n=2*mySize
52:
53: ! Set flag if we are doing a nonsymmetric problem; the default is symmetric.
54:
55: call PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,"-mat_nonsym",mat_nonsymmetric,flg,ierr)
56: CHKERRA(ierr)
57: call PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,"-test_scaledMat",testscaledMat,flg,ierr)
58: CHKERRA(ierr)
59:
60: ! Register two stages for separate profiling of the two linear solves.
61: ! Use the runtime option -log_view for a printout of performance
62: ! statistics at the program's conlusion.
64: call PetscLogStageRegister("Original Solve",stages(0),ierr)
65: CHKERRA(ierr)
66: call PetscLogStageRegister("Second Solve",stages(1),ierr)
67: CHKERRA(ierr)
68:
69: ! -------------- Stage 0: Solve Original System ----------------------
70: ! Indicate to PETSc profiling that we're beginning the first stage
71:
72: call PetscLogStagePush(stages(0),ierr)
73: CHKERRA(ierr)
74:
75: ! Create parallel matrix, specifying only its global dimensions.
76: ! When using MatCreate(), the matrix format can be specified at
77: ! runtime. Also, the parallel partitioning of the matrix is
78: ! determined by PETSc at runtime.
79:
80: call MatCreate(PETSC_COMM_WORLD,C,ierr)
81: CHKERRA(ierr)
82: call MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n,ierr)
83: CHKERRA(ierr)
84: call MatSetFromOptions(C,ierr)
85: CHKERRA(ierr)
86: call MatSetUp(C,ierr)
87: CHKERRA(ierr)
89: ! Currently, all PETSc parallel matrix formats are partitioned by
90: ! contiguous chunks of rows across the processors. Determine which
91: ! rows of the matrix are locally owned.
92:
93: call MatGetOwnershipRange(C,Istart,Iend,ierr)
94:
95: ! Set matrix entries matrix in parallel.
96: ! - Each processor needs to insert only elements that it owns
97: ! locally (but any non-local elements will be sent to the
98: ! appropriate processor during matrix assembly).
99: !- Always specify global row and columns of matrix entries.
101: intitializeC: do Ii=Istart,Iend-1
102: v =-1.0; i = Ii/n; j = Ii - i*n
103: if (i>0) then
104: JJ = Ii - n
105: call MatSetValues(C,one,Ii,one,JJ,v,ADD_VALUES,ierr)
106: CHKERRA(ierr)
107: endif
108:
109: if (i<m-1) then
110: JJ = Ii + n
111: call MatSetValues(C,one,Ii,one,JJ,v,ADD_VALUES,ierr)
112: CHKERRA(ierr)
113: endif
114:
115: if (j>0) then
116: JJ = Ii - 1
117: call MatSetValues(C,one,Ii,one,JJ,v,ADD_VALUES,ierr)
118: CHKERRA(ierr)
119: endif
120:
121: if (j<n-1) then
122: JJ = Ii + 1
123: call MatSetValues(C,one,Ii,one,JJ,v,ADD_VALUES,ierr)
124: CHKERRA(ierr)
125: endif
126:
127: v=4.0
128: call MatSetValues(C,one,Ii,one,Ii,v,ADD_VALUES,ierr)
129: CHKERRA(ierr)
130:
131: enddo intitializeC
132:
133: ! Make the matrix nonsymmetric if desired
134: if (mat_nonsymmetric) then
135: do Ii=Istart,Iend-1
136: v=-1.5; i=Ii/n
137: if (i>1) then
138: JJ=Ii-n-1
139: call MatSetValues(C,one,Ii,one,JJ,v,ADD_VALUES,ierr)
140: CHKERRA(ierr)
141: endif
142: enddo
143: else
144: call MatSetOption(C,MAT_SYMMETRIC,PETSC_TRUE,ierr)
145: CHKERRA(ierr)
146: call MatSetOption(C,MAT_SYMMETRY_ETERNAL,PETSC_TRUE,ierr)
147: CHKERRA(ierr)
148: endif
149:
150: ! Assemble matrix, using the 2-step process:
151: ! MatAssemblyBegin(), MatAssemblyEnd()
152: ! Computations can be done while messages are in transition
153: ! by placing code between these two statements.
154:
155: call MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY,ierr)
156: CHKERRA(ierr)
157: call MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY,ierr)
158: CHKERRA(ierr)
160: ! Create parallel vectors.
161: ! - When using VecSetSizes(), we specify only the vector's global
162: ! dimension; the parallel partitioning is determined at runtime.
163: ! - Note: We form 1 vector from scratch and then duplicate as needed.
164:
165: call VecCreate(PETSC_COMM_WORLD,u,ierr)
166: call VecSetSizes(u,PETSC_DECIDE,m*n,ierr)
167: call VecSetFromOptions(u,ierr)
168: call VecDuplicate(u,b,ierr)
169: call VecDuplicate(b,x,ierr)
171: ! Currently, all parallel PETSc vectors are partitioned by
172: ! contiguous chunks across the processors. Determine which
173: ! range of entries are locally owned.
174:
175: call VecGetOwnershipRange(x,low,high,ierr)
176: CHKERRA(ierr)
178: !Set elements within the exact solution vector in parallel.
179: ! - Each processor needs to insert only elements that it owns
180: ! locally (but any non-local entries will be sent to the
181: ! appropriate processor during vector assembly).
182: ! - Always specify global locations of vector entries.
184: call VecGetLocalSize(x,ldim,ierr)
185: CHKERRA(ierr)
186: do i=0,ldim-1
187: iglobal = i + low
188: v = real(i + 100*myRank)
189: call VecSetValues(u,one,iglobal,v,INSERT_VALUES,ierr)
190: CHKERRA(ierr)
191: enddo
193: ! Assemble vector, using the 2-step process:
194: ! VecAssemblyBegin(), VecAssemblyEnd()
195: ! Computations can be done while messages are in transition,
196: ! by placing code between these two statements.
197:
198: call VecAssemblyBegin(u,ierr)
199: CHKERRA(ierr)
200: call VecAssemblyEnd(u,ierr)
201: CHKERRA(ierr)
202:
203: ! Compute right-hand-side vector
204:
205: call MatMult(C,u,b,ierr)
207: CHKERRA(ierr)
208:
209: ! Create linear solver context
211: call KSPCreate(PETSC_COMM_WORLD,myKsp,ierr)
212: CHKERRA(ierr)
213: ! Set operators. Here the matrix that defines the linear system
214: ! also serves as the preconditioning matrix.
216: call KSPSetOperators(myKsp,C,C,ierr)
217: CHKERRA(ierr)
218: ! Set runtime options (e.g., -ksp_type <type> -pc_type <type>)
220: call KSPSetFromOptions(myKsp,ierr)
221: CHKERRA(ierr)
222: ! Solve linear system. Here we explicitly call KSPSetUp() for more
223: ! detailed performance monitoring of certain preconditioners, such
224: ! as ICC and ILU. This call is optional, as KSPSetUp() will
225: ! automatically be called within KSPSolve() if it hasn't been
226: ! called already.
228: call KSPSetUp(myKsp,ierr)
229: CHKERRA(ierr)
231: call KSPSolve(myKsp,b,x,ierr)
232: CHKERRA(ierr)
234: ! Check the error
235:
236: call VecAXPY(x,myNone,u,ierr)
237: call VecNorm(x,NORM_2,norm,ierr)
239: call KSPGetIterationNumber(myKsp,its,ierr)
240: if (.not. testscaledMat .or. norm > 1.e-7) then
241: write(outputString,'(a,f11.9,a,i2.2,a)') 'Norm of error ',norm,', Iterations ',its,'\n'
242: call PetscPrintf(PETSC_COMM_WORLD,outputString,ierr)
243: endif
244:
245: ! -------------- Stage 1: Solve Second System ----------------------
246:
247: ! Solve another linear system with the same method. We reuse the KSP
248: ! context, matrix and vector data structures, and hence save the
249: ! overhead of creating new ones.
251: ! Indicate to PETSc profiling that we're concluding the first
252: ! stage with PetscLogStagePop(), and beginning the second stage with
253: ! PetscLogStagePush().
254:
255: call PetscLogStagePop(ierr)
256: CHKERRA(ierr)
257: call PetscLogStagePush(stages(1),ierr)
258: CHKERRA(ierr)
259:
260: ! Initialize all matrix entries to zero. MatZeroEntries() retains the
261: ! nonzero structure of the matrix for sparse formats.
262:
263: call MatZeroEntries(C,ierr)
264: CHKERRA(ierr)
265:
266: ! Assemble matrix again. Note that we retain the same matrix data
267: ! structure and the same nonzero pattern; we just change the values
268: ! of the matrix entries.
269:
270: do i=0,m-1
271: do j=2*myRank,2*myRank+1
272: v =-1.0; Ii=j + n*i
273: if (i>0) then
274: JJ = Ii - n
275: call MatSetValues(C,one,Ii,one,JJ,v,ADD_VALUES,ierr)
276: CHKERRA(ierr)
277: endif
278:
279: if (i<m-1) then
280: JJ = Ii + n
281: call MatSetValues(C,one,Ii,one,JJ,v,ADD_VALUES,ierr)
282: CHKERRA(ierr)
283: endif
284:
285: if (j>0) then
286: JJ = Ii - 1
287: call MatSetValues(C,one,Ii,one,JJ,v,ADD_VALUES,ierr)
288: CHKERRA(ierr)
289: endif
290:
291: if (j<n-1) then
292: JJ = Ii + 1
293: call MatSetValues(C,one,Ii,one,JJ,v,ADD_VALUES,ierr)
294: CHKERRA(ierr)
295: endif
296:
297: v=6.0
298: call MatSetValues(C,one,Ii,one,Ii,v,ADD_VALUES,ierr)
299: CHKERRA(ierr)
300:
301: enddo
302: enddo
303:
304: ! Make the matrix nonsymmetric if desired
305:
306: if (mat_nonsymmetric) then
307: do Ii=Istart,Iend-1
308: v=-1.5; i=Ii/n
309: if (i>1) then
310: JJ=Ii-n-1
311: call MatSetValues(C,one,Ii,one,JJ,v,ADD_VALUES,ierr)
312: CHKERRA(ierr)
313: endif
314: enddo
315: endif
316:
317: ! Assemble matrix, using the 2-step process:
318: ! MatAssemblyBegin(), MatAssemblyEnd()
319: ! Computations can be done while messages are in transition
320: ! by placing code between these two statements.
321:
322: call MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY,ierr)
323: CHKERRA(ierr)
324: call MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY,ierr)
325: CHKERRA(ierr)
326:
327: if (testscaledMat) then
328: ! Scale a(0,0) and a(M-1,M-1)
330: if (myRank /= 0) then
331: v = 6.0*0.00001; Ii = 0; JJ = 0
332: call MatSetValues(C,one,Ii,one,JJ,v,INSERT_VALUES,ierr)
333: CHKERRA(ierr)
334: elseif (myRank == mySize -1) then
335: v = 6.0*0.00001; Ii = m*n-1; JJ = m*n-1
336: call MatSetValues(C,one,Ii,one,JJ,v,INSERT_VALUES,ierr)
337:
338: endif
339:
340: call MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY,ierr)
341: call MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY,ierr)
342:
343: ! Compute a new right-hand-side vector
344:
345: call VecDestroy(u,ierr)
346: call VecCreate(PETSC_COMM_WORLD,u,ierr)
347: call VecSetSizes(u,PETSC_DECIDE,m*n,ierr)
348: call VecSetFromOptions(u,ierr)
350: call PetscRandomCreate(PETSC_COMM_WORLD,rctx,ierr)
351: call PetscRandomSetFromOptions(rctx,ierr)
352: call VecSetRandom(u,rctx,ierr)
353: call PetscRandomDestroy(rctx,ierr)
354: call VecAssemblyBegin(u,ierr)
355: call VecAssemblyEnd(u,ierr)
356:
357: endif
358:
359: call PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,"-test_newMat",testnewC,flg,ierr)
360: CHKERRA(ierr)
361:
362: if (testnewC) then
363: ! User may use a new matrix C with same nonzero pattern, e.g.
364: ! ex5 -ksp_monitor -mat_type sbaij -pc_type cholesky -pc_factor_mat_solver_type mumps -test_newMat
366: call MatDuplicate(C,MAT_COPY_VALUES,Ctmp,ierr)
367: call MatDestroy(C,ierr)
368: call MatDuplicate(Ctmp,MAT_COPY_VALUES,C,ierr)
369: call MatDestroy(Ctmp,ierr)
370: endif
371:
372: call MatMult(C,u,b,ierr);CHKERRA(ierr)
374: ! Set operators. Here the matrix that defines the linear system
375: ! also serves as the preconditioning matrix.
376:
377: call KSPSetOperators(myKsp,C,C,ierr);CHKERRA(ierr)
378:
379: ! Solve linear system
380:
381: call KSPSetUp(myKsp,ierr); CHKERRA(ierr)
382: call KSPSolve(myKsp,b,x,ierr); CHKERRA(ierr)
383:
384: ! Check the error
385:
386: call VecAXPY(x,myNone,u,ierr); CHKERRA(ierr)
387: call VecNorm(x,NORM_2,norm,ierr); CHKERRA(ierr)
388: call KSPGetIterationNumber(myKsp,its,ierr); CHKERRA(ierr)
389: if (.not. testscaledMat .or. norm > 1.e-7) then
390: write(outputString,'(a,f11.9,a,i2.2,a)') 'Norm of error ',norm,', Iterations ',its,'\n'
391: call PetscPrintf(PETSC_COMM_WORLD,outputString,ierr)
392: endif
393:
394: ! Free work space. All PETSc objects should be destroyed when they
395: ! are no longer needed.
396:
397: call KSPDestroy(myKsp,ierr); CHKERRA(ierr)
398: call VecDestroy(u,ierr); CHKERRA(ierr)
399: call VecDestroy(x,ierr); CHKERRA(ierr)
400: call VecDestroy(b,ierr); CHKERRA(ierr)
401: call MatDestroy(C,ierr); CHKERRA(ierr)
402:
403: ! Indicate to PETSc profiling that we're concluding the second stage
404:
405: call PetscLogStagePop(ierr)
406: CHKERRA(ierr)
407:
408: call PetscFinalize(ierr)
409:
410: !/*TEST
411: !
412: ! test:
413: ! args: -pc_type jacobi -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always
414: !
415: ! test:
416: ! suffix: 2
417: ! nsize: 2
418: ! args: -pc_type jacobi -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always -ksp_rtol .000001
419: !
420: ! test:
421: ! suffix: 5
422: ! nsize: 2
423: ! args: -ksp_gmres_cgs_refinement_type refine_always -ksp_monitor_lg_residualnorm -ksp_monitor_lg_true_residualnorm
424: ! output_file: output/ex5f_5.out
425: !
426: ! test:
427: ! suffix: asm
428: ! nsize: 4
429: ! args: -pc_type asm
430: ! output_file: output/ex5f_asm.out
431: !
432: ! test:
433: ! suffix: asm_baij
434: ! nsize: 4
435: ! args: -pc_type asm -mat_type baij
436: ! output_file: output/ex5f_asm.out
437: !
438: ! test:
439: ! suffix: redundant_0
440: ! args: -m 1000 -pc_type redundant -pc_redundant_number 1 -redundant_ksp_type gmres -redundant_pc_type jacobi
441: !
442: ! test:
443: ! suffix: redundant_1
444: ! nsize: 5
445: ! args: -pc_type redundant -pc_redundant_number 1 -redundant_ksp_type gmres -redundant_pc_type jacobi
446: !
447: ! test:
448: ! suffix: redundant_2
449: ! nsize: 5
450: ! args: -pc_type redundant -pc_redundant_number 3 -redundant_ksp_type gmres -redundant_pc_type jacobi
451: !
452: ! test:
453: ! suffix: redundant_3
454: ! nsize: 5
455: ! args: -pc_type redundant -pc_redundant_number 5 -redundant_ksp_type gmres -redundant_pc_type jacobi
456: !
457: ! test:
458: ! suffix: redundant_4
459: ! nsize: 5
460: ! args: -pc_type redundant -pc_redundant_number 3 -redundant_ksp_type gmres -redundant_pc_type jacobi -psubcomm_type interlaced
461: !
462: ! test:
463: ! suffix: superlu_dist
464: ! nsize: 15
465: ! requires: superlu_dist
466: ! 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
467: !
468: ! test:
469: ! suffix: superlu_dist_2
470: ! nsize: 15
471: ! requires: superlu_dist
472: ! 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
473: ! output_file: output/ex5f_superlu_dist.out
474: !
475: ! test:
476: ! suffix: superlu_dist_3
477: ! nsize: 15
478: ! requires: superlu_dist
479: ! 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
480: ! output_file: output/ex5f_superlu_dist.out
481: !
482: ! test:
483: ! suffix: superlu_dist_0
484: ! nsize: 1
485: ! requires: superlu_dist
486: ! args: -pc_type lu -pc_factor_mat_solver_type superlu_dist -test_scaledMat
487: ! output_file: output/ex5f_superlu_dist.out
488: !
489: !TEST*/
491: end program main