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