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