Actual source code: ex9.c
petsc-3.13.6 2020-09-29
2: static char help[] = "The solution of 2 different linear systems with different linear solvers.\n\
3: Also, this example illustrates the repeated\n\
4: solution of linear systems, while reusing matrix, vector, and solver data\n\
5: structures throughout the process. Note the various stages of event logging.\n\n";
7: /*T
8: Concepts: KSP^repeatedly solving linear systems;
9: Concepts: PetscLog^profiling multiple stages of code;
10: Concepts: PetscLog^user-defined event profiling;
11: Processors: n
12: T*/
16: /*
17: Include "petscksp.h" so that we can use KSP solvers. Note that this file
18: automatically includes:
19: petscsys.h - base PETSc routines petscvec.h - vectors
20: petscmat.h - matrices
21: petscis.h - index sets petscksp.h - Krylov subspace methods
22: petscviewer.h - viewers petscpc.h - preconditioners
23: */
24: #include <petscksp.h>
26: /*
27: Declare user-defined routines
28: */
29: extern PetscErrorCode CheckError(Vec,Vec,Vec,PetscInt,PetscReal,PetscLogEvent);
30: extern PetscErrorCode MyKSPMonitor(KSP,PetscInt,PetscReal,void*);
32: int main(int argc,char **args)
33: {
34: Vec x1,b1,x2,b2; /* solution and RHS vectors for systems #1 and #2 */
35: Vec u; /* exact solution vector */
36: Mat C1,C2; /* matrices for systems #1 and #2 */
37: KSP ksp1,ksp2; /* KSP contexts for systems #1 and #2 */
38: PetscInt ntimes = 3; /* number of times to solve the linear systems */
39: PetscLogEvent CHECK_ERROR; /* event number for error checking */
40: PetscInt ldim,low,high,iglobal,Istart,Iend,Istart2,Iend2;
41: PetscInt Ii,J,i,j,m = 3,n = 2,its,t;
43: PetscBool flg = PETSC_FALSE, unsym = PETSC_TRUE;
44: PetscScalar v;
45: PetscMPIInt rank,size;
46: #if defined(PETSC_USE_LOG)
47: PetscLogStage stages[3];
48: #endif
50: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
51: PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
52: PetscOptionsGetInt(NULL,NULL,"-t",&ntimes,NULL);
53: PetscOptionsGetBool(NULL,NULL,"-unsym",&unsym,NULL);
54: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
55: MPI_Comm_size(PETSC_COMM_WORLD,&size);
56: n = 2*size;
58: /*
59: Register various stages for profiling
60: */
61: PetscLogStageRegister("Prelim setup",&stages[0]);
62: PetscLogStageRegister("Linear System 1",&stages[1]);
63: PetscLogStageRegister("Linear System 2",&stages[2]);
65: /*
66: Register a user-defined event for profiling (error checking).
67: */
68: CHECK_ERROR = 0;
69: PetscLogEventRegister("Check Error",KSP_CLASSID,&CHECK_ERROR);
71: /* - - - - - - - - - - - - Stage 0: - - - - - - - - - - - - - -
72: Preliminary Setup
73: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
75: PetscLogStagePush(stages[0]);
77: /*
78: Create data structures for first linear system.
79: - Create parallel matrix, specifying only its global dimensions.
80: When using MatCreate(), the matrix format can be specified at
81: runtime. Also, the parallel partitioning of the matrix is
82: determined by PETSc at runtime.
83: - Create parallel vectors.
84: - When using VecSetSizes(), we specify only the vector's global
85: dimension; the parallel partitioning is determined at runtime.
86: - Note: We form 1 vector from scratch and then duplicate as needed.
87: */
88: MatCreate(PETSC_COMM_WORLD,&C1);
89: MatSetSizes(C1,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);
90: MatSetFromOptions(C1);
91: MatSetUp(C1);
92: MatGetOwnershipRange(C1,&Istart,&Iend);
93: VecCreate(PETSC_COMM_WORLD,&u);
94: VecSetSizes(u,PETSC_DECIDE,m*n);
95: VecSetFromOptions(u);
96: VecDuplicate(u,&b1);
97: VecDuplicate(u,&x1);
99: /*
100: Create first linear solver context.
101: Set runtime options (e.g., -pc_type <type>).
102: Note that the first linear system uses the default option
103: names, while the second linear system uses a different
104: options prefix.
105: */
106: KSPCreate(PETSC_COMM_WORLD,&ksp1);
107: KSPSetFromOptions(ksp1);
109: /*
110: Set user-defined monitoring routine for first linear system.
111: */
112: PetscOptionsGetBool(NULL,NULL,"-my_ksp_monitor",&flg,NULL);
113: if (flg) {KSPMonitorSet(ksp1,MyKSPMonitor,NULL,0);}
115: /*
116: Create data structures for second linear system.
117: */
118: MatCreate(PETSC_COMM_WORLD,&C2);
119: MatSetSizes(C2,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);
120: MatSetFromOptions(C2);
121: MatSetUp(C2);
122: MatGetOwnershipRange(C2,&Istart2,&Iend2);
123: VecDuplicate(u,&b2);
124: VecDuplicate(u,&x2);
126: /*
127: Create second linear solver context
128: */
129: KSPCreate(PETSC_COMM_WORLD,&ksp2);
131: /*
132: Set different options prefix for second linear system.
133: Set runtime options (e.g., -s2_pc_type <type>)
134: */
135: KSPAppendOptionsPrefix(ksp2,"s2_");
136: KSPSetFromOptions(ksp2);
138: /*
139: Assemble exact solution vector in parallel. Note that each
140: processor needs to set only its local part of the vector.
141: */
142: VecGetLocalSize(u,&ldim);
143: VecGetOwnershipRange(u,&low,&high);
144: for (i=0; i<ldim; i++) {
145: iglobal = i + low;
146: v = (PetscScalar)(i + 100*rank);
147: VecSetValues(u,1,&iglobal,&v,ADD_VALUES);
148: }
149: VecAssemblyBegin(u);
150: VecAssemblyEnd(u);
152: /*
153: Log the number of flops for computing vector entries
154: */
155: PetscLogFlops(2.0*ldim);
157: /*
158: End curent profiling stage
159: */
160: PetscLogStagePop();
162: /* --------------------------------------------------------------
163: Linear solver loop:
164: Solve 2 different linear systems several times in succession
165: -------------------------------------------------------------- */
167: for (t=0; t<ntimes; t++) {
169: /* - - - - - - - - - - - - Stage 1: - - - - - - - - - - - - - -
170: Assemble and solve first linear system
171: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
173: /*
174: Begin profiling stage #1
175: */
176: PetscLogStagePush(stages[1]);
178: /*
179: Initialize all matrix entries to zero. MatZeroEntries() retains
180: the nonzero structure of the matrix for sparse formats.
181: */
182: if (t > 0) {MatZeroEntries(C1);}
184: /*
185: Set matrix entries in parallel. Also, log the number of flops
186: for computing matrix entries.
187: - Each processor needs to insert only elements that it owns
188: locally (but any non-local elements will be sent to the
189: appropriate processor during matrix assembly).
190: - Always specify global row and columns of matrix entries.
191: */
192: for (Ii=Istart; Ii<Iend; Ii++) {
193: v = -1.0; i = Ii/n; j = Ii - i*n;
194: if (i>0) {J = Ii - n; MatSetValues(C1,1,&Ii,1,&J,&v,ADD_VALUES);}
195: if (i<m-1) {J = Ii + n; MatSetValues(C1,1,&Ii,1,&J,&v,ADD_VALUES);}
196: if (j>0) {J = Ii - 1; MatSetValues(C1,1,&Ii,1,&J,&v,ADD_VALUES);}
197: if (j<n-1) {J = Ii + 1; MatSetValues(C1,1,&Ii,1,&J,&v,ADD_VALUES);}
198: v = 4.0; MatSetValues(C1,1,&Ii,1,&Ii,&v,ADD_VALUES);
199: }
200: if (unsym) {
201: for (Ii=Istart; Ii<Iend; Ii++) { /* Make matrix nonsymmetric */
202: v = -1.0*(t+0.5); i = Ii/n;
203: if (i>0) {J = Ii - n; MatSetValues(C1,1,&Ii,1,&J,&v,ADD_VALUES);}
204: }
205: PetscLogFlops(2.0*(Iend-Istart));
206: }
208: /*
209: Assemble matrix, using the 2-step process:
210: MatAssemblyBegin(), MatAssemblyEnd()
211: Computations can be done while messages are in transition
212: by placing code between these two statements.
213: */
214: MatAssemblyBegin(C1,MAT_FINAL_ASSEMBLY);
215: MatAssemblyEnd(C1,MAT_FINAL_ASSEMBLY);
217: /*
218: Indicate same nonzero structure of successive linear system matrices
219: */
220: MatSetOption(C1,MAT_NEW_NONZERO_LOCATIONS,PETSC_TRUE);
222: /*
223: Compute right-hand-side vector
224: */
225: MatMult(C1,u,b1);
227: /*
228: Set operators. Here the matrix that defines the linear system
229: also serves as the preconditioning matrix.
230: */
231: KSPSetOperators(ksp1,C1,C1);
233: /*
234: Use the previous solution of linear system #1 as the initial
235: guess for the next solve of linear system #1. The user MUST
236: call KSPSetInitialGuessNonzero() in indicate use of an initial
237: guess vector; otherwise, an initial guess of zero is used.
238: */
239: if (t>0) {
240: KSPSetInitialGuessNonzero(ksp1,PETSC_TRUE);
241: }
243: /*
244: Solve the first linear system. Here we explicitly call
245: KSPSetUp() for more detailed performance monitoring of
246: certain preconditioners, such as ICC and ILU. This call
247: is optional, ase KSPSetUp() will automatically be called
248: within KSPSolve() if it hasn't been called already.
249: */
250: KSPSetUp(ksp1);
251: KSPSolve(ksp1,b1,x1);
252: KSPGetIterationNumber(ksp1,&its);
254: /*
255: Check error of solution to first linear system
256: */
257: CheckError(u,x1,b1,its,1.e-4,CHECK_ERROR);
259: /* - - - - - - - - - - - - Stage 2: - - - - - - - - - - - - - -
260: Assemble and solve second linear system
261: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
263: /*
264: Conclude profiling stage #1; begin profiling stage #2
265: */
266: PetscLogStagePop();
267: PetscLogStagePush(stages[2]);
269: /*
270: Initialize all matrix entries to zero
271: */
272: if (t > 0) {MatZeroEntries(C2);}
274: /*
275: Assemble matrix in parallel. Also, log the number of flops
276: for computing matrix entries.
277: - To illustrate the features of parallel matrix assembly, we
278: intentionally set the values differently from the way in
279: which the matrix is distributed across the processors. Each
280: entry that is not owned locally will be sent to the appropriate
281: processor during MatAssemblyBegin() and MatAssemblyEnd().
282: - For best efficiency the user should strive to set as many
283: entries locally as possible.
284: */
285: for (i=0; i<m; i++) {
286: for (j=2*rank; j<2*rank+2; j++) {
287: v = -1.0; Ii = j + n*i;
288: if (i>0) {J = Ii - n; MatSetValues(C2,1,&Ii,1,&J,&v,ADD_VALUES);}
289: if (i<m-1) {J = Ii + n; MatSetValues(C2,1,&Ii,1,&J,&v,ADD_VALUES);}
290: if (j>0) {J = Ii - 1; MatSetValues(C2,1,&Ii,1,&J,&v,ADD_VALUES);}
291: if (j<n-1) {J = Ii + 1; MatSetValues(C2,1,&Ii,1,&J,&v,ADD_VALUES);}
292: v = 6.0 + t*0.5; MatSetValues(C2,1,&Ii,1,&Ii,&v,ADD_VALUES);
293: }
294: }
295: if (unsym) {
296: for (Ii=Istart2; Ii<Iend2; Ii++) { /* Make matrix nonsymmetric */
297: v = -1.0*(t+0.5); i = Ii/n;
298: if (i>0) {J = Ii - n; MatSetValues(C2,1,&Ii,1,&J,&v,ADD_VALUES);}
299: }
300: }
301: MatAssemblyBegin(C2,MAT_FINAL_ASSEMBLY);
302: MatAssemblyEnd(C2,MAT_FINAL_ASSEMBLY);
303: PetscLogFlops(2.0*(Iend-Istart));
305: /*
306: Indicate same nonzero structure of successive linear system matrices
307: */
308: MatSetOption(C2,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);
310: /*
311: Compute right-hand-side vector
312: */
313: MatMult(C2,u,b2);
315: /*
316: Set operators. Here the matrix that defines the linear system
317: also serves as the preconditioning matrix. Indicate same nonzero
318: structure of successive preconditioner matrices by setting flag
319: SAME_NONZERO_PATTERN.
320: */
321: KSPSetOperators(ksp2,C2,C2);
323: /*
324: Solve the second linear system
325: */
326: KSPSetUp(ksp2);
327: KSPSolve(ksp2,b2,x2);
328: KSPGetIterationNumber(ksp2,&its);
330: /*
331: Check error of solution to second linear system
332: */
333: CheckError(u,x2,b2,its,1.e-4,CHECK_ERROR);
335: /*
336: Conclude profiling stage #2
337: */
338: PetscLogStagePop();
339: }
340: /* --------------------------------------------------------------
341: End of linear solver loop
342: -------------------------------------------------------------- */
344: /*
345: Free work space. All PETSc objects should be destroyed when they
346: are no longer needed.
347: */
348: KSPDestroy(&ksp1); KSPDestroy(&ksp2);
349: VecDestroy(&x1); VecDestroy(&x2);
350: VecDestroy(&b1); VecDestroy(&b2);
351: MatDestroy(&C1); MatDestroy(&C2);
352: VecDestroy(&u);
354: PetscFinalize();
355: return ierr;
356: }
357: /* ------------------------------------------------------------- */
358: /*
359: CheckError - Checks the error of the solution.
361: Input Parameters:
362: u - exact solution
363: x - approximate solution
364: b - work vector
365: its - number of iterations for convergence
366: tol - tolerance
367: CHECK_ERROR - the event number for error checking
368: (for use with profiling)
370: Notes:
371: In order to profile this section of code separately from the
372: rest of the program, we register it as an "event" with
373: PetscLogEventRegister() in the main program. Then, we indicate
374: the start and end of this event by respectively calling
375: PetscLogEventBegin(CHECK_ERROR,u,x,b,0);
376: PetscLogEventEnd(CHECK_ERROR,u,x,b,0);
377: Here, we specify the objects most closely associated with
378: the event (the vectors u,x,b). Such information is optional;
379: we could instead just use 0 instead for all objects.
380: */
381: PetscErrorCode CheckError(Vec u,Vec x,Vec b,PetscInt its,PetscReal tol,PetscLogEvent CHECK_ERROR)
382: {
383: PetscScalar none = -1.0;
384: PetscReal norm;
387: PetscLogEventBegin(CHECK_ERROR,u,x,b,0);
389: /*
390: Compute error of the solution, using b as a work vector.
391: */
392: VecCopy(x,b);
393: VecAXPY(b,none,u);
394: VecNorm(b,NORM_2,&norm);
395: if (norm > tol) {
396: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g, Iterations %D\n",(double)norm,its);
397: }
398: PetscLogEventEnd(CHECK_ERROR,u,x,b,0);
399: return 0;
400: }
401: /* ------------------------------------------------------------- */
402: /*
403: MyKSPMonitor - This is a user-defined routine for monitoring
404: the KSP iterative solvers.
406: Input Parameters:
407: ksp - iterative context
408: n - iteration number
409: rnorm - 2-norm (preconditioned) residual value (may be estimated)
410: dummy - optional user-defined monitor context (unused here)
411: */
412: PetscErrorCode MyKSPMonitor(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
413: {
414: Vec x;
417: /*
418: Build the solution vector
419: */
420: KSPBuildSolution(ksp,NULL,&x);
422: /*
423: Write the solution vector and residual norm to stdout.
424: - PetscPrintf() handles output for multiprocessor jobs
425: by printing from only one processor in the communicator.
426: - The parallel viewer PETSC_VIEWER_STDOUT_WORLD handles
427: data from multiple processors so that the output
428: is not jumbled.
429: */
430: PetscPrintf(PETSC_COMM_WORLD,"iteration %D solution vector:\n",n);
431: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
432: PetscPrintf(PETSC_COMM_WORLD,"iteration %D KSP Residual norm %14.12e \n",n,rnorm);
433: return 0;
434: }
438: /*TEST
440: test:
441: args: -t 2 -pc_type jacobi -ksp_monitor_short -ksp_type gmres -ksp_gmres_cgs_refinement_type refine_always -s2_ksp_type bcgs -s2_pc_type jacobi -s2_ksp_monitor_short
443: test:
444: requires: hpddm
445: suffix: hpddm
446: args: -t 2 -pc_type jacobi -ksp_monitor_short -ksp_type {{gmres hpddm}} -s2_ksp_type {{gmres hpddm}} -s2_pc_type jacobi -s2_ksp_monitor_short
448: test:
449: requires: hpddm
450: suffix: hpddm_2
451: args: -t 2 -pc_type jacobi -ksp_monitor_short -ksp_type gmres -s2_ksp_type hpddm -s2_ksp_hpddm_type gcrodr -s2_ksp_hpddm_recycle 10 -s2_pc_type jacobi -s2_ksp_monitor_short
453: testset:
454: requires: hpddm
455: output_file: output/ex9_hpddm_cg.out
456: args: -unsym 0 -t 2 -pc_type jacobi -ksp_monitor_short -s2_pc_type jacobi -s2_ksp_monitor_short -ksp_rtol 1.e-2 -s2_ksp_rtol 1.e-2
457: test:
458: suffix: hpddm_cg_p_p
459: args: -ksp_type cg -s2_ksp_type cg
460: test:
461: suffix: hpddm_cg_p_h
462: args: -ksp_type cg -s2_ksp_type hpddm -s2_ksp_hpddm_type cg
463: test:
464: suffix: hpddm_cg_h_h
465: args: -ksp_type hpddm -ksp_hpddm_type cg -s2_ksp_type hpddm -s2_ksp_hpddm_type cg
467: TEST*/