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