Actual source code: ex10.c

petsc-3.11.3 2019-06-26
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: /* 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*/