Actual source code: ex10.c
2: static char help[] = "Solve a small system and a large system through preloading\n\
3: Input arguments are:\n\
4: -permute <natural,rcm,nd,...> : solve system in permuted indexing\n\
5: -f0 <small_sys_binary> -f1 <large_sys_binary> \n\n";
7: /*
8: Include "petscksp.h" so that we can use KSP solvers. Note that this file
9: automatically includes:
10: petscsys.h - base PETSc routines petscvec.h - vectors
11: petscmat.h - matrices
12: petscis.h - index sets petscksp.h - Krylov subspace methods
13: petscviewer.h - viewers petscpc.h - preconditioners
14: */
15: #include <petscksp.h>
17: typedef enum {
18: RHS_FILE,
19: RHS_ONE,
20: RHS_RANDOM
21: } RHSType;
22: const char *const RHSTypes[] = {"FILE", "ONE", "RANDOM", "RHSType", "RHS_", NULL};
24: PetscErrorCode CheckResult(KSP *ksp, Mat *A, Vec *b, Vec *x, IS *rowperm)
25: {
26: PetscReal norm; /* norm of solution error */
27: PetscInt its;
28: KSPGetTotalIterations(*ksp, &its);
29: PetscPrintf(PETSC_COMM_WORLD, "Number of iterations = %" PetscInt_FMT "\n", its);
31: KSPGetResidualNorm(*ksp, &norm);
32: if (norm < 1.e-12) {
33: PetscPrintf(PETSC_COMM_WORLD, "Residual norm < 1.e-12\n");
34: } else {
35: PetscPrintf(PETSC_COMM_WORLD, "Residual norm %e\n", (double)norm);
36: }
38: KSPDestroy(ksp);
39: MatDestroy(A);
40: VecDestroy(x);
41: VecDestroy(b);
42: ISDestroy(rowperm);
43: return 0;
44: }
46: PetscErrorCode CreateSystem(const char filename[PETSC_MAX_PATH_LEN], RHSType rhstype, MatOrderingType ordering, PetscBool permute, IS *rowperm_out, Mat *A_out, Vec *b_out, Vec *x_out)
47: {
48: Vec x, b, b2;
49: Mat A; /* linear system matrix */
50: PetscViewer viewer; /* viewer */
51: PetscBool same;
52: PetscInt j, len, start, idx, n1, n2;
53: const PetscScalar *val;
54: IS rowperm = NULL, colperm = NULL;
56: /* open binary file. Note that we use FILE_MODE_READ to indicate reading from this file */
57: PetscViewerBinaryOpen(PETSC_COMM_WORLD, filename, FILE_MODE_READ, &viewer);
59: /* load the matrix and vector; then destroy the viewer */
60: MatCreate(PETSC_COMM_WORLD, &A);
61: MatSetFromOptions(A);
62: MatLoad(A, viewer);
63: switch (rhstype) {
64: case RHS_FILE:
65: /* Vectors in the file might a different size than the matrix so we need a
66: * Vec whose size hasn't been set yet. It'll get fixed below. Otherwise we
67: * can create the correct size Vec. */
68: VecCreate(PETSC_COMM_WORLD, &b);
69: VecLoad(b, viewer);
70: break;
71: case RHS_ONE:
72: MatCreateVecs(A, &b, NULL);
73: VecSet(b, 1.0);
74: break;
75: case RHS_RANDOM:
76: MatCreateVecs(A, &b, NULL);
77: VecSetRandom(b, NULL);
78: break;
79: }
80: PetscViewerDestroy(&viewer);
82: /* if the loaded matrix is larger than the vector (due to being padded
83: to match the block size of the system), then create a new padded vector
84: */
85: MatGetLocalSize(A, NULL, &n1);
86: VecGetLocalSize(b, &n2);
87: same = (n1 == n2) ? PETSC_TRUE : PETSC_FALSE;
88: MPIU_Allreduce(MPI_IN_PLACE, &same, 1, MPIU_BOOL, MPI_LAND, PETSC_COMM_WORLD);
90: if (!same) { /* create a new vector b by padding the old one */
91: VecCreate(PETSC_COMM_WORLD, &b2);
92: VecSetSizes(b2, n1, PETSC_DECIDE);
93: VecSetFromOptions(b2);
94: VecGetOwnershipRange(b, &start, NULL);
95: VecGetLocalSize(b, &len);
96: VecGetArrayRead(b, &val);
97: for (j = 0; j < len; j++) {
98: idx = start + j;
99: VecSetValues(b2, 1, &idx, val + j, INSERT_VALUES);
100: }
101: VecRestoreArrayRead(b, &val);
102: VecDestroy(&b);
103: VecAssemblyBegin(b2);
104: VecAssemblyEnd(b2);
105: b = b2;
106: }
107: VecDuplicate(b, &x);
109: if (permute) {
110: Mat Aperm;
111: MatGetOrdering(A, ordering, &rowperm, &colperm);
112: MatPermute(A, rowperm, colperm, &Aperm);
113: VecPermute(b, colperm, PETSC_FALSE);
114: MatDestroy(&A);
115: A = Aperm; /* Replace original operator with permuted version */
116: ISDestroy(&colperm);
117: }
119: *b_out = b;
120: *x_out = x;
121: *A_out = A;
122: *rowperm_out = rowperm;
124: return 0;
125: }
127: /* ATTENTION: this is the example used in the Profiling chaper of the PETSc manual,
128: where we referenced its profiling stages, preloading and output etc.
129: When you modify it, please make sure it is still consistent with the manual.
130: */
131: int main(int argc, char **args)
132: {
133: Vec x, b;
134: Mat A; /* linear system matrix */
135: KSP ksp; /* Krylov subspace method context */
136: char file[2][PETSC_MAX_PATH_LEN], ordering[256] = MATORDERINGRCM;
137: RHSType rhstype = RHS_FILE;
138: PetscBool flg, preload = PETSC_FALSE, trans = PETSC_FALSE, permute = PETSC_FALSE;
139: IS rowperm = NULL;
142: PetscInitialize(&argc, &args, (char *)0, help);
144: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Preloading example options", "");
145: {
146: /*
147: Determine files from which we read the two linear systems
148: (matrix and right-hand-side vector).
149: */
150: PetscOptionsBool("-trans", "Solve transpose system instead", "", trans, &trans, &flg);
151: PetscOptionsString("-f", "First file to load (small system)", "", file[0], file[0], sizeof(file[0]), &flg);
152: PetscOptionsFList("-permute", "Permute matrix and vector to solve in new ordering", "", MatOrderingList, ordering, ordering, sizeof(ordering), &permute);
154: if (flg) {
155: PetscStrcpy(file[1], file[0]);
156: preload = PETSC_FALSE;
157: } else {
158: PetscOptionsString("-f0", "First file to load (small system)", "", file[0], file[0], sizeof(file[0]), &flg);
160: PetscOptionsString("-f1", "Second file to load (larger system)", "", file[1], file[1], sizeof(file[1]), &flg);
161: if (!flg) preload = PETSC_FALSE; /* don't bother with second system */
162: }
164: PetscOptionsEnum("-rhs", "Right hand side", "", RHSTypes, (PetscEnum)rhstype, (PetscEnum *)&rhstype, NULL);
165: }
166: PetscOptionsEnd();
168: /*
169: To use preloading, one usually has code like the following:
171: PetscPreLoadBegin(preload,"first stage);
172: lines of code
173: PetscPreLoadStage("second stage");
174: lines of code
175: PetscPreLoadEnd();
177: The two macro PetscPreLoadBegin() and PetscPreLoadEnd() implicitly form a
178: loop with maximal two iterations, depending whether preloading is turned on or
179: not. If it is, either through the preload arg of PetscPreLoadBegin or through
180: -preload command line, the trip count is 2, otherwise it is 1. One can use the
181: predefined variable PetscPreLoadIt within the loop body to get the current
182: iteration number, which is 0 or 1. If preload is turned on, the runtime doesn't
183: do profiling for the first iteration, but it will do profiling for the second
184: iteration instead.
186: One can solve a small system in the first iteration and a large system in
187: the second iteration. This process preloads the instructions with the small
188: system so that more accurate performance monitoring (via -log_view) can be done
189: with the large one (that actually is the system of interest).
191: But in this example, we turned off preloading and duplicated the code for
192: the large system. In general, it is a bad practice and one should not duplicate
193: code. We do that because we want to show profiling stages for both the small
194: system and the large system.
195: */
197: /*=========================
198: solve a small system
199: =========================*/
201: PetscPreLoadBegin(preload, "Load System 0");
202: CreateSystem(file[0], rhstype, ordering, permute, &rowperm, &A, &b, &x);
204: PetscPreLoadStage("KSPSetUp 0");
205: KSPCreate(PETSC_COMM_WORLD, &ksp);
206: KSPSetOperators(ksp, A, A);
207: KSPSetFromOptions(ksp);
209: /*
210: Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
211: enable more precise profiling of setting up the preconditioner.
212: These calls are optional, since both will be called within
213: KSPSolve() if they haven't been called already.
214: */
215: KSPSetUp(ksp);
216: KSPSetUpOnBlocks(ksp);
218: PetscPreLoadStage("KSPSolve 0");
219: if (trans) KSPSolveTranspose(ksp, b, x);
220: else KSPSolve(ksp, b, x);
222: if (permute) VecPermute(x, rowperm, PETSC_TRUE);
224: CheckResult(&ksp, &A, &b, &x, &rowperm);
226: /*=========================
227: solve a large system
228: =========================*/
230: PetscPreLoadStage("Load System 1");
232: CreateSystem(file[1], rhstype, ordering, permute, &rowperm, &A, &b, &x);
234: PetscPreLoadStage("KSPSetUp 1");
235: KSPCreate(PETSC_COMM_WORLD, &ksp);
236: KSPSetOperators(ksp, A, A);
237: KSPSetFromOptions(ksp);
239: /*
240: Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
241: enable more precise profiling of setting up the preconditioner.
242: These calls are optional, since both will be called within
243: KSPSolve() if they haven't been called already.
244: */
245: KSPSetUp(ksp);
246: KSPSetUpOnBlocks(ksp);
248: PetscPreLoadStage("KSPSolve 1");
249: if (trans) KSPSolveTranspose(ksp, b, x);
250: else KSPSolve(ksp, b, x);
252: if (permute) VecPermute(x, rowperm, PETSC_TRUE);
254: CheckResult(&ksp, &A, &b, &x, &rowperm);
256: PetscPreLoadEnd();
257: /*
258: Always call PetscFinalize() before exiting a program. This routine
259: - finalizes the PETSc libraries as well as MPI
260: - provides summary and diagnostic information if certain runtime
261: options are chosen (e.g., -log_view).
262: */
263: PetscFinalize();
264: return 0;
265: }
267: /*TEST
269: test:
270: TODO: Matrix row/column sizes are not compatible with block size
271: suffix: 1
272: nsize: 4
273: output_file: output/ex10_1.out
274: requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
275: args: -f0 ${DATAFILESPATH}/matrices/medium -f1 ${DATAFILESPATH}/matrices/arco6 -ksp_gmres_classicalgramschmidt -mat_type baij -matload_block_size 3 -pc_type bjacobi
277: test:
278: TODO: Matrix row/column sizes are not compatible with block size
279: suffix: 2
280: nsize: 4
281: output_file: output/ex10_2.out
282: requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
283: args: -f0 ${DATAFILESPATH}/matrices/medium -f1 ${DATAFILESPATH}/matrices/arco6 -ksp_gmres_classicalgramschmidt -mat_type baij -matload_block_size 3 -pc_type bjacobi -trans
285: test:
286: suffix: 3
287: requires: double complex !defined(PETSC_USE_64BIT_INDICES)
288: args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/nh-complex-int32-float64 -ksp_type bicg
290: test:
291: suffix: 4
292: args: -f ${DATAFILESPATH}/matrices/medium -ksp_type bicg -permute rcm
293: requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
295: TEST*/