Actual source code: ex9.c
petsc-3.4.5 2014-06-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*/
14: /*
15: Include "petscksp.h" so that we can use KSP solvers. Note that this file
16: automatically includes:
17: petscsys.h - base PETSc routines petscvec.h - vectors
18: petscmat.h - matrices
19: petscis.h - index sets petscksp.h - Krylov subspace methods
20: petscviewer.h - viewers petscpc.h - preconditioners
21: */
22: #include <petscksp.h>
24: /*
25: Declare user-defined routines
26: */
27: extern PetscErrorCode CheckError(Vec,Vec,Vec,PetscInt,PetscReal,PetscLogEvent);
28: 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;
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);
51: PetscOptionsGetInt(NULL,"-m",&m,NULL);
52: PetscOptionsGetInt(NULL,"-t",&ntimes,NULL);
53: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
54: MPI_Comm_size(PETSC_COMM_WORLD,&size);
55: n = 2*size;
57: /*
58: Register various stages for profiling
59: */
60: PetscLogStageRegister("Prelim setup",&stages[0]);
61: PetscLogStageRegister("Linear System 1",&stages[1]);
62: PetscLogStageRegister("Linear System 2",&stages[2]);
64: /*
65: Register a user-defined event for profiling (error checking).
66: */
67: CHECK_ERROR = 0;
68: PetscLogEventRegister("Check Error",KSP_CLASSID,&CHECK_ERROR);
70: /* - - - - - - - - - - - - Stage 0: - - - - - - - - - - - - - -
71: Preliminary Setup
72: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
74: PetscLogStagePush(stages[0]);
76: /*
77: Create data structures for first linear system.
78: - Create parallel matrix, specifying only its global dimensions.
79: When using MatCreate(), the matrix format can be specified at
80: runtime. Also, the parallel partitioning of the matrix is
81: determined by PETSc at runtime.
82: - Create parallel vectors.
83: - When using VecSetSizes(), we specify only the vector's global
84: dimension; the parallel partitioning is determined at runtime.
85: - Note: We form 1 vector from scratch and then duplicate as needed.
86: */
87: MatCreate(PETSC_COMM_WORLD,&C1);
88: MatSetSizes(C1,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);
89: MatSetFromOptions(C1);
90: MatSetUp(C1);
91: MatGetOwnershipRange(C1,&Istart,&Iend);
92: VecCreate(PETSC_COMM_WORLD,&u);
93: VecSetSizes(u,PETSC_DECIDE,m*n);
94: VecSetFromOptions(u);
95: VecDuplicate(u,&b1);
96: VecDuplicate(u,&x1);
98: /*
99: Create first linear solver context.
100: Set runtime options (e.g., -pc_type <type>).
101: Note that the first linear system uses the default option
102: names, while the second linear systme uses a different
103: options prefix.
104: */
105: KSPCreate(PETSC_COMM_WORLD,&ksp1);
106: KSPSetFromOptions(ksp1);
108: /*
109: Set user-defined monitoring routine for first linear system.
110: */
111: PetscOptionsGetBool(NULL,"-my_ksp_monitor",&flg,NULL);
112: if (flg) {KSPMonitorSet(ksp1,MyKSPMonitor,NULL,0);}
114: /*
115: Create data structures for second linear system.
116: */
117: MatCreate(PETSC_COMM_WORLD,&C2);
118: MatSetSizes(C2,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);
119: MatSetFromOptions(C2);
120: MatSetUp(C2);
121: MatGetOwnershipRange(C2,&Istart2,&Iend2);
122: VecDuplicate(u,&b2);
123: VecDuplicate(u,&x2);
125: /*
126: Create second linear solver context
127: */
128: KSPCreate(PETSC_COMM_WORLD,&ksp2);
130: /*
131: Set different options prefix for second linear system.
132: Set runtime options (e.g., -s2_pc_type <type>)
133: */
134: KSPAppendOptionsPrefix(ksp2,"s2_");
135: KSPSetFromOptions(ksp2);
137: /*
138: Assemble exact solution vector in parallel. Note that each
139: processor needs to set only its local part of the vector.
140: */
141: VecGetLocalSize(u,&ldim);
142: VecGetOwnershipRange(u,&low,&high);
143: for (i=0; i<ldim; i++) {
144: iglobal = i + low;
145: v = (PetscScalar)(i + 100*rank);
146: VecSetValues(u,1,&iglobal,&v,ADD_VALUES);
147: }
148: VecAssemblyBegin(u);
149: VecAssemblyEnd(u);
151: /*
152: Log the number of flops for computing vector entries
153: */
154: PetscLogFlops(2.0*ldim);
156: /*
157: End curent profiling stage
158: */
159: PetscLogStagePop();
161: /* --------------------------------------------------------------
162: Linear solver loop:
163: Solve 2 different linear systems several times in succession
164: -------------------------------------------------------------- */
166: for (t=0; t<ntimes; t++) {
168: /* - - - - - - - - - - - - Stage 1: - - - - - - - - - - - - - -
169: Assemble and solve first linear system
170: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
172: /*
173: Begin profiling stage #1
174: */
175: PetscLogStagePush(stages[1]);
177: /*
178: Initialize all matrix entries to zero. MatZeroEntries() retains
179: the nonzero structure of the matrix for sparse formats.
180: */
181: if (t > 0) {MatZeroEntries(C1);}
183: /*
184: Set matrix entries in parallel. Also, log the number of flops
185: for computing matrix entries.
186: - Each processor needs to insert only elements that it owns
187: locally (but any non-local elements will be sent to the
188: appropriate processor during matrix assembly).
189: - Always specify global row and columns of matrix entries.
190: */
191: for (Ii=Istart; Ii<Iend; Ii++) {
192: v = -1.0; i = Ii/n; j = Ii - i*n;
193: if (i>0) {J = Ii - n; MatSetValues(C1,1,&Ii,1,&J,&v,ADD_VALUES);}
194: if (i<m-1) {J = Ii + n; MatSetValues(C1,1,&Ii,1,&J,&v,ADD_VALUES);}
195: if (j>0) {J = Ii - 1; MatSetValues(C1,1,&Ii,1,&J,&v,ADD_VALUES);}
196: if (j<n-1) {J = Ii + 1; MatSetValues(C1,1,&Ii,1,&J,&v,ADD_VALUES);}
197: v = 4.0; MatSetValues(C1,1,&Ii,1,&Ii,&v,ADD_VALUES);
198: }
199: for (Ii=Istart; Ii<Iend; Ii++) { /* Make matrix nonsymmetric */
200: v = -1.0*(t+0.5); i = Ii/n;
201: if (i>0) {J = Ii - n; MatSetValues(C1,1,&Ii,1,&J,&v,ADD_VALUES);}
202: }
203: PetscLogFlops(2.0*(Iend-Istart));
205: /*
206: Assemble matrix, using the 2-step process:
207: MatAssemblyBegin(), MatAssemblyEnd()
208: Computations can be done while messages are in transition
209: by placing code between these two statements.
210: */
211: MatAssemblyBegin(C1,MAT_FINAL_ASSEMBLY);
212: MatAssemblyEnd(C1,MAT_FINAL_ASSEMBLY);
214: /*
215: Indicate same nonzero structure of successive linear system matrices
216: */
217: MatSetOption(C1,MAT_NEW_NONZERO_LOCATIONS,PETSC_TRUE);
219: /*
220: Compute right-hand-side vector
221: */
222: MatMult(C1,u,b1);
224: /*
225: Set operators. Here the matrix that defines the linear system
226: also serves as the preconditioning matrix.
227: - The flag SAME_NONZERO_PATTERN indicates that the
228: preconditioning matrix has identical nonzero structure
229: as during the last linear solve (although the values of
230: the entries have changed). Thus, we can save some
231: work in setting up the preconditioner (e.g., no need to
232: redo symbolic factorization for ILU/ICC preconditioners).
233: - If the nonzero structure of the matrix is different during
234: the second linear solve, then the flag DIFFERENT_NONZERO_PATTERN
235: must be used instead. If you are unsure whether the
236: matrix structure has changed or not, use the flag
237: DIFFERENT_NONZERO_PATTERN.
238: - Caution: If you specify SAME_NONZERO_PATTERN, PETSc
239: believes your assertion and does not check the structure
240: of the matrix. If you erroneously claim that the structure
241: is the same when it actually is not, the new preconditioner
242: will not function correctly. Thus, use this optimization
243: feature with caution!
244: */
245: KSPSetOperators(ksp1,C1,C1,SAME_NONZERO_PATTERN);
247: /*
248: Use the previous solution of linear system #1 as the initial
249: guess for the next solve of linear system #1. The user MUST
250: call KSPSetInitialGuessNonzero() in indicate use of an initial
251: guess vector; otherwise, an initial guess of zero is used.
252: */
253: if (t>0) {
254: KSPSetInitialGuessNonzero(ksp1,PETSC_TRUE);
255: }
257: /*
258: Solve the first linear system. Here we explicitly call
259: KSPSetUp() for more detailed performance monitoring of
260: certain preconditioners, such as ICC and ILU. This call
261: is optional, ase KSPSetUp() will automatically be called
262: within KSPSolve() if it hasn't been called already.
263: */
264: KSPSetUp(ksp1);
265: KSPSolve(ksp1,b1,x1);
266: KSPGetIterationNumber(ksp1,&its);
268: /*
269: Check error of solution to first linear system
270: */
271: CheckError(u,x1,b1,its,1.e-4,CHECK_ERROR);
273: /* - - - - - - - - - - - - Stage 2: - - - - - - - - - - - - - -
274: Assemble and solve second linear system
275: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
277: /*
278: Conclude profiling stage #1; begin profiling stage #2
279: */
280: PetscLogStagePop();
281: PetscLogStagePush(stages[2]);
283: /*
284: Initialize all matrix entries to zero
285: */
286: if (t > 0) {MatZeroEntries(C2);}
288: /*
289: Assemble matrix in parallel. Also, log the number of flops
290: for computing matrix entries.
291: - To illustrate the features of parallel matrix assembly, we
292: intentionally set the values differently from the way in
293: which the matrix is distributed across the processors. Each
294: entry that is not owned locally will be sent to the appropriate
295: processor during MatAssemblyBegin() and MatAssemblyEnd().
296: - For best efficiency the user should strive to set as many
297: entries locally as possible.
298: */
299: for (i=0; i<m; i++) {
300: for (j=2*rank; j<2*rank+2; j++) {
301: v = -1.0; Ii = j + n*i;
302: if (i>0) {J = Ii - n; MatSetValues(C2,1,&Ii,1,&J,&v,ADD_VALUES);}
303: if (i<m-1) {J = Ii + n; MatSetValues(C2,1,&Ii,1,&J,&v,ADD_VALUES);}
304: if (j>0) {J = Ii - 1; MatSetValues(C2,1,&Ii,1,&J,&v,ADD_VALUES);}
305: if (j<n-1) {J = Ii + 1; MatSetValues(C2,1,&Ii,1,&J,&v,ADD_VALUES);}
306: v = 6.0 + t*0.5; MatSetValues(C2,1,&Ii,1,&Ii,&v,ADD_VALUES);
307: }
308: }
309: for (Ii=Istart2; Ii<Iend2; Ii++) { /* Make matrix nonsymmetric */
310: v = -1.0*(t+0.5); i = Ii/n;
311: if (i>0) {J = Ii - n; MatSetValues(C2,1,&Ii,1,&J,&v,ADD_VALUES);}
312: }
313: MatAssemblyBegin(C2,MAT_FINAL_ASSEMBLY);
314: MatAssemblyEnd(C2,MAT_FINAL_ASSEMBLY);
315: PetscLogFlops(2.0*(Iend-Istart));
317: /*
318: Indicate same nonzero structure of successive linear system matrices
319: */
320: MatSetOption(C2,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);
322: /*
323: Compute right-hand-side vector
324: */
325: MatMult(C2,u,b2);
327: /*
328: Set operators. Here the matrix that defines the linear system
329: also serves as the preconditioning matrix. Indicate same nonzero
330: structure of successive preconditioner matrices by setting flag
331: SAME_NONZERO_PATTERN.
332: */
333: KSPSetOperators(ksp2,C2,C2,SAME_NONZERO_PATTERN);
335: /*
336: Solve the second linear system
337: */
338: KSPSetUp(ksp2);
339: KSPSolve(ksp2,b2,x2);
340: KSPGetIterationNumber(ksp2,&its);
342: /*
343: Check error of solution to second linear system
344: */
345: CheckError(u,x2,b2,its,1.e-4,CHECK_ERROR);
347: /*
348: Conclude profiling stage #2
349: */
350: PetscLogStagePop();
351: }
352: /* --------------------------------------------------------------
353: End of linear solver loop
354: -------------------------------------------------------------- */
356: /*
357: Free work space. All PETSc objects should be destroyed when they
358: are no longer needed.
359: */
360: KSPDestroy(&ksp1); KSPDestroy(&ksp2);
361: VecDestroy(&x1); VecDestroy(&x2);
362: VecDestroy(&b1); VecDestroy(&b2);
363: MatDestroy(&C1); MatDestroy(&C2);
364: VecDestroy(&u);
366: PetscFinalize();
367: return 0;
368: }
371: /* ------------------------------------------------------------- */
372: /*
373: CheckError - Checks the error of the solution.
375: Input Parameters:
376: u - exact solution
377: x - approximate solution
378: b - work vector
379: its - number of iterations for convergence
380: tol - tolerance
381: CHECK_ERROR - the event number for error checking
382: (for use with profiling)
384: Notes:
385: In order to profile this section of code separately from the
386: rest of the program, we register it as an "event" with
387: PetscLogEventRegister() in the main program. Then, we indicate
388: the start and end of this event by respectively calling
389: PetscLogEventBegin(CHECK_ERROR,u,x,b,0);
390: PetscLogEventEnd(CHECK_ERROR,u,x,b,0);
391: Here, we specify the objects most closely associated with
392: the event (the vectors u,x,b). Such information is optional;
393: we could instead just use 0 instead for all objects.
394: */
395: PetscErrorCode CheckError(Vec u,Vec x,Vec b,PetscInt its,PetscReal tol,PetscLogEvent CHECK_ERROR)
396: {
397: PetscScalar none = -1.0;
398: PetscReal norm;
401: PetscLogEventBegin(CHECK_ERROR,u,x,b,0);
403: /*
404: Compute error of the solution, using b as a work vector.
405: */
406: VecCopy(x,b);
407: VecAXPY(b,none,u);
408: VecNorm(b,NORM_2,&norm);
409: if (norm > tol) {
410: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %G, Iterations %D\n",norm,its);
411: }
412: PetscLogEventEnd(CHECK_ERROR,u,x,b,0);
413: return 0;
414: }
415: /* ------------------------------------------------------------- */
418: /*
419: MyKSPMonitor - This is a user-defined routine for monitoring
420: the KSP iterative solvers.
422: Input Parameters:
423: ksp - iterative context
424: n - iteration number
425: rnorm - 2-norm (preconditioned) residual value (may be estimated)
426: dummy - optional user-defined monitor context (unused here)
427: */
428: PetscErrorCode MyKSPMonitor(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
429: {
430: Vec x;
433: /*
434: Build the solution vector
435: */
436: KSPBuildSolution(ksp,NULL,&x);
438: /*
439: Write the solution vector and residual norm to stdout.
440: - PetscPrintf() handles output for multiprocessor jobs
441: by printing from only one processor in the communicator.
442: - The parallel viewer PETSC_VIEWER_STDOUT_WORLD handles
443: data from multiple processors so that the output
444: is not jumbled.
445: */
446: PetscPrintf(PETSC_COMM_WORLD,"iteration %D solution vector:\n",n);
447: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
448: PetscPrintf(PETSC_COMM_WORLD,"iteration %D KSP Residual norm %14.12e \n",n,rnorm);
449: return 0;
450: }