Actual source code: ex10.c

petsc-3.9.2 2018-05-20
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;
 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);

217:   KSPGetTotalIterations(ksp,&its);
218:   PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %d\n",its);

220:   KSPGetResidualNorm(ksp,&norm);
221:   if (norm < 1.e-12) {
222:     PetscPrintf(PETSC_COMM_WORLD,"Residual norm < 1.e-12\n");
223:   } else {
224:     PetscPrintf(PETSC_COMM_WORLD,"Residual norm %e\n",(double)norm);
225:   }

227:   KSPDestroy(&ksp);
228:   MatDestroy(&A);
229:   VecDestroy(&x);
230:   VecDestroy(&b);

232:   PetscPreLoadEnd();
233:   /*
234:      Always call PetscFinalize() before exiting a program.  This routine
235:        - finalizes the PETSc libraries as well as MPI
236:        - provides summary and diagnostic information if certain runtime
237:          options are chosen (e.g., -log_view).
238:   */
239:   PetscFinalize();
240:   return ierr;
241: }

243: /*TEST

245:    test:
246:       nsize: 4
247:       output_file: output/ex10_1.out
248:       requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
249:       args: -f0 ${DATAFILESPATH}/matrices/medium -f1 ${DATAFILESPATH}/matrices/arco6 -ksp_gmres_classicalgramschmidt -mat_type baij -matload_block_size 3 -pc_type bjacobi
250: TEST*/