Actual source code: ex5f.F90
petsc-3.11.4 2019-09-28
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, parameter :: m = 3
27: PetscMPIInt :: mySize,myRank
28: PetscBool, parameter :: &
29: testnewC = PETSC_FALSE, &
30: testscaledMat = PETSC_FALSE, &
31: mat_nonsymmetric = PETSC_FALSE
32: PetscBool :: flg
33: PetscRandom :: rctx
34: PetscLogStage,dimension(2) :: stages
35: character(len=80) :: outputString
36:
37: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
38: if (ierr /= 0) then
39: write(6,*)'Unable to initialize PETSc'
40: stop
41: endif
42:
43: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-m',m,flg,ierr)
44: CHKERRA(ierr)
45: call MPI_Comm_rank(PETSC_COMM_WORLD,myRank,ierr)
46: CHKERRA(ierr)
47: call MPI_Comm_size(PETSC_COMM_WORLD,mySize,ierr)
48: CHKERRA(ierr)
49: n=2*mySize
51:
52: ! Set flag if we are doing a nonsymmetric problem; the default is symmetric.
53:
54: call PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,"-mat_nonsym",mat_nonsymmetric,flg,ierr)
55: CHKERRA(ierr)
56: call PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,"-test_scaledMat",testscaledMat,flg,ierr)
57: CHKERRA(ierr)
58:
59: ! Register two stages for separate profiling of the two linear solves.
60: ! Use the runtime option -log_view for a printout of performance
61: ! statistics at the program's conlusion.
63: call PetscLogStageRegister("Original Solve",stages(1),ierr)
64: CHKERRA(ierr)
65: call PetscLogStageRegister("Second Solve",stages(2),ierr)
66: CHKERRA(ierr)
67:
68: ! -------------- Stage 0: Solve Original System ----------------------
69: ! Indicate to PETSc profiling that we're beginning the first stage
70:
71: call PetscLogStagePush(stages(1),ierr)
72: CHKERRA(ierr)
73:
74: ! Create parallel matrix, specifying only its global dimensions.
75: ! When using MatCreate(), the matrix format can be specified at
76: ! runtime. Also, the parallel partitioning of the matrix is
77: ! determined by PETSc at runtime.
78:
79: call MatCreate(PETSC_COMM_WORLD,C,ierr)
80: CHKERRA(ierr)
81: call MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n,ierr)
82: CHKERRA(ierr)
83: call MatSetFromOptions(C,ierr)
84: CHKERRA(ierr)
85: call MatSetUp(C,ierr)
86: CHKERRA(ierr)
88: ! Currently, all PETSc parallel matrix formats are partitioned by
89: ! contiguous chunks of rows across the processors. Determine which
90: ! rows of the matrix are locally owned.
91:
92: call MatGetOwnershipRange(C,Istart,Iend,ierr)
93:
94: ! Set matrix entries matrix in parallel.
95: ! - Each processor needs to insert only elements that it owns
96: ! locally (but any non-local elements will be sent to the
97: ! appropriate processor during matrix assembly).
98: !- Always specify global row and columns of matrix entries.
100: intitializeC: do Ii=Istart,Iend-1
101: v =-1.0; i = Ii/n; j = Ii - i*n
102: if (i>0) then
103: JJ = Ii - n
104: call MatSetValues(C,1,Ii,1,JJ,v,ADD_VALUES,ierr)
105: CHKERRA(ierr)
106: endif
107:
108: if (i<m-1) then
109: JJ = Ii + n
110: call MatSetValues(C,1,Ii,1,JJ,v,ADD_VALUES,ierr)
111: CHKERRA(ierr)
112: endif
113:
114: if (j>0) then
115: JJ = Ii - 1
116: call MatSetValues(C,1,Ii,1,JJ,v,ADD_VALUES,ierr)
117: CHKERRA(ierr)
118: endif
119:
120: if (j<n-1) then
121: JJ = Ii + 1
122: call MatSetValues(C,1,Ii,1,JJ,v,ADD_VALUES,ierr)
123: CHKERRA(ierr)
124: endif
125:
126: v=4.0
127: call MatSetValues(C,1,Ii,1,Ii,v,ADD_VALUES,ierr)
128: CHKERRA(ierr)
129:
130: enddo intitializeC
131:
132: ! Make the matrix nonsymmetric if desired
133: if (mat_nonsymmetric) then
134: do Ii=Istart,Iend-1
135: v=-1.5; i=Ii/n
136: if (i>1) then
137: JJ=Ii-n-1
138: call MatSetValues(C,1,Ii,1,JJ,v,ADD_VALUES,ierr)
139: CHKERRA(ierr)
140: endif
141: enddo
142: else
143: call MatSetOption(C,MAT_SYMMETRIC,PETSC_TRUE,ierr)
144: CHKERRA(ierr)
145: call MatSetOption(C,MAT_SYMMETRY_ETERNAL,PETSC_TRUE,ierr)
146: CHKERRA(ierr)
147: endif
148:
149: ! Assemble matrix, using the 2-step process:
150: ! MatAssemblyBegin(), MatAssemblyEnd()
151: ! Computations can be done while messages are in transition
152: ! by placing code between these two statements.
153:
154: call MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY,ierr)
155: CHKERRA(ierr)
156: call MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY,ierr)
157: CHKERRA(ierr)
159: ! Create parallel vectors.
160: ! - When using VecSetSizes(), we specify only the vector's global
161: ! dimension; the parallel partitioning is determined at runtime.
162: ! - Note: We form 1 vector from scratch and then duplicate as needed.
163:
164: call VecCreate(PETSC_COMM_WORLD,u,ierr)
165: call VecSetSizes(u,PETSC_DECIDE,m*n,ierr)
166: call VecSetFromOptions(u,ierr)
167: call VecDuplicate(u,b,ierr)
168: call VecDuplicate(b,x,ierr)
170: ! Currently, all parallel PETSc vectors are partitioned by
171: ! contiguous chunks across the processors. Determine which
172: ! range of entries are locally owned.
173:
174: call VecGetOwnershipRange(x,low,high,ierr)
175: CHKERRA(ierr)
177: !Set elements within the exact solution vector in parallel.
178: ! - Each processor needs to insert only elements that it owns
179: ! locally (but any non-local entries will be sent to the
180: ! appropriate processor during vector assembly).
181: ! - Always specify global locations of vector entries.
183: call VecGetLocalSize(x,ldim,ierr)
184: CHKERRA(ierr)
185: do i=0,ldim-1
186: iglobal = i + low
187: v = real(i + 100*myRank)
188: call VecSetValues(u,1,iglobal,v,INSERT_VALUES,ierr)
189: CHKERRA(ierr)
190: enddo
192: ! Assemble vector, using the 2-step process:
193: ! VecAssemblyBegin(), VecAssemblyEnd()
194: ! Computations can be done while messages are in transition,
195: ! by placing code between these two statements.
196:
197: call VecAssemblyBegin(u,ierr)
198: CHKERRA(ierr)
199: call VecAssemblyEnd(u,ierr)
200: CHKERRA(ierr)
201:
202: ! Compute right-hand-side vector
203:
204: call MatMult(C,u,b,ierr)
206: CHKERRA(ierr)
207:
208: ! Create linear solver context
210: call KSPCreate(PETSC_COMM_WORLD,myKsp,ierr)
211: CHKERRA(ierr)
212: ! Set operators. Here the matrix that defines the linear system
213: ! also serves as the preconditioning matrix.
215: call KSPSetOperators(myKsp,C,C,ierr)
216: CHKERRA(ierr)
217: ! Set runtime options (e.g., -ksp_type <type> -pc_type <type>)
219: call KSPSetFromOptions(myKsp,ierr)
220: CHKERRA(ierr)
221: ! Solve linear system. Here we explicitly call KSPSetUp() for more
222: ! detailed performance monitoring of certain preconditioners, such
223: ! as ICC and ILU. This call is optional, as KSPSetUp() will
224: ! automatically be called within KSPSolve() if it hasn't been
225: ! called already.
227: call KSPSetUp(myKsp,ierr)
228: CHKERRA(ierr)
230: call KSPSolve(myKsp,b,x,ierr)
231: CHKERRA(ierr)
233: ! Check the error
234:
235: call VecAXPY(x,myNone,u,ierr)
236: call VecNorm(x,NORM_2,norm,ierr)
238: call KSPGetIterationNumber(myKsp,its,ierr)
239: if (.not. testscaledMat .or. norm > 1.e-7) then
240: write(outputString,*)'Norm of error',norm,'Iterations',its,'\n'
241: call PetscPrintf(PETSC_COMM_WORLD,outputString,ierr)
242: endif
243:
244: ! -------------- Stage 1: Solve Second System ----------------------
245:
246: ! Solve another linear system with the same method. We reuse the KSP
247: ! context, matrix and vector data structures, and hence save the
248: ! overhead of creating new ones.
250: ! Indicate to PETSc profiling that we're concluding the first
251: ! stage with PetscLogStagePop(), and beginning the second stage with
252: ! PetscLogStagePush().
253:
254: call PetscLogStagePop(ierr)
255: CHKERRA(ierr)
256: call PetscLogStagePush(stages(2),ierr)
257: CHKERRA(ierr)
258:
259: ! Initialize all matrix entries to zero. MatZeroEntries() retains the
260: ! nonzero structure of the matrix for sparse formats.
261:
262: call MatZeroEntries(C,ierr)
263: CHKERRA(ierr)
264:
265: ! Assemble matrix again. Note that we retain the same matrix data
266: ! structure and the same nonzero pattern; we just change the values
267: ! of the matrix entries.
268:
269: do i=0,m-1
270: do j=2*myRank,2*myRank+1
271: v =-1.0; Ii=j + n*i
272: if (i>0) then
273: JJ = Ii - n
274: call MatSetValues(C,1,Ii,1,JJ,v,ADD_VALUES,ierr)
275: CHKERRA(ierr)
276: endif
277:
278: if (i<m-1) then
279: JJ = Ii + n
280: call MatSetValues(C,1,Ii,1,JJ,v,ADD_VALUES,ierr)
281: CHKERRA(ierr)
282: endif
283:
284: if (j>0) then
285: JJ = Ii - 1
286: call MatSetValues(C,1,Ii,1,JJ,v,ADD_VALUES,ierr)
287: CHKERRA(ierr)
288: endif
289:
290: if (j<n-1) then
291: JJ = Ii + 1
292: call MatSetValues(C,1,Ii,1,JJ,v,ADD_VALUES,ierr)
293: CHKERRA(ierr)
294: endif
295:
296: v=6.0
297: call MatSetValues(C,1,Ii,1,JJ,v,ADD_VALUES,ierr)
298: CHKERRA(ierr)
299:
300: enddo
301: enddo
302:
303: ! Make the matrix nonsymmetric if desired
304:
305: if (mat_nonsymmetric) then
306: do Ii=Istart,Iend-1
307: v=-1.5; i=Ii/n
308: if (i>1) then
309: JJ=Ii-n-1
310: call MatSetValues(C,1,Ii,1,JJ,v,ADD_VALUES,ierr)
311: CHKERRA(ierr)
312: endif
313: enddo
314: else
315: call MatSetOption(C,MAT_SYMMETRIC,PETSC_TRUE,ierr)
316: CHKERRA(ierr)
317: call MatSetOption(C,MAT_SYMMETRY_ETERNAL,PETSC_TRUE,ierr)
318: CHKERRA(ierr)
319: endif
320:
321: ! Assemble matrix, using the 2-step process:
322: ! MatAssemblyBegin(), MatAssemblyEnd()
323: ! Computations can be done while messages are in transition
324: ! by placing code between these two statements.
325:
326: call MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY,ierr)
327: CHKERRA(ierr)
328: call MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY,ierr)
329: CHKERRA(ierr)
330:
331: if (testscaledMat) then
332: ! Scale a(0,0) and a(M-1,M-1)
334: if (myRank /= 0) then
335: v = 6.0*0.00001; Ii = 0; JJ = 0
336: call MatSetValues(C,1,Ii,1,JJ,v,INSERT_VALUES,ierr)
337: CHKERRA(ierr)
338: elseif (myRank == mySize -1) then
339: v = 6.0*0.00001; Ii = m*n-1; JJ = m*n-1
340: call MatSetValues(C,1,Ii,1,JJ,v,INSERT_VALUES,ierr)
341:
342: endif
343:
344: call MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY,ierr)
345: call MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY,ierr)
346:
347: ! Compute a new right-hand-side vector
348:
349: call VecDestroy(u,ierr)
350: call VecCreate(PETSC_COMM_WORLD,u,ierr)
351: call VecSetSizes(u,PETSC_DECIDE,m*n,ierr)
352: call VecSetFromOptions(u,ierr)
354: call PetscRandomCreate(PETSC_COMM_WORLD,rctx,ierr)
355: call PetscRandomSetFromOptions(rctx,ierr)
356: call VecSetRandom(u,rctx,ierr)
357: call PetscRandomDestroy(rctx,ierr)
358: call VecAssemblyBegin(u,ierr)
359: call VecAssemblyEnd(u,ierr)
360:
361: endif
362:
363: call PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,"-test_newMat",testnewC,flg,ierr)
364: CHKERRA(ierr)
365:
366: if (testnewC) then
367: ! User may use a new matrix C with same nonzero pattern, e.g.
368: ! ex5 -ksp_monitor -mat_type sbaij -pc_type cholesky -pc_factor_mat_solver_type mumps -test_newMat
370: call MatDuplicate(C,MAT_COPY_VALUES,Ctmp,ierr)
371: call MatDestroy(C,ierr)
372: call MatDuplicate(Ctmp,MAT_COPY_VALUES,C,ierr)
373: call MatDestroy(Ctmp,ierr)
374: endif
375:
376: call MatMult(C,u,b,ierr);CHKERRA(ierr)
378: ! Set operators. Here the matrix that defines the linear system
379: ! also serves as the preconditioning matrix.
380:
381: call KSPSetOperators(myKsp,C,C,ierr);CHKERRA(ierr)
382:
383: ! Solve linear system
384:
385: call KSPSetUp(myKsp,ierr); CHKERRA(ierr)
386: call KSPSolve(myKsp,b,x,ierr); CHKERRA(ierr)
387:
388: ! Check the error
389:
390: call VecAXPY(x,myNone,u,ierr); CHKERRA(ierr)
391: call VecNorm(x,NORM_2,norm,ierr); CHKERRA(ierr)
392: call KSPGetIterationNumber(myKsp,its,ierr); CHKERRA(ierr)
393: if (.not. testscaledMat .or. norm > 1.e-7) then
394: write(outputString,*)'Norm of error',norm,'Iterations',its,'\n'
395: call PetscPrintf(PETSC_COMM_WORLD,outputString,ierr)
396: endif
397:
398: ! Free work space. All PETSc objects should be destroyed when they
399: ! are no longer needed.
400:
401: call KSPDestroy(myKsp,ierr); CHKERRA(ierr)
402: call VecDestroy(u,ierr); CHKERRA(ierr)
403: call VecDestroy(x,ierr); CHKERRA(ierr)
404: call VecDestroy(b,ierr); CHKERRA(ierr)
405: call MatDestroy(C,ierr); CHKERRA(ierr)
406:
407: ! Indicate to PETSc profiling that we're concluding the second stage
408:
409: call PetscLogStagePop(ierr)
410: CHKERRA(ierr)
411:
412: call PetscFinalize(ierr)
413:
415: end program main