Actual source code: ex10.c


  2: static char help[] = "Solve a small system and a large system through preloading\n\
  3:   Input arguments are:\n\
  4:   -permute <natural,rcm,nd,...> : solve system in permuted indexing\n\
  5:   -f0 <small_sys_binary> -f1 <large_sys_binary> \n\n";

  7: /*
  8:   Include "petscksp.h" so that we can use KSP solvers.  Note that this file
  9:   automatically includes:
 10:      petscsys.h       - base PETSc routines   petscvec.h - vectors
 11:      petscmat.h - matrices
 12:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 13:      petscviewer.h - viewers               petscpc.h  - preconditioners
 14: */
 15: #include <petscksp.h>

 17: typedef enum {
 18:   RHS_FILE,
 19:   RHS_ONE,
 20:   RHS_RANDOM
 21: } RHSType;
 22: const char *const RHSTypes[] = {"FILE", "ONE", "RANDOM", "RHSType", "RHS_", NULL};

 24: PetscErrorCode CheckResult(KSP *ksp, Mat *A, Vec *b, Vec *x, IS *rowperm)
 25: {
 26:   PetscReal norm; /* norm of solution error */
 27:   PetscInt  its;
 28:   KSPGetTotalIterations(*ksp, &its);
 29:   PetscPrintf(PETSC_COMM_WORLD, "Number of iterations = %" PetscInt_FMT "\n", its);

 31:   KSPGetResidualNorm(*ksp, &norm);
 32:   if (norm < 1.e-12) {
 33:     PetscPrintf(PETSC_COMM_WORLD, "Residual norm < 1.e-12\n");
 34:   } else {
 35:     PetscPrintf(PETSC_COMM_WORLD, "Residual norm %e\n", (double)norm);
 36:   }

 38:   KSPDestroy(ksp);
 39:   MatDestroy(A);
 40:   VecDestroy(x);
 41:   VecDestroy(b);
 42:   ISDestroy(rowperm);
 43:   return 0;
 44: }

 46: PetscErrorCode CreateSystem(const char filename[PETSC_MAX_PATH_LEN], RHSType rhstype, MatOrderingType ordering, PetscBool permute, IS *rowperm_out, Mat *A_out, Vec *b_out, Vec *x_out)
 47: {
 48:   Vec                x, b, b2;
 49:   Mat                A;      /* linear system matrix */
 50:   PetscViewer        viewer; /* viewer */
 51:   PetscBool          same;
 52:   PetscInt           j, len, start, idx, n1, n2;
 53:   const PetscScalar *val;
 54:   IS                 rowperm = NULL, colperm = NULL;

 56:   /* open binary file. Note that we use FILE_MODE_READ to indicate reading from this file */
 57:   PetscViewerBinaryOpen(PETSC_COMM_WORLD, filename, FILE_MODE_READ, &viewer);

 59:   /* load the matrix and vector; then destroy the viewer */
 60:   MatCreate(PETSC_COMM_WORLD, &A);
 61:   MatSetFromOptions(A);
 62:   MatLoad(A, viewer);
 63:   switch (rhstype) {
 64:   case RHS_FILE:
 65:     /* Vectors in the file might a different size than the matrix so we need a
 66:      * Vec whose size hasn't been set yet.  It'll get fixed below.  Otherwise we
 67:      * can create the correct size Vec. */
 68:     VecCreate(PETSC_COMM_WORLD, &b);
 69:     VecLoad(b, viewer);
 70:     break;
 71:   case RHS_ONE:
 72:     MatCreateVecs(A, &b, NULL);
 73:     VecSet(b, 1.0);
 74:     break;
 75:   case RHS_RANDOM:
 76:     MatCreateVecs(A, &b, NULL);
 77:     VecSetRandom(b, NULL);
 78:     break;
 79:   }
 80:   PetscViewerDestroy(&viewer);

 82:   /* if the loaded matrix is larger than the vector (due to being padded
 83:      to match the block size of the system), then create a new padded vector
 84:    */
 85:   MatGetLocalSize(A, NULL, &n1);
 86:   VecGetLocalSize(b, &n2);
 87:   same = (n1 == n2) ? PETSC_TRUE : PETSC_FALSE;
 88:   MPIU_Allreduce(MPI_IN_PLACE, &same, 1, MPIU_BOOL, MPI_LAND, PETSC_COMM_WORLD);

 90:   if (!same) { /* create a new vector b by padding the old one */
 91:     VecCreate(PETSC_COMM_WORLD, &b2);
 92:     VecSetSizes(b2, n1, PETSC_DECIDE);
 93:     VecSetFromOptions(b2);
 94:     VecGetOwnershipRange(b, &start, NULL);
 95:     VecGetLocalSize(b, &len);
 96:     VecGetArrayRead(b, &val);
 97:     for (j = 0; j < len; j++) {
 98:       idx = start + j;
 99:       VecSetValues(b2, 1, &idx, val + j, INSERT_VALUES);
100:     }
101:     VecRestoreArrayRead(b, &val);
102:     VecDestroy(&b);
103:     VecAssemblyBegin(b2);
104:     VecAssemblyEnd(b2);
105:     b = b2;
106:   }
107:   VecDuplicate(b, &x);

109:   if (permute) {
110:     Mat Aperm;
111:     MatGetOrdering(A, ordering, &rowperm, &colperm);
112:     MatPermute(A, rowperm, colperm, &Aperm);
113:     VecPermute(b, colperm, PETSC_FALSE);
114:     MatDestroy(&A);
115:     A = Aperm; /* Replace original operator with permuted version */
116:     ISDestroy(&colperm);
117:   }

119:   *b_out       = b;
120:   *x_out       = x;
121:   *A_out       = A;
122:   *rowperm_out = rowperm;

124:   return 0;
125: }

127: /* ATTENTION: this is the example used in the Profiling chaper of the PETSc manual,
128:    where we referenced its profiling stages, preloading and output etc.
129:    When you modify it, please make sure it is still consistent with the manual.
130:  */
131: int main(int argc, char **args)
132: {
133:   Vec       x, b;
134:   Mat       A;   /* linear system matrix */
135:   KSP       ksp; /* Krylov subspace method context */
136:   char      file[2][PETSC_MAX_PATH_LEN], ordering[256] = MATORDERINGRCM;
137:   RHSType   rhstype = RHS_FILE;
138:   PetscBool flg, preload = PETSC_FALSE, trans = PETSC_FALSE, permute = PETSC_FALSE;
139:   IS        rowperm = NULL;

142:   PetscInitialize(&argc, &args, (char *)0, help);

144:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Preloading example options", "");
145:   {
146:     /*
147:        Determine files from which we read the two linear systems
148:        (matrix and right-hand-side vector).
149:     */
150:     PetscOptionsBool("-trans", "Solve transpose system instead", "", trans, &trans, &flg);
151:     PetscOptionsString("-f", "First file to load (small system)", "", file[0], file[0], sizeof(file[0]), &flg);
152:     PetscOptionsFList("-permute", "Permute matrix and vector to solve in new ordering", "", MatOrderingList, ordering, ordering, sizeof(ordering), &permute);

154:     if (flg) {
155:       PetscStrcpy(file[1], file[0]);
156:       preload = PETSC_FALSE;
157:     } else {
158:       PetscOptionsString("-f0", "First file to load (small system)", "", file[0], file[0], sizeof(file[0]), &flg);
160:       PetscOptionsString("-f1", "Second file to load (larger system)", "", file[1], file[1], sizeof(file[1]), &flg);
161:       if (!flg) preload = PETSC_FALSE; /* don't bother with second system */
162:     }

164:     PetscOptionsEnum("-rhs", "Right hand side", "", RHSTypes, (PetscEnum)rhstype, (PetscEnum *)&rhstype, NULL);
165:   }
166:   PetscOptionsEnd();

168:   /*
169:     To use preloading, one usually has code like the following:

171:     PetscPreLoadBegin(preload,"first stage);
172:       lines of code
173:     PetscPreLoadStage("second stage");
174:       lines of code
175:     PetscPreLoadEnd();

177:     The two macro PetscPreLoadBegin() and PetscPreLoadEnd() implicitly form a
178:     loop with maximal two iterations, depending whether preloading is turned on or
179:     not. If it is, either through the preload arg of PetscPreLoadBegin or through
180:     -preload command line, the trip count is 2, otherwise it is 1. One can use the
181:     predefined variable PetscPreLoadIt within the loop body to get the current
182:     iteration number, which is 0 or 1. If preload is turned on, the runtime doesn't
183:     do profiling for the first iteration, but it will do profiling for the second
184:     iteration instead.

186:     One can solve a small system in the first iteration and a large system in
187:     the second iteration. This process preloads the instructions with the small
188:     system so that more accurate performance monitoring (via -log_view) can be done
189:     with the large one (that actually is the system of interest).

191:     But in this example, we turned off preloading and duplicated the code for
192:     the large system. In general, it is a bad practice and one should not duplicate
193:     code. We do that because we want to show profiling stages for both the small
194:     system and the large system.
195:   */

197:   /*=========================
198:       solve a small system
199:     =========================*/

201:   PetscPreLoadBegin(preload, "Load System 0");
202:   CreateSystem(file[0], rhstype, ordering, permute, &rowperm, &A, &b, &x);

204:   PetscPreLoadStage("KSPSetUp 0");
205:   KSPCreate(PETSC_COMM_WORLD, &ksp);
206:   KSPSetOperators(ksp, A, A);
207:   KSPSetFromOptions(ksp);

209:   /*
210:     Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
211:     enable more precise profiling of setting up the preconditioner.
212:     These calls are optional, since both will be called within
213:     KSPSolve() if they haven't been called already.
214:   */
215:   KSPSetUp(ksp);
216:   KSPSetUpOnBlocks(ksp);

218:   PetscPreLoadStage("KSPSolve 0");
219:   if (trans) KSPSolveTranspose(ksp, b, x);
220:   else KSPSolve(ksp, b, x);

222:   if (permute) VecPermute(x, rowperm, PETSC_TRUE);

224:   CheckResult(&ksp, &A, &b, &x, &rowperm);

226:   /*=========================
227:     solve a large system
228:     =========================*/

230:   PetscPreLoadStage("Load System 1");

232:   CreateSystem(file[1], rhstype, ordering, permute, &rowperm, &A, &b, &x);

234:   PetscPreLoadStage("KSPSetUp 1");
235:   KSPCreate(PETSC_COMM_WORLD, &ksp);
236:   KSPSetOperators(ksp, A, A);
237:   KSPSetFromOptions(ksp);

239:   /*
240:     Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
241:     enable more precise profiling of setting up the preconditioner.
242:     These calls are optional, since both will be called within
243:     KSPSolve() if they haven't been called already.
244:   */
245:   KSPSetUp(ksp);
246:   KSPSetUpOnBlocks(ksp);

248:   PetscPreLoadStage("KSPSolve 1");
249:   if (trans) KSPSolveTranspose(ksp, b, x);
250:   else KSPSolve(ksp, b, x);

252:   if (permute) VecPermute(x, rowperm, PETSC_TRUE);

254:   CheckResult(&ksp, &A, &b, &x, &rowperm);

256:   PetscPreLoadEnd();
257:   /*
258:      Always call PetscFinalize() before exiting a program.  This routine
259:        - finalizes the PETSc libraries as well as MPI
260:        - provides summary and diagnostic information if certain runtime
261:          options are chosen (e.g., -log_view).
262:   */
263:   PetscFinalize();
264:   return 0;
265: }

267: /*TEST

269:    test:
270:       TODO: Matrix row/column sizes are not compatible with block size
271:       suffix: 1
272:       nsize: 4
273:       output_file: output/ex10_1.out
274:       requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
275:       args: -f0 ${DATAFILESPATH}/matrices/medium -f1 ${DATAFILESPATH}/matrices/arco6 -ksp_gmres_classicalgramschmidt -mat_type baij -matload_block_size 3 -pc_type bjacobi

277:    test:
278:       TODO: Matrix row/column sizes are not compatible with block size
279:       suffix: 2
280:       nsize: 4
281:       output_file: output/ex10_2.out
282:       requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
283:       args: -f0 ${DATAFILESPATH}/matrices/medium -f1 ${DATAFILESPATH}/matrices/arco6 -ksp_gmres_classicalgramschmidt -mat_type baij -matload_block_size 3 -pc_type bjacobi -trans

285:    test:
286:       suffix: 3
287:       requires: double complex !defined(PETSC_USE_64BIT_INDICES)
288:       args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/nh-complex-int32-float64 -ksp_type bicg

290:    test:
291:       suffix: 4
292:       args: -f ${DATAFILESPATH}/matrices/medium -ksp_type bicg -permute rcm
293:       requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)

295: TEST*/