Actual source code: ex10.c
petsc-3.13.6 2020-09-29
2: static char help[] = "Solve a small system and a large system through preloading\n\
3: Input arguments are:\n\
4: -f0 <small_sys_binary> -f1 <large_sys_binary> \n\n";
6: /*T
7: Concepts: KSP^basic parallel example
8: Concepts: Mat^loading a binary matrix and vector;
9: Concepts: PetscLog^preloading executable
10: Processors: n
11: T*/
13: /*
14: Include "petscksp.h" so that we can use KSP solvers. Note that this file
15: automatically includes:
16: petscsys.h - base PETSc routines petscvec.h - vectors
17: petscmat.h - matrices
18: petscis.h - index sets petscksp.h - Krylov subspace methods
19: petscviewer.h - viewers petscpc.h - preconditioners
20: */
21: #include <petscksp.h>
23: typedef enum {
24: RHS_FILE,
25: RHS_ONE,
26: RHS_RANDOM
27: } RHSType;
28: const char *const RHSTypes[] = {"FILE", "ONE", "RANDOM", "RHSType", "RHS_", NULL};
30: /* ATTENTION: this is the example used in the Profiling chaper of the PETSc manual,
31: where we referenced its profiling stages, preloading and output etc.
32: When you modify it, please make sure it is still consistent with the manual.
33: */
34: int main(int argc,char **args)
35: {
36: PetscErrorCode ierr;
37: Vec x,b,b2;
38: Mat A; /* linear system matrix */
39: KSP ksp; /* Krylov subspace method context */
40: PetscReal norm; /* norm of solution error */
41: char file[2][PETSC_MAX_PATH_LEN];
42: PetscViewer viewer; /* viewer */
43: PetscBool flg,preload=PETSC_FALSE,same,trans=PETSC_FALSE;
44: RHSType rhstype = RHS_FILE;
45: PetscInt its,j,len,start,idx,n1,n2;
46: const PetscScalar *val;
48: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
50: /*
51: Determine files from which we read the two linear systems
52: (matrix and right-hand-side vector).
53: */
54: PetscOptionsGetBool(NULL,NULL,"-trans",&trans,&flg);
55: PetscOptionsGetString(NULL,NULL,"-f",file[0],PETSC_MAX_PATH_LEN,&flg);
56: if (flg) {
57: PetscStrcpy(file[1],file[0]);
58: preload = PETSC_FALSE;
59: } else {
60: PetscOptionsGetString(NULL,NULL,"-f0",file[0],PETSC_MAX_PATH_LEN,&flg);
61: if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Must indicate binary file with the -f0 or -f option");
62: PetscOptionsGetString(NULL,NULL,"-f1",file[1],PETSC_MAX_PATH_LEN,&flg);
63: if (!flg) preload = PETSC_FALSE; /* don't bother with second system */
64: }
65: PetscOptionsGetEnum(NULL,NULL,"-rhs",RHSTypes,(PetscEnum*)&rhstype,NULL);
67: /*
68: To use preloading, one usually has code like the following:
70: PetscPreLoadBegin(preload,"first stage);
71: lines of code
72: PetscPreLoadStage("second stage");
73: lines of code
74: PetscPreLoadEnd();
76: The two macro PetscPreLoadBegin() and PetscPreLoadEnd() implicitly form a
77: loop with maximal two iterations, depending whether preloading is turned on or
78: not. If it is, either through the preload arg of PetscPreLoadBegin or through
79: -preload command line, the trip count is 2, otherwise it is 1. One can use the
80: predefined variable PetscPreLoadIt within the loop body to get the current
81: iteration number, which is 0 or 1. If preload is turned on, the runtime doesn't
82: do profiling for the first iteration, but it will do profiling for the second
83: iteration instead.
85: One can solve a small system in the first iteration and a large system in
86: the second iteration. This process preloads the instructions with the small
87: system so that more accurate performance monitoring (via -log_view) can be done
88: with the large one (that actually is the system of interest).
90: But in this example, we turned off preloading and duplicated the code for
91: the large system. In general, it is a bad practice and one should not duplicate
92: code. We do that because we want to show profiling stages for both the small
93: system and the large system.
94: */
95: PetscPreLoadBegin(preload,"Load System 0");
97: /*=========================
98: solve a small system
99: =========================*/
101: /* open binary file. Note that we use FILE_MODE_READ to indicate reading from this file */
102: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[0],FILE_MODE_READ,&viewer);
104: /* load the matrix and vector; then destroy the viewer */
105: MatCreate(PETSC_COMM_WORLD,&A);
106: MatSetFromOptions(A);
107: MatLoad(A,viewer);
108: switch (rhstype) {
109: case RHS_FILE:
110: /* Vectors in the file might a different size than the matrix so we need a
111: * Vec whose size hasn't been set yet. It'll get fixed below. Otherwise we
112: * can create the correct size Vec. */
113: VecCreate(PETSC_COMM_WORLD,&b);
114: VecLoad(b,viewer);
115: break;
116: case RHS_ONE:
117: MatCreateVecs(A,&b,NULL);
118: VecSet(b,1.0);
119: break;
120: case RHS_RANDOM:
121: MatCreateVecs(A,&b,NULL);
122: VecSetRandom(b,NULL);
123: break;
124: }
125: PetscViewerDestroy(&viewer);
127: /* if the loaded matrix is larger than the vector (due to being padded
128: to match the block size of the system), then create a new padded vector
129: */
130: MatGetLocalSize(A,NULL,&n1);
131: VecGetLocalSize(b,&n2);
132: same = (n1 == n2)? PETSC_TRUE : PETSC_FALSE;
133: MPIU_Allreduce(MPI_IN_PLACE,&same,1,MPIU_BOOL,MPI_LAND,PETSC_COMM_WORLD);
135: if (!same) { /* create a new vector b by padding the old one */
136: VecCreate(PETSC_COMM_WORLD,&b2);
137: VecSetSizes(b2,n1,PETSC_DECIDE);
138: VecSetFromOptions(b2);
139: VecGetOwnershipRange(b,&start,NULL);
140: VecGetLocalSize(b,&len);
141: VecGetArrayRead(b,&val);
142: for (j=0; j<len; j++) {
143: idx = start+j;
144: VecSetValues(b2,1,&idx,val+j,INSERT_VALUES);
145: }
146: VecRestoreArrayRead(b,&val);
147: VecDestroy(&b);
148: VecAssemblyBegin(b2);
149: VecAssemblyEnd(b2);
150: b = b2;
151: }
152: VecDuplicate(b,&x);
154: PetscPreLoadStage("KSPSetUp 0");
155: KSPCreate(PETSC_COMM_WORLD,&ksp);
156: KSPSetOperators(ksp,A,A);
157: KSPSetFromOptions(ksp);
159: /*
160: Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
161: enable more precise profiling of setting up the preconditioner.
162: These calls are optional, since both will be called within
163: KSPSolve() if they haven't been called already.
164: */
165: KSPSetUp(ksp);
166: KSPSetUpOnBlocks(ksp);
168: PetscPreLoadStage("KSPSolve 0");
169: if (trans) {KSPSolveTranspose(ksp,b,x);}
170: else {KSPSolve(ksp,b,x);}
172: KSPGetTotalIterations(ksp,&its);
173: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %d\n",its);
175: KSPGetResidualNorm(ksp,&norm);
176: if (norm < 1.e-12) {
177: PetscPrintf(PETSC_COMM_WORLD,"Residual norm < 1.e-12\n");
178: } else {
179: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %e\n",(double)norm);
180: }
182: KSPDestroy(&ksp);
183: MatDestroy(&A);
184: VecDestroy(&x);
185: VecDestroy(&b);
187: /*=========================
188: solve a large system
189: =========================*/
190: /* the code is duplicated. Bad practice. See comments above */
191: PetscPreLoadStage("Load System 1");
192: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[1],FILE_MODE_READ,&viewer);
194: /* load the matrix and vector; then destroy the viewer */
195: MatCreate(PETSC_COMM_WORLD,&A);
196: MatSetFromOptions(A);
197: MatLoad(A,viewer);
198: switch (rhstype) {
199: case RHS_FILE:
200: /* Vectors in the file might a different size than the matrix so we need a
201: * Vec whose size hasn't been set yet. It'll get fixed below. Otherwise we
202: * can create the correct size Vec. */
203: VecCreate(PETSC_COMM_WORLD,&b);
204: VecLoad(b,viewer);
205: break;
206: case RHS_ONE:
207: MatCreateVecs(A,&b,NULL);
208: VecSet(b,1.0);
209: break;
210: case RHS_RANDOM:
211: MatCreateVecs(A,&b,NULL);
212: VecSetRandom(b,NULL);
213: break;
214: }
215: PetscViewerDestroy(&viewer);
217: MatGetLocalSize(A,NULL,&n1);
218: VecGetLocalSize(b,&n2);
219: same = (n1 == n2)? PETSC_TRUE : PETSC_FALSE;
220: MPIU_Allreduce(MPI_IN_PLACE,&same,1,MPIU_BOOL,MPI_LAND,PETSC_COMM_WORLD);
222: if (!same) { /* create a new vector b by padding the old one */
223: VecCreate(PETSC_COMM_WORLD,&b2);
224: VecSetSizes(b2,n1,PETSC_DECIDE);
225: VecSetFromOptions(b2);
226: VecGetOwnershipRange(b,&start,NULL);
227: VecGetLocalSize(b,&len);
228: VecGetArrayRead(b,&val);
229: for (j=0; j<len; j++) {
230: idx = start+j;
231: VecSetValues(b2,1,&idx,val+j,INSERT_VALUES);
232: }
233: VecRestoreArrayRead(b,&val);
234: VecDestroy(&b);
235: VecAssemblyBegin(b2);
236: VecAssemblyEnd(b2);
237: b = b2;
238: }
239: VecDuplicate(b,&x);
241: PetscPreLoadStage("KSPSetUp 1");
242: KSPCreate(PETSC_COMM_WORLD,&ksp);
243: KSPSetOperators(ksp,A,A);
244: KSPSetFromOptions(ksp);
245: /*
246: Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
247: enable more precise profiling of setting up the preconditioner.
248: These calls are optional, since both will be called within
249: KSPSolve() if they haven't been called already.
250: */
251: KSPSetUp(ksp);
252: KSPSetUpOnBlocks(ksp);
254: PetscPreLoadStage("KSPSolve 1");
255: if (trans) {KSPSolveTranspose(ksp,b,x);}
256: else {KSPSolve(ksp,b,x);}
258: KSPGetTotalIterations(ksp,&its);
259: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %d\n",its);
261: KSPGetResidualNorm(ksp,&norm);
262: if (norm < 1.e-12) {
263: PetscPrintf(PETSC_COMM_WORLD,"Residual norm < 1.e-12\n");
264: } else {
265: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %e\n",(double)norm);
266: }
268: KSPDestroy(&ksp);
269: MatDestroy(&A);
270: VecDestroy(&x);
271: VecDestroy(&b);
272: PetscPreLoadEnd();
273: /*
274: Always call PetscFinalize() before exiting a program. This routine
275: - finalizes the PETSc libraries as well as MPI
276: - provides summary and diagnostic information if certain runtime
277: options are chosen (e.g., -log_view).
278: */
279: PetscFinalize();
280: return ierr;
281: }
283: /*TEST
285: test:
286: TODO: Matrix row/column sizes are not compatible with block size
287: suffix: 1
288: nsize: 4
289: output_file: output/ex10_1.out
290: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
291: args: -f0 ${DATAFILESPATH}/matrices/medium -f1 ${DATAFILESPATH}/matrices/arco6 -ksp_gmres_classicalgramschmidt -mat_type baij -matload_block_size 3 -pc_type bjacobi
293: test:
294: TODO: Matrix row/column sizes are not compatible with block size
295: suffix: 2
296: nsize: 4
297: output_file: output/ex10_2.out
298: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
299: args: -f0 ${DATAFILESPATH}/matrices/medium -f1 ${DATAFILESPATH}/matrices/arco6 -ksp_gmres_classicalgramschmidt -mat_type baij -matload_block_size 3 -pc_type bjacobi -trans
301: TEST*/