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