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