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