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