Actual source code: ex10.c
petsc-master 2019-04-20
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: /* ATTENTION: this is the example used in the Profiling chaper of the PETSc manual,
24: where we referenced its profiling stages, preloading and output etc.
25: When you modify it, please make sure it is still consistent with the manual.
26: */
27: int main(int argc,char **args)
28: {
29: PetscErrorCode ierr;
30: Vec x,b,b2;
31: Mat A; /* linear system matrix */
32: KSP ksp; /* Krylov subspace method context */
33: PetscReal norm; /* norm of solution error */
34: char file[2][PETSC_MAX_PATH_LEN];
35: PetscViewer viewer; /* viewer */
36: PetscBool flg,preload=PETSC_FALSE,same,trans=PETSC_FALSE;
37: PetscInt its,j,len,start,idx,n1,n2;
38: const PetscScalar *val;
40: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
42: /*
43: Determine files from which we read the two linear systems
44: (matrix and right-hand-side vector).
45: */
46: PetscOptionsGetBool(NULL,NULL,"-trans",&trans,&flg);
47: PetscOptionsGetString(NULL,NULL,"-f",file[0],PETSC_MAX_PATH_LEN,&flg);
48: if (flg) {
49: PetscStrcpy(file[1],file[0]);
50: preload = PETSC_FALSE;
51: } else {
52: PetscOptionsGetString(NULL,NULL,"-f0",file[0],PETSC_MAX_PATH_LEN,&flg);
53: if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate binary file with the -f0 or -f option");
54: PetscOptionsGetString(NULL,NULL,"-f1",file[1],PETSC_MAX_PATH_LEN,&flg);
55: if (!flg) preload = PETSC_FALSE; /* don't bother with second system */
56: }
58: /*
59: To use preloading, one usually has code like the following:
61: PetscPreLoadBegin(preload,"first stage);
62: lines of code
63: PetscPreLoadStage("second stage");
64: lines of code
65: PetscPreLoadEnd();
67: The two macro PetscPreLoadBegin() and PetscPreLoadEnd() implicitly form a
68: loop with maximal two iterations, depending whether preloading is turned on or
69: not. If it is, either through the preload arg of PetscPreLoadBegin or through
70: -preload command line, the trip count is 2, otherwise it is 1. One can use the
71: predefined variable PetscPreLoadIt within the loop body to get the current
72: iteration number, which is 0 or 1. If preload is turned on, the runtime doesn't
73: do profiling for the first iteration, but it will do profiling for the second
74: iteration instead.
76: One can solve a small system in the first iteration and a large system in
77: the second iteration. This process preloads the instructions with the small
78: system so that more accurate performance monitoring (via -log_view) can be done
79: with the large one (that actually is the system of interest).
81: But in this example, we turned off preloading and duplicated the code for
82: the large system. In general, it is a bad practice and one should not duplicate
83: code. We do that because we want to show profiling stages for both the small
84: system and the large system.
85: */
86: PetscPreLoadBegin(preload,"Load System 0");
88: /*=========================
89: solve a small system
90: =========================*/
92: /* open binary file. Note that we use FILE_MODE_READ to indicate reading from this file */
93: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[0],FILE_MODE_READ,&viewer);
95: /* load the matrix and vector; then destroy the viewer */
96: MatCreate(PETSC_COMM_WORLD,&A);
97: VecCreate(PETSC_COMM_WORLD,&b);
98: MatSetFromOptions(A);
99: MatLoad(A,viewer);
100: VecLoad(b,viewer);
101: PetscViewerDestroy(&viewer);
103: /* if the loaded matrix is larger than the vector (due to being padded
104: to match the block size of the system), then create a new padded vector
105: */
106: MatGetLocalSize(A,NULL,&n1);
107: VecGetLocalSize(b,&n2);
108: same = (n1 == n2)? PETSC_TRUE : PETSC_FALSE;
109: MPIU_Allreduce(MPI_IN_PLACE,&same,1,MPIU_BOOL,MPI_LAND,PETSC_COMM_WORLD);
111: if (!same) { /* create a new vector b by padding the old one */
112: VecCreate(PETSC_COMM_WORLD,&b2);
113: VecSetSizes(b2,n1,PETSC_DECIDE);
114: VecSetFromOptions(b2);
115: VecGetOwnershipRange(b,&start,NULL);
116: VecGetLocalSize(b,&len);
117: VecGetArrayRead(b,&val);
118: for (j=0; j<len; j++) {
119: idx = start+j;
120: VecSetValues(b2,1,&idx,val+j,INSERT_VALUES);
121: }
122: VecRestoreArrayRead(b,&val);
123: VecDestroy(&b);
124: VecAssemblyBegin(b2);
125: VecAssemblyEnd(b2);
126: b = b2;
127: }
128: VecDuplicate(b,&x);
130: PetscPreLoadStage("KSPSetUp 0");
131: KSPCreate(PETSC_COMM_WORLD,&ksp);
132: KSPSetOperators(ksp,A,A);
133: KSPSetFromOptions(ksp);
135: /*
136: Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
137: enable more precise profiling of setting up the preconditioner.
138: These calls are optional, since both will be called within
139: KSPSolve() if they haven't been called already.
140: */
141: KSPSetUp(ksp);
142: KSPSetUpOnBlocks(ksp);
144: PetscPreLoadStage("KSPSolve 0");
145: if (trans) {KSPSolveTranspose(ksp,b,x);}
146: else {KSPSolve(ksp,b,x);}
148: KSPGetTotalIterations(ksp,&its);
149: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %d\n",its);
151: KSPGetResidualNorm(ksp,&norm);
152: if (norm < 1.e-12) {
153: PetscPrintf(PETSC_COMM_WORLD,"Residual norm < 1.e-12\n");
154: } else {
155: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %e\n",(double)norm);
156: }
158: KSPDestroy(&ksp);
159: MatDestroy(&A);
160: VecDestroy(&x);
161: VecDestroy(&b);
163: /*=========================
164: solve a large system
165: =========================*/
166: /* the code is duplicated. Bad practice. See comments above */
167: PetscPreLoadStage("Load System 1");
168: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[1],FILE_MODE_READ,&viewer);
170: /* load the matrix and vector; then destroy the viewer */
171: MatCreate(PETSC_COMM_WORLD,&A);
172: VecCreate(PETSC_COMM_WORLD,&b);
173: MatSetFromOptions(A);
174: MatLoad(A,viewer);
175: VecLoad(b,viewer);
176: PetscViewerDestroy(&viewer);
178: MatGetLocalSize(A,NULL,&n1);
179: VecGetLocalSize(b,&n2);
180: same = (n1 == n2)? PETSC_TRUE : PETSC_FALSE;
181: MPIU_Allreduce(MPI_IN_PLACE,&same,1,MPIU_BOOL,MPI_LAND,PETSC_COMM_WORLD);
183: if (!same) { /* create a new vector b by padding the old one */
184: VecCreate(PETSC_COMM_WORLD,&b2);
185: VecSetSizes(b2,n1,PETSC_DECIDE);
186: VecSetFromOptions(b2);
187: VecGetOwnershipRange(b,&start,NULL);
188: VecGetLocalSize(b,&len);
189: VecGetArrayRead(b,&val);
190: for (j=0; j<len; j++) {
191: idx = start+j;
192: VecSetValues(b2,1,&idx,val+j,INSERT_VALUES);
193: }
194: VecRestoreArrayRead(b,&val);
195: VecDestroy(&b);
196: VecAssemblyBegin(b2);
197: VecAssemblyEnd(b2);
198: b = b2;
199: }
200: VecDuplicate(b,&x);
202: PetscPreLoadStage("KSPSetUp 1");
203: KSPCreate(PETSC_COMM_WORLD,&ksp);
204: KSPSetOperators(ksp,A,A);
205: KSPSetFromOptions(ksp);
206: /*
207: Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
208: enable more precise profiling of setting up the preconditioner.
209: These calls are optional, since both will be called within
210: KSPSolve() if they haven't been called already.
211: */
212: KSPSetUp(ksp);
213: KSPSetUpOnBlocks(ksp);
215: PetscPreLoadStage("KSPSolve 1");
216: if (trans) {KSPSolveTranspose(ksp,b,x);}
217: else {KSPSolve(ksp,b,x);}
219: KSPGetTotalIterations(ksp,&its);
220: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %d\n",its);
222: KSPGetResidualNorm(ksp,&norm);
223: if (norm < 1.e-12) {
224: PetscPrintf(PETSC_COMM_WORLD,"Residual norm < 1.e-12\n");
225: } else {
226: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %e\n",(double)norm);
227: }
229: KSPDestroy(&ksp);
230: MatDestroy(&A);
231: VecDestroy(&x);
232: VecDestroy(&b);
233: PetscPreLoadEnd();
234: /*
235: Always call PetscFinalize() before exiting a program. This routine
236: - finalizes the PETSc libraries as well as MPI
237: - provides summary and diagnostic information if certain runtime
238: options are chosen (e.g., -log_view).
239: */
240: PetscFinalize();
241: return ierr;
242: }
244: /*TEST
246: test:
247: nsize: 4
248: output_file: output/ex10_1.out
249: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
250: args: -f0 ${DATAFILESPATH}/matrices/medium -f1 ${DATAFILESPATH}/matrices/arco6 -ksp_gmres_classicalgramschmidt -mat_type baij -matload_block_size 3 -pc_type bjacobi
252: test:
253: suffix: 2
254: nsize: 4
255: output_file: output/ex10_2.out
256: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
257: args: -f0 ${DATAFILESPATH}/matrices/medium -f1 ${DATAFILESPATH}/matrices/arco6 -ksp_gmres_classicalgramschmidt -mat_type baij -matload_block_size 3 -pc_type bjacobi -trans
259: TEST*/