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 repeatedn
5: solution of linear systems, while reusing matrix, vector, and solver datan
6: structures throughout the process. Note the various stages of event logging.nn";
8: /*T
9: Concepts: SLES^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 "petscsles.h" so that we can use SLES 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 petscsles.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*);
31: #undef __FUNCT__
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: SLES sles1,sles2; /* SLES contexts for systems #1 and #2 */
39: KSP ksp1; /* KSP context for system #1 */
40: int ntimes = 3; /* number of times to solve the linear systems */
41: int CHECK_ERROR; /* event number for error checking */
42: int ldim,ierr,low,high,iglobal,Istart,Iend,Istart2,Iend2;
43: int I,J,i,j,m = 3,n = 2,rank,size,its,t;
44: int stages[3];
45: PetscTruth flg;
46: PetscScalar v;
48: PetscInitialize(&argc,&args,(char *)0,help);
49: PetscOptionsGetInt(PETSC_NULL,"-m",&m,PETSC_NULL);
50: PetscOptionsGetInt(PETSC_NULL,"-t",&ntimes,PETSC_NULL);
51: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
52: MPI_Comm_size(PETSC_COMM_WORLD,&size);
53: n = 2*size;
55: /*
56: Register various stages for profiling
57: */
58: PetscLogStageRegister(&stages[0],"Prelim setup");
59: PetscLogStageRegister(&stages[1],"Linear System 1");
60: PetscLogStageRegister(&stages[2],"Linear System 2");
62: /*
63: Register a user-defined event for profiling (error checking).
64: */
65: CHECK_ERROR = 0;
66: PetscLogEventRegister(&CHECK_ERROR,"Check Error",SLES_COOKIE);
68: /* - - - - - - - - - - - - Stage 0: - - - - - - - - - - - - - -
69: Preliminary Setup
70: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
72: PetscLogStagePush(stages[0]);
74: /*
75: Create data structures for first linear system.
76: - Create parallel matrix, specifying only its global dimensions.
77: When using MatCreate(), the matrix format can be specified at
78: runtime. Also, the parallel partitioning of the matrix is
79: determined by PETSc at runtime.
80: - Create parallel vectors.
81: - When using VecSetSizes(), we specify only the vector's global
82: dimension; the parallel partitioning is determined at runtime.
83: - Note: We form 1 vector from scratch and then duplicate as needed.
84: */
85: MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n,&C1);
86: MatSetFromOptions(C1);
87: MatGetOwnershipRange(C1,&Istart,&Iend);
88: VecCreate(PETSC_COMM_WORLD,&u);
89: VecSetSizes(u,PETSC_DECIDE,m*n);
90: VecSetFromOptions(u);
91: VecDuplicate(u,&b1);
92: VecDuplicate(u,&x1);
94: /*
95: Create first linear solver context.
96: Set runtime options (e.g., -pc_type <type>).
97: Note that the first linear system uses the default option
98: names, while the second linear systme uses a different
99: options prefix.
100: */
101: SLESCreate(PETSC_COMM_WORLD,&sles1);
102: SLESSetFromOptions(sles1);
104: /*
105: Set user-defined monitoring routine for first linear system.
106: */
107: SLESGetKSP(sles1,&ksp1);
108: PetscOptionsHasName(PETSC_NULL,"-my_ksp_monitor",&flg);
109: if (flg) {KSPSetMonitor(ksp1,MyKSPMonitor,PETSC_NULL,0);}
111: /*
112: Create data structures for second linear system.
113: */
114: MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n,&C2);
115: MatSetFromOptions(C2);
116: MatGetOwnershipRange(C2,&Istart2,&Iend2);
117: VecDuplicate(u,&b2);
118: VecDuplicate(u,&x2);
120: /*
121: Create second linear solver context
122: */
123: SLESCreate(PETSC_COMM_WORLD,&sles2);
125: /*
126: Set different options prefix for second linear system.
127: Set runtime options (e.g., -s2_pc_type <type>)
128: */
129: SLESAppendOptionsPrefix(sles2,"s2_");
130: SLESSetFromOptions(sles2);
132: /*
133: Assemble exact solution vector in parallel. Note that each
134: processor needs to set only its local part of the vector.
135: */
136: VecGetLocalSize(u,&ldim);
137: VecGetOwnershipRange(u,&low,&high);
138: for (i=0; i<ldim; i++) {
139: iglobal = i + low;
140: v = (PetscScalar)(i + 100*rank);
141: VecSetValues(u,1,&iglobal,&v,ADD_VALUES);
142: }
143: VecAssemblyBegin(u);
144: VecAssemblyEnd(u);
146: /*
147: Log the number of flops for computing vector entries
148: */
149: PetscLogFlops(2*ldim);
151: /*
152: End curent profiling stage
153: */
154: PetscLogStagePop();
156: /* --------------------------------------------------------------
157: Linear solver loop:
158: Solve 2 different linear systems several times in succession
159: -------------------------------------------------------------- */
161: for (t=0; t<ntimes; t++) {
163: /* - - - - - - - - - - - - Stage 1: - - - - - - - - - - - - - -
164: Assemble and solve first linear system
165: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
167: /*
168: Begin profiling stage #1
169: */
170: PetscLogStagePush(stages[1]);
172: /*
173: Initialize all matrix entries to zero. MatZeroEntries() retains
174: the nonzero structure of the matrix for sparse formats.
175: */
176: MatZeroEntries(C1);
178: /*
179: Set matrix entries in parallel. Also, log the number of flops
180: for computing matrix entries.
181: - Each processor needs to insert only elements that it owns
182: locally (but any non-local elements will be sent to the
183: appropriate processor during matrix assembly).
184: - Always specify global row and columns of matrix entries.
185: */
186: for (I=Istart; I<Iend; I++) {
187: v = -1.0; i = I/n; j = I - i*n;
188: if (i>0) {J = I - n; MatSetValues(C1,1,&I,1,&J,&v,ADD_VALUES);}
189: if (i<m-1) {J = I + n; MatSetValues(C1,1,&I,1,&J,&v,ADD_VALUES);}
190: if (j>0) {J = I - 1; MatSetValues(C1,1,&I,1,&J,&v,ADD_VALUES);}
191: if (j<n-1) {J = I + 1; MatSetValues(C1,1,&I,1,&J,&v,ADD_VALUES);}
192: v = 4.0; MatSetValues(C1,1,&I,1,&I,&v,ADD_VALUES);
193: }
194: for (I=Istart; I<Iend; I++) { /* Make matrix nonsymmetric */
195: v = -1.0*(t+0.5); i = I/n;
196: if (i>0) {J = I - n; MatSetValues(C1,1,&I,1,&J,&v,ADD_VALUES);}
197: }
198: PetscLogFlops(2*(Istart-Iend));
200: /*
201: Assemble matrix, using the 2-step process:
202: MatAssemblyBegin(), MatAssemblyEnd()
203: Computations can be done while messages are in transition
204: by placing code between these two statements.
205: */
206: MatAssemblyBegin(C1,MAT_FINAL_ASSEMBLY);
207: MatAssemblyEnd(C1,MAT_FINAL_ASSEMBLY);
209: /*
210: Indicate same nonzero structure of successive linear system matrices
211: */
212: MatSetOption(C1,MAT_NO_NEW_NONZERO_LOCATIONS);
214: /*
215: Compute right-hand-side vector
216: */
217: MatMult(C1,u,b1);
219: /*
220: Set operators. Here the matrix that defines the linear system
221: also serves as the preconditioning matrix.
222: - The flag SAME_NONZERO_PATTERN indicates that the
223: preconditioning matrix has identical nonzero structure
224: as during the last linear solve (although the values of
225: the entries have changed). Thus, we can save some
226: work in setting up the preconditioner (e.g., no need to
227: redo symbolic factorization for ILU/ICC preconditioners).
228: - If the nonzero structure of the matrix is different during
229: the second linear solve, then the flag DIFFERENT_NONZERO_PATTERN
230: must be used instead. If you are unsure whether the
231: matrix structure has changed or not, use the flag
232: DIFFERENT_NONZERO_PATTERN.
233: - Caution: If you specify SAME_NONZERO_PATTERN, PETSc
234: believes your assertion and does not check the structure
235: of the matrix. If you erroneously claim that the structure
236: is the same when it actually is not, the new preconditioner
237: will not function correctly. Thus, use this optimization
238: feature with caution!
239: */
240: SLESSetOperators(sles1,C1,C1,SAME_NONZERO_PATTERN);
242: /*
243: Use the previous solution of linear system #1 as the initial
244: guess for the next solve of linear system #1. The user MUST
245: call KSPSetInitialGuessNonzero() in indicate use of an initial
246: guess vector; otherwise, an initial guess of zero is used.
247: */
248: if (t>0) {
249: KSPSetInitialGuessNonzero(ksp1,PETSC_TRUE);
250: }
252: /*
253: Solve the first linear system. Here we explicitly call
254: SLESSetUp() for more detailed performance monitoring of
255: certain preconditioners, such as ICC and ILU. This call
256: is optional, ase SLESSetUp() will automatically be called
257: within SLESSolve() if it hasn't been called already.
258: */
259: SLESSetUp(sles1,b1,x1);
260: SLESSolve(sles1,b1,x1,&its);
262: /*
263: Check error of solution to first linear system
264: */
265: CheckError(u,x1,b1,its,CHECK_ERROR);
267: /* - - - - - - - - - - - - Stage 2: - - - - - - - - - - - - - -
268: Assemble and solve second linear system
269: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
271: /*
272: Conclude profiling stage #1; begin profiling stage #2
273: */
274: PetscLogStagePop();
275: PetscLogStagePush(stages[2]);
277: /*
278: Initialize all matrix entries to zero
279: */
280: MatZeroEntries(C2);
282: /*
283: Assemble matrix in parallel. Also, log the number of flops
284: for computing matrix entries.
285: - To illustrate the features of parallel matrix assembly, we
286: intentionally set the values differently from the way in
287: which the matrix is distributed across the processors. Each
288: entry that is not owned locally will be sent to the appropriate
289: processor during MatAssemblyBegin() and MatAssemblyEnd().
290: - For best efficiency the user should strive to set as many
291: entries locally as possible.
292: */
293: for (i=0; i<m; i++) {
294: for (j=2*rank; j<2*rank+2; j++) {
295: v = -1.0; I = j + n*i;
296: if (i>0) {J = I - n; MatSetValues(C2,1,&I,1,&J,&v,ADD_VALUES);}
297: if (i<m-1) {J = I + n; MatSetValues(C2,1,&I,1,&J,&v,ADD_VALUES);}
298: if (j>0) {J = I - 1; MatSetValues(C2,1,&I,1,&J,&v,ADD_VALUES);}
299: if (j<n-1) {J = I + 1; MatSetValues(C2,1,&I,1,&J,&v,ADD_VALUES);}
300: v = 6.0 + t*0.5; MatSetValues(C2,1,&I,1,&I,&v,ADD_VALUES);
301: }
302: }
303: for (I=Istart2; I<Iend2; I++) { /* Make matrix nonsymmetric */
304: v = -1.0*(t+0.5); i = I/n;
305: if (i>0) {J = I - n; MatSetValues(C2,1,&I,1,&J,&v,ADD_VALUES);}
306: }
307: MatAssemblyBegin(C2,MAT_FINAL_ASSEMBLY);
308: MatAssemblyEnd(C2,MAT_FINAL_ASSEMBLY);
309: PetscLogFlops(2*(Istart-Iend));
311: /*
312: Indicate same nonzero structure of successive linear system matrices
313: */
314: MatSetOption(C2,MAT_NO_NEW_NONZERO_LOCATIONS);
316: /*
317: Compute right-hand-side vector
318: */
319: MatMult(C2,u,b2);
321: /*
322: Set operators. Here the matrix that defines the linear system
323: also serves as the preconditioning matrix. Indicate same nonzero
324: structure of successive preconditioner matrices by setting flag
325: SAME_NONZERO_PATTERN.
326: */
327: SLESSetOperators(sles2,C2,C2,SAME_NONZERO_PATTERN);
329: /*
330: Solve the second linear system
331: */
332: SLESSetUp(sles2,b2,x2);
333: SLESSolve(sles2,b2,x2,&its);
335: /*
336: Check error of solution to second linear system
337: */
338: CheckError(u,x2,b2,its,CHECK_ERROR);
340: /*
341: Conclude profiling stage #2
342: */
343: PetscLogStagePop();
344: }
345: /* --------------------------------------------------------------
346: End of linear solver loop
347: -------------------------------------------------------------- */
349: /*
350: Free work space. All PETSc objects should be destroyed when they
351: are no longer needed.
352: */
353: SLESDestroy(sles1); SLESDestroy(sles2);
354: VecDestroy(x1); VecDestroy(x2);
355: VecDestroy(b1); VecDestroy(b2);
356: MatDestroy(C1); MatDestroy(C2);
357: VecDestroy(u);
359: PetscFinalize();
360: return 0;
361: }
362: #undef __FUNCT__
364: /* ------------------------------------------------------------- */
365: /*
366: CheckError - Checks the error of the solution.
368: Input Parameters:
369: u - exact solution
370: x - approximate solution
371: b - work vector
372: its - number of iterations for convergence
373: CHECK_ERROR - the event number for error checking
374: (for use with profiling)
376: Notes:
377: In order to profile this section of code separately from the
378: rest of the program, we register it as an "event" with
379: PetscLogEventRegister() in the main program. Then, we indicate
380: the start and end of this event by respectively calling
381: PetscLogEventBegin(CHECK_ERROR,u,x,b,0);
382: PetscLogEventEnd(CHECK_ERROR,u,x,b,0);
383: Here, we specify the objects most closely associated with
384: the event (the vectors u,x,b). Such information is optional;
385: we could instead just use 0 instead for all objects.
386: */
387: int CheckError(Vec u,Vec x,Vec b,int its,int CHECK_ERROR)
388: {
389: PetscScalar none = -1.0;
390: PetscReal norm;
391: int ierr;
393: PetscLogEventBegin(CHECK_ERROR,u,x,b,0);
395: /*
396: Compute error of the solution, using b as a work vector.
397: */
398: VecCopy(x,b);
399: VecAXPY(&none,u,b);
400: VecNorm(b,NORM_2,&norm);
401: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %A, Iterations %dn",norm,its);
402: PetscLogEventEnd(CHECK_ERROR,u,x,b,0);
403: return 0;
404: }
405: /* ------------------------------------------------------------- */
406: #undef __FUNCT__
408: /*
409: MyKSPMonitor - This is a user-defined routine for monitoring
410: the SLES iterative solvers.
412: Input Parameters:
413: ksp - iterative context
414: n - iteration number
415: rnorm - 2-norm (preconditioned) residual value (may be estimated)
416: dummy - optional user-defined monitor context (unused here)
417: */
418: int MyKSPMonitor(KSP ksp,int n,PetscReal rnorm,void *dummy)
419: {
420: Vec x;
421: int ierr;
423: /*
424: Build the solution vector
425: */
426: KSPBuildSolution(ksp,PETSC_NULL,&x);
428: /*
429: Write the solution vector and residual norm to stdout.
430: - PetscPrintf() handles output for multiprocessor jobs
431: by printing from only one processor in the communicator.
432: - The parallel viewer PETSC_VIEWER_STDOUT_WORLD handles
433: data from multiple processors so that the output
434: is not jumbled.
435: */
436: PetscPrintf(PETSC_COMM_WORLD,"iteration %d solution vector:n",n);
437: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
438: PetscPrintf(PETSC_COMM_WORLD,"iteration %d KSP Residual norm %14.12e n",n,rnorm);
439: return 0;
440: }