Actual source code: ex10.c
petsc-3.10.5 2019-03-28
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;
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: PetscOptionsGetString(NULL,NULL,"-f",file[0],PETSC_MAX_PATH_LEN,&flg);
47: if (flg) {
48: PetscStrcpy(file[1],file[0]);
49: preload = PETSC_FALSE;
50: } else {
51: PetscOptionsGetString(NULL,NULL,"-f0",file[0],PETSC_MAX_PATH_LEN,&flg);
52: if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate binary file with the -f0 or -f option");
53: PetscOptionsGetString(NULL,NULL,"-f1",file[1],PETSC_MAX_PATH_LEN,&flg);
54: if (!flg) preload = PETSC_FALSE; /* don't bother with second system */
55: }
57: /*
58: To use preloading, one usually has code like the following:
60: PetscPreLoadBegin(preload,"first stage);
61: lines of code
62: PetscPreLoadStage("second stage");
63: lines of code
64: PetscPreLoadEnd();
66: The two macro PetscPreLoadBegin() and PetscPreLoadEnd() implicitly form a
67: loop with maximal two iterations, depending whether preloading is turned on or
68: not. If it is, either through the preload arg of PetscPreLoadBegin or through
69: -preload command line, the trip count is 2, otherwise it is 1. One can use the
70: predefined variable PetscPreLoadIt within the loop body to get the current
71: iteration number, which is 0 or 1. If preload is turned on, the runtime doesn't
72: do profiling for the first iteration, but it will do profiling for the second
73: iteration instead.
75: One can solve a small system in the first iteration and a large system in
76: the second iteration. This process preloads the instructions with the small
77: system so that more accurate performance monitoring (via -log_view) can be done
78: with the large one (that actually is the system of interest).
80: But in this example, we turned off preloading and duplicated the code for
81: the large system. In general, it is a bad practice and one should not duplicate
82: code. We do that because we want to show profiling stages for both the small
83: system and the large system.
84: */
85: PetscPreLoadBegin(preload,"Load System 0");
87: /*=========================
88: solve a samll system
89: =========================*/
91: /* open binary file. Note that we use FILE_MODE_READ to indicate reading from this file */
92: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[0],FILE_MODE_READ,&viewer);
94: /* load the matrix and vector; then destroy the viewer */
95: MatCreate(PETSC_COMM_WORLD,&A);
96: VecCreate(PETSC_COMM_WORLD,&b);
97: MatSetFromOptions(A);
98: MatLoad(A,viewer);
99: VecLoad(b,viewer);
100: PetscViewerDestroy(&viewer);
102: /* if the loaded matrix is larger than the vector (due to being padded
103: to match the block size of the system), then create a new padded vector
104: */
105: MatGetLocalSize(A,NULL,&n1);
106: VecGetLocalSize(b,&n2);
107: same = (n1 == n2)? PETSC_TRUE : PETSC_FALSE;
108: MPIU_Allreduce(MPI_IN_PLACE,&same,1,MPIU_BOOL,MPI_LAND,PETSC_COMM_WORLD);
110: if (!same) { /* create a new vector b by padding the old one */
111: VecCreate(PETSC_COMM_WORLD,&b2);
112: VecSetSizes(b2,n1,PETSC_DECIDE);
113: VecSetFromOptions(b2);
114: VecGetOwnershipRange(b,&start,NULL);
115: VecGetLocalSize(b,&len);
116: VecGetArrayRead(b,&val);
117: for (j=0; j<len; j++) {
118: idx = start+j;
119: VecSetValues(b2,1,&idx,val+j,INSERT_VALUES);
120: }
121: VecRestoreArrayRead(b,&val);
122: VecDestroy(&b);
123: VecAssemblyBegin(b2);
124: VecAssemblyEnd(b2);
125: b = b2;
126: }
127: VecDuplicate(b,&x);
129: PetscPreLoadStage("KSPSetUp 0");
130: KSPCreate(PETSC_COMM_WORLD,&ksp);
131: KSPSetOperators(ksp,A,A);
132: KSPSetFromOptions(ksp);
134: /*
135: Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
136: enable more precise profiling of setting up the preconditioner.
137: These calls are optional, since both will be called within
138: KSPSolve() if they haven't been called already.
139: */
140: KSPSetUp(ksp);
141: KSPSetUpOnBlocks(ksp);
143: PetscPreLoadStage("KSPSolve 0");
144: KSPSolve(ksp,b,x);
146: KSPGetTotalIterations(ksp,&its);
147: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %d\n",its);
149: KSPGetResidualNorm(ksp,&norm);
150: if (norm < 1.e-12) {
151: PetscPrintf(PETSC_COMM_WORLD,"Residual norm < 1.e-12\n");
152: } else {
153: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %e\n",(double)norm);
154: }
156: KSPDestroy(&ksp);
157: MatDestroy(&A);
158: VecDestroy(&x);
159: VecDestroy(&b);
161: /*=========================
162: solve a large system
163: =========================*/
165: /* the code is duplicated. Bad practice. See comments above */
166: PetscPreLoadStage("Load System 1");
167: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[1],FILE_MODE_READ,&viewer);
169: /* load the matrix and vector; then destroy the viewer */
170: MatCreate(PETSC_COMM_WORLD,&A);
171: VecCreate(PETSC_COMM_WORLD,&b);
172: MatSetFromOptions(A);
173: MatLoad(A,viewer);
174: VecLoad(b,viewer);
175: PetscViewerDestroy(&viewer);
177: MatGetLocalSize(A,NULL,&n1);
178: VecGetLocalSize(b,&n2);
179: same = (n1 == n2)? PETSC_TRUE : PETSC_FALSE;
180: MPIU_Allreduce(MPI_IN_PLACE,&same,1,MPIU_BOOL,MPI_LAND,PETSC_COMM_WORLD);
182: if (!same) { /* create a new vector b by padding the old one */
183: VecCreate(PETSC_COMM_WORLD,&b2);
184: VecSetSizes(b2,n1,PETSC_DECIDE);
185: VecSetFromOptions(b2);
186: VecGetOwnershipRange(b,&start,NULL);
187: VecGetLocalSize(b,&len);
188: VecGetArrayRead(b,&val);
189: for (j=0; j<len; j++) {
190: idx = start+j;
191: VecSetValues(b2,1,&idx,val+j,INSERT_VALUES);
192: }
193: VecRestoreArrayRead(b,&val);
194: VecDestroy(&b);
195: VecAssemblyBegin(b2);
196: VecAssemblyEnd(b2);
197: b = b2;
198: }
199: VecDuplicate(b,&x);
201: PetscPreLoadStage("KSPSetUp 1");
202: KSPCreate(PETSC_COMM_WORLD,&ksp);
203: KSPSetOperators(ksp,A,A);
204: KSPSetFromOptions(ksp);
205: /*
206: Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
207: enable more precise profiling of setting up the preconditioner.
208: These calls are optional, since both will be called within
209: KSPSolve() if they haven't been called already.
210: */
211: KSPSetUp(ksp);
212: KSPSetUpOnBlocks(ksp);
214: PetscPreLoadStage("KSPSolve 1");
215: KSPSolve(ksp,b,x);
216: KSPGetTotalIterations(ksp,&its);
217: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %d\n",its);
219: KSPGetResidualNorm(ksp,&norm);
220: if (norm < 1.e-12) {
221: PetscPrintf(PETSC_COMM_WORLD,"Residual norm < 1.e-12\n");
222: } else {
223: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %e\n",(double)norm);
224: }
226: KSPDestroy(&ksp);
227: MatDestroy(&A);
228: VecDestroy(&x);
229: VecDestroy(&b);
231: PetscPreLoadEnd();
232: /*
233: Always call PetscFinalize() before exiting a program. This routine
234: - finalizes the PETSc libraries as well as MPI
235: - provides summary and diagnostic information if certain runtime
236: options are chosen (e.g., -log_view).
237: */
238: PetscFinalize();
239: return ierr;
240: }
242: /*TEST
244: test:
245: nsize: 4
246: output_file: output/ex10_1.out
247: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
248: args: -f0 ${DATAFILESPATH}/matrices/medium -f1 ${DATAFILESPATH}/matrices/arco6 -ksp_gmres_classicalgramschmidt -mat_type baij -matload_block_size 3 -pc_type bjacobi
250: TEST*/