Actual source code: ex10.c

petsc-3.12.2 2019-11-22
Report Typos and Errors

  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,1,"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:       nsize: 4
287:       output_file: output/ex10_1.out
288:       requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
289:       args: -f0 ${DATAFILESPATH}/matrices/medium -f1 ${DATAFILESPATH}/matrices/arco6 -ksp_gmres_classicalgramschmidt -mat_type baij -matload_block_size 3 -pc_type bjacobi

291:    test:
292:       suffix: 2
293:       nsize: 4
294:       output_file: output/ex10_2.out
295:       requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
296:       args: -f0 ${DATAFILESPATH}/matrices/medium -f1 ${DATAFILESPATH}/matrices/arco6 -ksp_gmres_classicalgramschmidt -mat_type baij -matload_block_size 3 -pc_type bjacobi -trans

298: TEST*/