Actual source code: ex5.c
1: /*$Id: ex5.c,v 1.93 2001/09/11 16:33:29 bsmith Exp $*/
3: static char help[] = "Solves two linear systems in parallel with SLES. The coden
4: illustrates repeated solution of linear systems with the same preconditionern
5: method but different matrices (having the same nonzero structure). The coden
6: also uses multiple profiling stages. Input arguments aren
7: -m <size> : problem sizen
8: -mat_nonsym : use nonsymmetric matrix (default is symmetric)nn";
10: /*T
11: Concepts: SLES^repeatedly solving linear systems;
12: Concepts: PetscLog^profiling multiple stages of code;
13: Processors: n
14: T*/
16: /*
17: Include "petscsles.h" so that we can use SLES solvers. Note that this file
18: automatically includes:
19: petsc.h - base PETSc routines petscvec.h - vectors
20: petscsys.h - system routines petscmat.h - matrices
21: petscis.h - index sets petscksp.h - Krylov subspace methods
22: petscviewer.h - viewers petscpc.h - preconditioners
23: */
24: #include petscsles.h
26: #undef __FUNCT__
28: int main(int argc,char **args)
29: {
30: SLES sles; /* linear solver context */
31: Mat C; /* matrix */
32: Vec x,u,b; /* approx solution, RHS, exact solution */
33: PetscReal norm; /* norm of solution error */
34: PetscScalar v,none = -1.0;
35: int I,J,ldim,ierr,low,high,iglobal,Istart,Iend;
36: int i,j,m = 3,n = 2,rank,size,its;
37: PetscTruth mat_nonsymmetric;
38: int stages[2];
40: PetscInitialize(&argc,&args,(char *)0,help);
41: PetscOptionsGetInt(PETSC_NULL,"-m",&m,PETSC_NULL);
42: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
43: MPI_Comm_size(PETSC_COMM_WORLD,&size);
44: n = 2*size;
46: /*
47: Set flag if we are doing a nonsymmetric problem; the default is symmetric.
48: */
49: PetscOptionsHasName(PETSC_NULL,"-mat_nonsym",&mat_nonsymmetric);
51: /*
52: Register two stages for separate profiling of the two linear solves.
53: Use the runtime option -log_summary for a printout of performance
54: statistics at the program's conlusion.
55: */
56: PetscLogStageRegister(&stages[0],"Original Solve");
57: PetscLogStageRegister(&stages[1],"Second Solve");
59: /* -------------- Stage 0: Solve Original System ---------------------- */
60: /*
61: Indicate to PETSc profiling that we're beginning the first stage
62: */
63: PetscLogStagePush(stages[0]);
65: /*
66: Create parallel matrix, specifying only its global dimensions.
67: When using MatCreate(), the matrix format can be specified at
68: runtime. Also, the parallel partitioning of the matrix is
69: determined by PETSc at runtime.
70: */
71: MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n,&C);
72: MatSetFromOptions(C);
74: /*
75: Currently, all PETSc parallel matrix formats are partitioned by
76: contiguous chunks of rows across the processors. Determine which
77: rows of the matrix are locally owned.
78: */
79: MatGetOwnershipRange(C,&Istart,&Iend);
81: /*
82: Set matrix entries matrix in parallel.
83: - Each processor needs to insert only elements that it owns
84: locally (but any non-local elements will be sent to the
85: appropriate processor during matrix assembly).
86: - Always specify global row and columns of matrix entries.
87: */
88: for (I=Istart; I<Iend; I++) {
89: v = -1.0; i = I/n; j = I - i*n;
90: if (i>0) {J = I - n; MatSetValues(C,1,&I,1,&J,&v,ADD_VALUES);}
91: if (i<m-1) {J = I + n; MatSetValues(C,1,&I,1,&J,&v,ADD_VALUES);}
92: if (j>0) {J = I - 1; MatSetValues(C,1,&I,1,&J,&v,ADD_VALUES);}
93: if (j<n-1) {J = I + 1; MatSetValues(C,1,&I,1,&J,&v,ADD_VALUES);}
94: v = 4.0; MatSetValues(C,1,&I,1,&I,&v,ADD_VALUES);
95: }
97: /*
98: Make the matrix nonsymmetric if desired
99: */
100: if (mat_nonsymmetric) {
101: for (I=Istart; I<Iend; I++) {
102: v = -1.5; i = I/n;
103: if (i>1) {J = I-n-1; MatSetValues(C,1,&I,1,&J,&v,ADD_VALUES);}
104: }
105: } else {
106: MatSetOption(C,MAT_SYMMETRIC);
107: }
109: /*
110: Assemble matrix, using the 2-step process:
111: MatAssemblyBegin(), MatAssemblyEnd()
112: Computations can be done while messages are in transition
113: by placing code between these two statements.
114: */
115: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
116: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
118: /*
119: Create parallel vectors.
120: - When using VecSetSizes(), we specify only the vector's global
121: dimension; the parallel partitioning is determined at runtime.
122: - Note: We form 1 vector from scratch and then duplicate as needed.
123: */
124: VecCreate(PETSC_COMM_WORLD,&u);
125: VecSetSizes(u,PETSC_DECIDE,m*n);
126: VecSetFromOptions(u);
127: VecDuplicate(u,&b);
128: VecDuplicate(b,&x);
130: /*
131: Currently, all parallel PETSc vectors are partitioned by
132: contiguous chunks across the processors. Determine which
133: range of entries are locally owned.
134: */
135: VecGetOwnershipRange(x,&low,&high);
137: /*
138: Set elements within the exact solution vector in parallel.
139: - Each processor needs to insert only elements that it owns
140: locally (but any non-local entries will be sent to the
141: appropriate processor during vector assembly).
142: - Always specify global locations of vector entries.
143: */
144: VecGetLocalSize(x,&ldim);
145: for (i=0; i<ldim; i++) {
146: iglobal = i + low;
147: v = (PetscScalar)(i + 100*rank);
148: VecSetValues(u,1,&iglobal,&v,INSERT_VALUES);
149: }
151: /*
152: Assemble vector, using the 2-step process:
153: VecAssemblyBegin(), VecAssemblyEnd()
154: Computations can be done while messages are in transition,
155: by placing code between these two statements.
156: */
157: VecAssemblyBegin(u);
158: VecAssemblyEnd(u);
160: /*
161: Compute right-hand-side vector
162: */
163: MatMult(C,u,b);
164:
165: /*
166: Create linear solver context
167: */
168: SLESCreate(PETSC_COMM_WORLD,&sles);
170: /*
171: Set operators. Here the matrix that defines the linear system
172: also serves as the preconditioning matrix.
173: */
174: SLESSetOperators(sles,C,C,DIFFERENT_NONZERO_PATTERN);
176: /*
177: Set runtime options (e.g., -ksp_type <type> -pc_type <type>)
178: */
180: SLESSetFromOptions(sles);
182: /*
183: Solve linear system. Here we explicitly call SLESSetUp() for more
184: detailed performance monitoring of certain preconditioners, such
185: as ICC and ILU. This call is optional, as SLESSetUp() will
186: automatically be called within SLESSolve() if it hasn't been
187: called already.
188: */
189: SLESSetUp(sles,b,x);
190: SLESSolve(sles,b,x,&its);
191:
192: /*
193: Check the error
194: */
195: VecAXPY(&none,u,x);
196: VecNorm(x,NORM_2,&norm);
197: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %A, Iterations %dn",norm,its);
199: /* -------------- Stage 1: Solve Second System ---------------------- */
200: /*
201: Solve another linear system with the same method. We reuse the SLES
202: context, matrix and vector data structures, and hence save the
203: overhead of creating new ones.
205: Indicate to PETSc profiling that we're concluding the first
206: stage with PetscLogStagePop(), and beginning the second stage with
207: PetscLogStagePush().
208: */
209: PetscLogStagePop();
210: PetscLogStagePush(stages[1]);
212: /*
213: Initialize all matrix entries to zero. MatZeroEntries() retains the
214: nonzero structure of the matrix for sparse formats.
215: */
216: MatZeroEntries(C);
218: /*
219: Assemble matrix again. Note that we retain the same matrix data
220: structure and the same nonzero pattern; we just change the values
221: of the matrix entries.
222: */
223: for (i=0; i<m; i++) {
224: for (j=2*rank; j<2*rank+2; j++) {
225: v = -1.0; I = j + n*i;
226: if (i>0) {J = I - n; MatSetValues(C,1,&I,1,&J,&v,ADD_VALUES);}
227: if (i<m-1) {J = I + n; MatSetValues(C,1,&I,1,&J,&v,ADD_VALUES);}
228: if (j>0) {J = I - 1; MatSetValues(C,1,&I,1,&J,&v,ADD_VALUES);}
229: if (j<n-1) {J = I + 1; MatSetValues(C,1,&I,1,&J,&v,ADD_VALUES);}
230: v = 6.0; MatSetValues(C,1,&I,1,&I,&v,ADD_VALUES);
231: }
232: }
233: if (mat_nonsymmetric) {
234: for (I=Istart; I<Iend; I++) {
235: v = -1.5; i = I/n;
236: if (i>1) {J = I-n-1; MatSetValues(C,1,&I,1,&J,&v,ADD_VALUES);}
237: }
238: }
239: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
240: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
242: /*
243: Compute another right-hand-side vector
244: */
245: MatMult(C,u,b);
247: /*
248: Set operators. Here the matrix that defines the linear system
249: also serves as the preconditioning matrix.
250: - The flag SAME_NONZERO_PATTERN indicates that the
251: preconditioning matrix has identical nonzero structure
252: as during the last linear solve (although the values of
253: the entries have changed). Thus, we can save some
254: work in setting up the preconditioner (e.g., no need to
255: redo symbolic factorization for ILU/ICC preconditioners).
256: - If the nonzero structure of the matrix is different during
257: the second linear solve, then the flag DIFFERENT_NONZERO_PATTERN
258: must be used instead. If you are unsure whether the
259: matrix structure has changed or not, use the flag
260: DIFFERENT_NONZERO_PATTERN.
261: - Caution: If you specify SAME_NONZERO_PATTERN, PETSc
262: believes your assertion and does not check the structure
263: of the matrix. If you erroneously claim that the structure
264: is the same when it actually is not, the new preconditioner
265: will not function correctly. Thus, use this optimization
266: feature with caution!
267: */
268: SLESSetOperators(sles,C,C,SAME_NONZERO_PATTERN);
270: /*
271: Solve linear system
272: */
273: SLESSetUp(sles,b,x);
274: SLESSolve(sles,b,x,&its);
276: /*
277: Check the error
278: */
279: VecAXPY(&none,u,x);
280: VecNorm(x,NORM_2,&norm);
281: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %A, Iterations %dn",norm,its);
283: /*
284: Free work space. All PETSc objects should be destroyed when they
285: are no longer needed.
286: */
287: SLESDestroy(sles);
288: VecDestroy(u);
289: VecDestroy(x);
290: VecDestroy(b);
291: MatDestroy(C);
293: /*
294: Indicate to PETSc profiling that we're concluding the second stage
295: */
296: PetscLogStagePop();
298: PetscFinalize();
299: return 0;
300: }