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