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*/
17: /*
18: Include "petscksp.h" so that we can use KSP solvers. Note that this file
19: automatically includes:
20: petscsys.h - base PETSc routines petscvec.h - vectors
21: petscmat.h - matrices
22: petscis.h - index sets petscksp.h - Krylov subspace methods
23: petscviewer.h - viewers petscpc.h - preconditioners
24: */
25: #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: PetscBool mat_nonsymmetric = PETSC_FALSE;
39: PetscBool testnewC = PETSC_FALSE,testscaledMat=PETSC_FALSE;
40: #if defined(PETSC_USE_LOG)
41: PetscLogStage stages[2];
42: #endif
44: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
45: PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
46: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
47: MPI_Comm_size(PETSC_COMM_WORLD,&size);
48: n = 2*size;
50: /*
51: Set flag if we are doing a nonsymmetric problem; the default is symmetric.
52: */
53: PetscOptionsGetBool(NULL,NULL,"-mat_nonsym",&mat_nonsymmetric,NULL);
55: PetscOptionsGetBool(NULL,NULL,"-test_scaledMat",&testscaledMat,NULL);
57: /*
58: Register two stages for separate profiling of the two linear solves.
59: Use the runtime option -log_view for a printout of performance
60: statistics at the program's conlusion.
61: */
62: PetscLogStageRegister("Original Solve",&stages[0]);
63: PetscLogStageRegister("Second Solve",&stages[1]);
65: /* -------------- Stage 0: Solve Original System ---------------------- */
66: /*
67: Indicate to PETSc profiling that we're beginning the first stage
68: */
69: PetscLogStagePush(stages[0]);
71: /*
72: Create parallel matrix, specifying only its global dimensions.
73: When using MatCreate(), the matrix format can be specified at
74: runtime. Also, the parallel partitioning of the matrix is
75: determined by PETSc at runtime.
76: */
77: MatCreate(PETSC_COMM_WORLD,&C);
78: MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);
79: MatSetFromOptions(C);
80: MatSetUp(C);
82: /*
83: Currently, all PETSc parallel matrix formats are partitioned by
84: contiguous chunks of rows across the processors. Determine which
85: rows of the matrix are locally owned.
86: */
87: MatGetOwnershipRange(C,&Istart,&Iend);
89: /*
90: Set matrix entries matrix in parallel.
91: - Each processor needs to insert only elements that it owns
92: locally (but any non-local elements will be sent to the
93: appropriate processor during matrix assembly).
94: - Always specify global row and columns of matrix entries.
95: */
96: for (Ii=Istart; Ii<Iend; Ii++) {
97: v = -1.0; i = Ii/n; j = Ii - i*n;
98: if (i>0) {J = Ii - n; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
99: if (i<m-1) {J = Ii + n; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
100: if (j>0) {J = Ii - 1; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
101: if (j<n-1) {J = Ii + 1; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
102: v = 4.0; MatSetValues(C,1,&Ii,1,&Ii,&v,ADD_VALUES);
103: }
105: /*
106: Make the matrix nonsymmetric if desired
107: */
108: if (mat_nonsymmetric) {
109: for (Ii=Istart; Ii<Iend; Ii++) {
110: v = -1.5; i = Ii/n;
111: if (i>1) {J = Ii-n-1; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
112: }
113: } else {
114: MatSetOption(C,MAT_SYMMETRIC,PETSC_TRUE);
115: MatSetOption(C,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);
116: }
118: /*
119: Assemble matrix, using the 2-step process:
120: MatAssemblyBegin(), MatAssemblyEnd()
121: Computations can be done while messages are in transition
122: by placing code between these two statements.
123: */
124: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
125: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
127: /*
128: Create parallel vectors.
129: - When using VecSetSizes(), we specify only the vector's global
130: dimension; the parallel partitioning is determined at runtime.
131: - Note: We form 1 vector from scratch and then duplicate as needed.
132: */
133: VecCreate(PETSC_COMM_WORLD,&u);
134: VecSetSizes(u,PETSC_DECIDE,m*n);
135: VecSetFromOptions(u);
136: VecDuplicate(u,&b);
137: VecDuplicate(b,&x);
139: /*
140: Currently, all parallel PETSc vectors are partitioned by
141: contiguous chunks across the processors. Determine which
142: range of entries are locally owned.
143: */
144: VecGetOwnershipRange(x,&low,&high);
146: /*
147: Set elements within the exact solution vector in parallel.
148: - Each processor needs to insert only elements that it owns
149: locally (but any non-local entries will be sent to the
150: appropriate processor during vector assembly).
151: - Always specify global locations of vector entries.
152: */
153: VecGetLocalSize(x,&ldim);
154: for (i=0; i<ldim; i++) {
155: iglobal = i + low;
156: v = (PetscScalar)(i + 100*rank);
157: VecSetValues(u,1,&iglobal,&v,INSERT_VALUES);
158: }
160: /*
161: Assemble vector, using the 2-step process:
162: VecAssemblyBegin(), VecAssemblyEnd()
163: Computations can be done while messages are in transition,
164: by placing code between these two statements.
165: */
166: VecAssemblyBegin(u);
167: VecAssemblyEnd(u);
169: /*
170: Compute right-hand-side vector
171: */
172: MatMult(C,u,b);
174: /*
175: Create linear solver context
176: */
177: KSPCreate(PETSC_COMM_WORLD,&ksp);
179: /*
180: Set operators. Here the matrix that defines the linear system
181: also serves as the preconditioning matrix.
182: */
183: KSPSetOperators(ksp,C,C);
185: /*
186: Set runtime options (e.g., -ksp_type <type> -pc_type <type>)
187: */
188: KSPSetFromOptions(ksp);
190: /*
191: Solve linear system. Here we explicitly call KSPSetUp() for more
192: detailed performance monitoring of certain preconditioners, such
193: as ICC and ILU. This call is optional, as KSPSetUp() will
194: automatically be called within KSPSolve() if it hasn't been
195: called already.
196: */
197: KSPSetUp(ksp);
198: KSPSolve(ksp,b,x);
200: /*
201: Check the error
202: */
203: VecAXPY(x,none,u);
204: VecNorm(x,NORM_2,&norm);
205: KSPGetIterationNumber(ksp,&its);
206: if (!testscaledMat || norm > 1.e-7) {
207: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g, Iterations %D\n",(double)norm,its);
208: }
210: /* -------------- Stage 1: Solve Second System ---------------------- */
211: /*
212: Solve another linear system with the same method. We reuse the KSP
213: context, matrix and vector data structures, and hence save the
214: overhead of creating new ones.
216: Indicate to PETSc profiling that we're concluding the first
217: stage with PetscLogStagePop(), and beginning the second stage with
218: PetscLogStagePush().
219: */
220: PetscLogStagePop();
221: PetscLogStagePush(stages[1]);
223: /*
224: Initialize all matrix entries to zero. MatZeroEntries() retains the
225: nonzero structure of the matrix for sparse formats.
226: */
227: MatZeroEntries(C);
229: /*
230: Assemble matrix again. Note that we retain the same matrix data
231: structure and the same nonzero pattern; we just change the values
232: of the matrix entries.
233: */
234: for (i=0; i<m; i++) {
235: for (j=2*rank; j<2*rank+2; j++) {
236: v = -1.0; Ii = j + n*i;
237: if (i>0) {J = Ii - n; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
238: if (i<m-1) {J = Ii + n; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
239: if (j>0) {J = Ii - 1; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
240: if (j<n-1) {J = Ii + 1; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
241: v = 6.0; MatSetValues(C,1,&Ii,1,&Ii,&v,ADD_VALUES);
242: }
243: }
244: if (mat_nonsymmetric) {
245: for (Ii=Istart; Ii<Iend; Ii++) {
246: v = -1.5; i = Ii/n;
247: if (i>1) {J = Ii-n-1; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
248: }
249: }
250: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
251: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
253: if (testscaledMat) {
254: PetscRandom rctx;
256: /* Scale a(0,0) and a(M-1,M-1) */
257: if (!rank) {
258: v = 6.0*0.00001; Ii = 0; J = 0;
259: MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);
260: } else if (rank == size -1) {
261: v = 6.0*0.00001; Ii = m*n-1; J = m*n-1;
262: MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);
263: }
264: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
265: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
267: /* Compute a new right-hand-side vector */
268: VecDestroy(&u);
269: VecCreate(PETSC_COMM_WORLD,&u);
270: VecSetSizes(u,PETSC_DECIDE,m*n);
271: VecSetFromOptions(u);
273: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
274: PetscRandomSetFromOptions(rctx);
275: VecSetRandom(u,rctx);
276: PetscRandomDestroy(&rctx);
277: VecAssemblyBegin(u);
278: VecAssemblyEnd(u);
279: }
281: PetscOptionsGetBool(NULL,NULL,"-test_newMat",&testnewC,NULL);
282: if (testnewC) {
283: /*
284: User may use a new matrix C with same nonzero pattern, e.g.
285: ./ex5 -ksp_monitor -mat_type sbaij -pc_type cholesky -pc_factor_mat_solver_type mumps -test_newMat
286: */
287: Mat Ctmp;
288: MatDuplicate(C,MAT_COPY_VALUES,&Ctmp);
289: MatDestroy(&C);
290: MatDuplicate(Ctmp,MAT_COPY_VALUES,&C);
291: MatDestroy(&Ctmp);
292: }
294: MatMult(C,u,b);
296: /*
297: Set operators. Here the matrix that defines the linear system
298: also serves as the preconditioning matrix.
299: */
300: KSPSetOperators(ksp,C,C);
302: /*
303: Solve linear system
304: */
305: KSPSetUp(ksp);
306: KSPSolve(ksp,b,x);
308: /*
309: Check the error
310: */
311: VecAXPY(x,none,u);
312: VecNorm(x,NORM_2,&norm);
313: KSPGetIterationNumber(ksp,&its);
314: if (!testscaledMat || norm > 1.e-7) {
315: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g, Iterations %D\n",(double)norm,its);
316: }
318: /*
319: Free work space. All PETSc objects should be destroyed when they
320: are no longer needed.
321: */
322: KSPDestroy(&ksp);
323: VecDestroy(&u);
324: VecDestroy(&x);
325: VecDestroy(&b);
326: MatDestroy(&C);
328: /*
329: Indicate to PETSc profiling that we're concluding the second stage
330: */
331: PetscLogStagePop();
333: PetscFinalize();
334: return ierr;
335: }
340: /*TEST
342: test:
343: args: -pc_type jacobi -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always
345: test:
346: suffix: 2
347: nsize: 2
348: args: -pc_type jacobi -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always -ksp_rtol .000001
350: test:
351: suffix: 5
352: nsize: 2
353: args: -ksp_gmres_cgs_refinement_type refine_always -ksp_monitor draw::draw_lg -ksp_monitor_true_residual draw::draw_lg
355: test:
356: suffix: asm
357: nsize: 4
358: args: -pc_type asm
360: test:
361: suffix: asm_baij
362: nsize: 4
363: args: -pc_type asm -mat_type baij
364: output_file: output/ex5_asm.out
366: test:
367: suffix: redundant_0
368: args: -m 1000 -pc_type redundant -pc_redundant_number 1 -redundant_ksp_type gmres -redundant_pc_type jacobi
370: test:
371: suffix: redundant_1
372: nsize: 5
373: args: -pc_type redundant -pc_redundant_number 1 -redundant_ksp_type gmres -redundant_pc_type jacobi
375: test:
376: suffix: redundant_2
377: nsize: 5
378: args: -pc_type redundant -pc_redundant_number 3 -redundant_ksp_type gmres -redundant_pc_type jacobi
380: test:
381: suffix: redundant_3
382: nsize: 5
383: args: -pc_type redundant -pc_redundant_number 5 -redundant_ksp_type gmres -redundant_pc_type jacobi
385: test:
386: suffix: redundant_4
387: nsize: 5
388: args: -pc_type redundant -pc_redundant_number 3 -redundant_ksp_type gmres -redundant_pc_type jacobi -psubcomm_type interlaced
390: test:
391: suffix: superlu_dist
392: nsize: 15
393: requires: superlu_dist
394: args: -pc_type lu -pc_factor_mat_solver_type superlu_dist -mat_superlu_dist_equil false -m 150 -mat_superlu_dist_r 3 -mat_superlu_dist_c 5 -test_scaledMat
396: test:
397: suffix: superlu_dist_2
398: nsize: 15
399: requires: superlu_dist
400: args: -pc_type lu -pc_factor_mat_solver_type superlu_dist -mat_superlu_dist_equil false -m 150 -mat_superlu_dist_r 3 -mat_superlu_dist_c 5 -test_scaledMat -mat_superlu_dist_fact SamePattern_SameRowPerm
401: output_file: output/ex5_superlu_dist.out
403: test:
404: suffix: superlu_dist_3
405: nsize: 15
406: requires: superlu_dist
407: args: -pc_type lu -pc_factor_mat_solver_type superlu_dist -mat_superlu_dist_equil false -m 500 -mat_superlu_dist_r 3 -mat_superlu_dist_c 5 -test_scaledMat -mat_superlu_dist_fact DOFACT
408: output_file: output/ex5_superlu_dist.out
410: test:
411: suffix: superlu_dist_0
412: nsize: 1
413: requires: superlu_dist
414: args: -pc_type lu -pc_factor_mat_solver_type superlu_dist -test_scaledMat
415: output_file: output/ex5_superlu_dist.out
417: TEST*/