Actual source code: ex30.c


  2: static char help[] = "Reads a PETSc matrix and vector from a file and solves a linear system.\n\
  3: It is copied and intended to move dirty codes from ksp/tutorials/ex10.c and simplify ex10.c.\n\
  4:   Input parameters include\n\
  5:   -f0 <input_file> : first file to load (small system)\n\
  6:   -f1 <input_file> : second file to load (larger system)\n\n\
  7:   -trans  : solve transpose system instead\n\n";
  8: /*
  9:   This code  can be used to test PETSc interface to other packages.\n\
 10:   Examples of command line options:       \n\
 11:    ex30 -f0 <datafile> -ksp_type preonly  \n\
 12:         -help -ksp_view                  \n\
 13:         -num_numfac <num_numfac> -num_rhs <num_rhs> \n\
 14:         -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type superlu or superlu_dist or mumps \n\
 15:         -ksp_type preonly -pc_type cholesky -pc_factor_mat_solver_type mumps \n\
 16:    mpiexec -n <np> ex30 -f0 <datafile> -ksp_type cg -pc_type asm -pc_asm_type basic -sub_pc_type icc -mat_type sbaij

 18:    ./ex30 -f0 $D/small -mat_sigma -3.999999999999999 -ksp_type fgmres -pc_type lu -pc_factor_mat_solver_type superlu -mat_superlu_conditionnumber -ckerror -mat_superlu_diagpivotthresh 0
 19:    ./ex30 -f0 $D/small -mat_sigma -3.999999999999999 -ksp_type fgmres -pc_type hypre -pc_hypre_type boomeramg -ksp_type fgmres -ckError
 20:    ./ex30 -f0 $D/small -mat_sigma -3.999999999999999 -ksp_type fgmres -pc_type lu -pc_factor_mat_solver_type petsc -pc_factor_shift_type NONZERO -pc_factor_shift_amount 1.e-5 -ckerror
 21:  \n\n";
 22: */

 24: #include <petscksp.h>

 26: int main(int argc, char **args)
 27: {
 28:   KSP         ksp;
 29:   Mat         A, B;
 30:   Vec         x, b, u, b2;                 /* approx solution, RHS, exact solution */
 31:   PetscViewer fd;                          /* viewer */
 32:   char        file[4][PETSC_MAX_PATH_LEN]; /* input file name */
 33:   PetscBool   table = PETSC_FALSE, flg, flgB = PETSC_FALSE, trans = PETSC_FALSE, partition = PETSC_FALSE, initialguess = PETSC_FALSE;
 34:   PetscBool   outputSoln = PETSC_FALSE;
 35:   PetscInt    its, num_numfac;
 36:   PetscReal   rnorm, enorm;
 37:   PetscBool   preload = PETSC_TRUE, diagonalscale, isSymmetric, ckrnorm = PETSC_TRUE, Test_MatDuplicate = PETSC_FALSE, ckerror = PETSC_FALSE;
 38:   PetscMPIInt rank;
 39:   PetscScalar sigma;
 40:   PetscInt    m;

 43:   PetscInitialize(&argc, &args, (char *)0, help);
 44:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
 45:   PetscOptionsGetBool(NULL, NULL, "-table", &table, NULL);
 46:   PetscOptionsGetBool(NULL, NULL, "-trans", &trans, NULL);
 47:   PetscOptionsGetBool(NULL, NULL, "-partition", &partition, NULL);
 48:   PetscOptionsGetBool(NULL, NULL, "-initialguess", &initialguess, NULL);
 49:   PetscOptionsGetBool(NULL, NULL, "-output_solution", &outputSoln, NULL);
 50:   PetscOptionsGetBool(NULL, NULL, "-ckrnorm", &ckrnorm, NULL);
 51:   PetscOptionsGetBool(NULL, NULL, "-ckerror", &ckerror, NULL);

 53:   /*
 54:      Determine files from which we read the two linear systems
 55:      (matrix and right-hand-side vector).
 56:   */
 57:   PetscOptionsGetString(NULL, NULL, "-f", file[0], sizeof(file[0]), &flg);
 58:   if (flg) {
 59:     PetscStrcpy(file[1], file[0]);
 60:     preload = PETSC_FALSE;
 61:   } else {
 62:     PetscOptionsGetString(NULL, NULL, "-f0", file[0], sizeof(file[0]), &flg);
 64:     PetscOptionsGetString(NULL, NULL, "-f1", file[1], sizeof(file[1]), &flg);
 65:     if (!flg) preload = PETSC_FALSE; /* don't bother with second system */
 66:   }

 68:   /* -----------------------------------------------------------
 69:                   Beginning of linear solver loop
 70:      ----------------------------------------------------------- */
 71:   /*
 72:      Loop through the linear solve 2 times.
 73:       - The intention here is to preload and solve a small system;
 74:         then load another (larger) system and solve it as well.
 75:         This process preloads the instructions with the smaller
 76:         system so that more accurate performance monitoring (via
 77:         -log_view) can be done with the larger one (that actually
 78:         is the system of interest).
 79:   */
 80:   PetscPreLoadBegin(preload, "Load system");

 82:   /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
 83:                          Load system
 84:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 86:   /*
 87:      Open binary file.  Note that we use FILE_MODE_READ to indicate
 88:      reading from this file.
 89:   */
 90:   PetscViewerBinaryOpen(PETSC_COMM_WORLD, file[PetscPreLoadIt], FILE_MODE_READ, &fd);

 92:   /*
 93:      Load the matrix and vector; then destroy the viewer.
 94:   */
 95:   MatCreate(PETSC_COMM_WORLD, &A);
 96:   MatSetFromOptions(A);
 97:   MatLoad(A, fd);

 99:   flg = PETSC_FALSE;
100:   PetscOptionsGetString(NULL, NULL, "-rhs", file[2], sizeof(file[2]), &flg);
101:   if (flg) { /* rhs is stored in a separate file */
102:     PetscViewerDestroy(&fd);
103:     PetscViewerBinaryOpen(PETSC_COMM_WORLD, file[2], FILE_MODE_READ, &fd);
104:     VecCreate(PETSC_COMM_WORLD, &b);
105:     VecLoad(b, fd);
106:   } else {
107:     /* if file contains no RHS, then use a vector of all ones */
108:     PetscInfo(0, "Using vector of ones for RHS\n");
109:     MatGetLocalSize(A, &m, NULL);
110:     VecCreate(PETSC_COMM_WORLD, &b);
111:     VecSetSizes(b, m, PETSC_DECIDE);
112:     VecSetFromOptions(b);
113:     VecSet(b, 1.0);
114:     PetscObjectSetName((PetscObject)b, "Rhs vector");
115:   }
116:   PetscViewerDestroy(&fd);

118:   /* Test MatDuplicate() */
119:   if (Test_MatDuplicate) {
120:     MatDuplicate(A, MAT_COPY_VALUES, &B);
121:     MatEqual(A, B, &flg);
122:     if (!flg) PetscPrintf(PETSC_COMM_WORLD, "  A != B \n");
123:     MatDestroy(&B);
124:   }

126:   /* Add a shift to A */
127:   PetscOptionsGetScalar(NULL, NULL, "-mat_sigma", &sigma, &flg);
128:   if (flg) {
129:     PetscOptionsGetString(NULL, NULL, "-fB", file[2], sizeof(file[2]), &flgB);
130:     if (flgB) {
131:       /* load B to get A = A + sigma*B */
132:       PetscViewerBinaryOpen(PETSC_COMM_WORLD, file[2], FILE_MODE_READ, &fd);
133:       MatCreate(PETSC_COMM_WORLD, &B);
134:       MatSetOptionsPrefix(B, "B_");
135:       MatLoad(B, fd);
136:       PetscViewerDestroy(&fd);
137:       MatAXPY(A, sigma, B, DIFFERENT_NONZERO_PATTERN); /* A <- sigma*B + A */
138:     } else {
139:       MatShift(A, sigma);
140:     }
141:   }

143:   /* Make A singular for testing zero-pivot of ilu factorization        */
144:   /* Example: ./ex30 -f0 <datafile> -test_zeropivot -set_row_zero -pc_factor_shift_nonzero */
145:   flg = PETSC_FALSE;
146:   PetscOptionsGetBool(NULL, NULL, "-test_zeropivot", &flg, NULL);
147:   if (flg) {
148:     PetscInt           row, ncols;
149:     const PetscInt    *cols;
150:     const PetscScalar *vals;
151:     PetscBool          flg1 = PETSC_FALSE;
152:     PetscScalar       *zeros;
153:     row = 0;
154:     MatGetRow(A, row, &ncols, &cols, &vals);
155:     PetscCalloc1(ncols + 1, &zeros);
156:     flg1 = PETSC_FALSE;
157:     PetscOptionsGetBool(NULL, NULL, "-set_row_zero", &flg1, NULL);
158:     if (flg1) { /* set entire row as zero */
159:       MatSetValues(A, 1, &row, ncols, cols, zeros, INSERT_VALUES);
160:     } else { /* only set (row,row) entry as zero */
161:       MatSetValues(A, 1, &row, 1, &row, zeros, INSERT_VALUES);
162:     }
163:     MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
164:     MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
165:   }

167:   /* Check whether A is symmetric */
168:   flg = PETSC_FALSE;
169:   PetscOptionsGetBool(NULL, NULL, "-check_symmetry", &flg, NULL);
170:   if (flg) {
171:     Mat Atrans;
172:     MatTranspose(A, MAT_INITIAL_MATRIX, &Atrans);
173:     MatEqual(A, Atrans, &isSymmetric);
174:     if (isSymmetric) {
175:       PetscPrintf(PETSC_COMM_WORLD, "A is symmetric \n");
176:     } else {
177:       PetscPrintf(PETSC_COMM_WORLD, "A is non-symmetric \n");
178:     }
179:     MatDestroy(&Atrans);
180:   }

182:   VecDuplicate(b, &b2);
183:   VecDuplicate(b, &x);
184:   PetscObjectSetName((PetscObject)x, "Solution vector");
185:   VecDuplicate(b, &u);
186:   PetscObjectSetName((PetscObject)u, "True Solution vector");
187:   VecSet(x, 0.0);

189:   if (ckerror) { /* Set true solution */
190:     VecSet(u, 1.0);
191:     MatMult(A, u, b);
192:   }

194:   /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
195:                     Setup solve for system
196:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

198:   if (partition) {
199:     MatPartitioning mpart;
200:     IS              mis, nis, is;
201:     PetscInt       *count;
202:     PetscMPIInt     size;
203:     Mat             BB;
204:     MPI_Comm_size(PETSC_COMM_WORLD, &size);
205:     MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
206:     PetscMalloc1(size, &count);
207:     MatPartitioningCreate(PETSC_COMM_WORLD, &mpart);
208:     MatPartitioningSetAdjacency(mpart, A);
209:     /* MatPartitioningSetVertexWeights(mpart, weight); */
210:     MatPartitioningSetFromOptions(mpart);
211:     MatPartitioningApply(mpart, &mis);
212:     MatPartitioningDestroy(&mpart);
213:     ISPartitioningToNumbering(mis, &nis);
214:     ISPartitioningCount(mis, size, count);
215:     ISDestroy(&mis);
216:     ISInvertPermutation(nis, count[rank], &is);
217:     PetscFree(count);
218:     ISDestroy(&nis);
219:     ISSort(is);
220:     MatCreateSubMatrix(A, is, is, MAT_INITIAL_MATRIX, &BB);

222:     /* need to move the vector also */
223:     ISDestroy(&is);
224:     MatDestroy(&A);
225:     A = BB;
226:   }

228:   /*
229:      Create linear solver; set operators; set runtime options.
230:   */
231:   KSPCreate(PETSC_COMM_WORLD, &ksp);
232:   KSPSetInitialGuessNonzero(ksp, initialguess);
233:   num_numfac = 1;
234:   PetscOptionsGetInt(NULL, NULL, "-num_numfac", &num_numfac, NULL);
235:   while (num_numfac--) {
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:     /*
249:      Tests "diagonal-scaling of preconditioned residual norm" as used
250:      by many ODE integrator codes including SUNDIALS. Note this is different
251:      than diagonally scaling the matrix before computing the preconditioner
252:     */
253:     diagonalscale = PETSC_FALSE;
254:     PetscOptionsGetBool(NULL, NULL, "-diagonal_scale", &diagonalscale, NULL);
255:     if (diagonalscale) {
256:       PC       pc;
257:       PetscInt j, start, end, n;
258:       Vec      scale;

260:       KSPGetPC(ksp, &pc);
261:       VecGetSize(x, &n);
262:       VecDuplicate(x, &scale);
263:       VecGetOwnershipRange(scale, &start, &end);
264:       for (j = start; j < end; j++) VecSetValue(scale, j, ((PetscReal)(j + 1)) / ((PetscReal)n), INSERT_VALUES);
265:       VecAssemblyBegin(scale);
266:       VecAssemblyEnd(scale);
267:       PCSetDiagonalScale(pc, scale);
268:       VecDestroy(&scale);
269:     }

271:     /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
272:                          Solve system
273:       - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
274:     /*
275:      Solve linear system;
276:     */
277:     if (trans) {
278:       KSPSolveTranspose(ksp, b, x);
279:       KSPGetIterationNumber(ksp, &its);
280:     } else {
281:       PetscInt num_rhs = 1;
282:       PetscOptionsGetInt(NULL, NULL, "-num_rhs", &num_rhs, NULL);

284:       while (num_rhs--) KSPSolve(ksp, b, x);
285:       KSPGetIterationNumber(ksp, &its);
286:       if (ckrnorm) { /* Check residual for each rhs */
287:         if (trans) {
288:           MatMultTranspose(A, x, b2);
289:         } else {
290:           MatMult(A, x, b2);
291:         }
292:         VecAXPY(b2, -1.0, b);
293:         VecNorm(b2, NORM_2, &rnorm);
294:         PetscPrintf(PETSC_COMM_WORLD, "  Number of iterations = %3" PetscInt_FMT "\n", its);
295:         PetscPrintf(PETSC_COMM_WORLD, "  Residual norm %g\n", (double)rnorm);
296:       }
297:       if (ckerror && !trans) { /* Check error for each rhs */
298:         /* VecView(x,PETSC_VIEWER_STDOUT_WORLD); */
299:         VecAXPY(u, -1.0, x);
300:         VecNorm(u, NORM_2, &enorm);
301:         PetscPrintf(PETSC_COMM_WORLD, "  Error norm %g\n", (double)enorm);
302:       }

304:     } /* while (num_rhs--) */

306:     /*
307:      Write output (optionally using table for solver details).
308:       - PetscPrintf() handles output for multiprocessor jobs
309:         by printing from only one processor in the communicator.
310:       - KSPView() prints information about the linear solver.
311:     */
312:     if (table && ckrnorm) {
313:       char       *matrixname, kspinfo[120];
314:       PetscViewer viewer;

316:       /*
317:         Open a string viewer; then write info to it.
318:       */
319:       PetscViewerStringOpen(PETSC_COMM_WORLD, kspinfo, sizeof(kspinfo), &viewer);
320:       KSPView(ksp, viewer);
321:       PetscStrrchr(file[PetscPreLoadIt], '/', &matrixname);
322:       PetscPrintf(PETSC_COMM_WORLD, "%-8.8s %3" PetscInt_FMT " %2.0e %s \n", matrixname, its, (double)rnorm, kspinfo);

324:       /*
325:         Destroy the viewer
326:       */
327:       PetscViewerDestroy(&viewer);
328:     }

330:     PetscOptionsGetString(NULL, NULL, "-solution", file[3], sizeof(file[3]), &flg);
331:     if (flg) {
332:       PetscViewer viewer;
333:       Vec         xstar;

335:       PetscViewerBinaryOpen(PETSC_COMM_WORLD, file[3], FILE_MODE_READ, &viewer);
336:       VecCreate(PETSC_COMM_WORLD, &xstar);
337:       VecLoad(xstar, viewer);
338:       VecAXPY(xstar, -1.0, x);
339:       VecNorm(xstar, NORM_2, &enorm);
340:       PetscPrintf(PETSC_COMM_WORLD, "Error norm %g\n", (double)enorm);
341:       VecDestroy(&xstar);
342:       PetscViewerDestroy(&viewer);
343:     }
344:     if (outputSoln) {
345:       PetscViewer viewer;

347:       PetscViewerBinaryOpen(PETSC_COMM_WORLD, "solution.petsc", FILE_MODE_WRITE, &viewer);
348:       VecView(x, viewer);
349:       PetscViewerDestroy(&viewer);
350:     }

352:     flg = PETSC_FALSE;
353:     PetscOptionsGetBool(NULL, NULL, "-ksp_reason", &flg, NULL);
354:     if (flg) {
355:       KSPConvergedReason reason;
356:       KSPGetConvergedReason(ksp, &reason);
357:       PetscPrintf(PETSC_COMM_WORLD, "KSPConvergedReason: %s\n", KSPConvergedReasons[reason]);
358:     }

360:   } /* while (num_numfac--) */

362:   /*
363:      Free work space.  All PETSc objects should be destroyed when they
364:      are no longer needed.
365:   */
366:   MatDestroy(&A);
367:   VecDestroy(&b);
368:   VecDestroy(&u);
369:   VecDestroy(&x);
370:   VecDestroy(&b2);
371:   KSPDestroy(&ksp);
372:   if (flgB) MatDestroy(&B);
373:   PetscPreLoadEnd();
374:   /* -----------------------------------------------------------
375:                       End of linear solver loop
376:      ----------------------------------------------------------- */

378:   PetscFinalize();
379:   return 0;
380: }

382: /*TEST

384:     test:
385:       args: -f ${DATAFILESPATH}/matrices/small -ksp_type preonly -pc_type ilu -pc_factor_mat_ordering_type natural -num_numfac 2 -pc_factor_reuse_fill
386:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
387:       output_file: output/ex30.out

389:     test:
390:       TODO: Matrix row/column sizes are not compatible with block size
391:       suffix: 2
392:       args: -f ${DATAFILESPATH}/matrices/arco1 -mat_type baij -matload_block_size 3 -ksp_type preonly -pc_type ilu -pc_factor_mat_ordering_type natural -num_numfac 2
393:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)

395:     test:
396:       suffix: shiftnz
397:       args: -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type lu -pc_factor_shift_type NONZERO -pc_factor_shift_amount 1.e-5
398:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)

400:     test:
401:       suffix: shiftpd
402:       args: -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type lu -pc_factor_shift_type POSITIVE_DEFINITE
403:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)

405:     test:
406:       suffix: shift_cholesky_aij
407:       args: -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type cholesky -pc_factor_shift_type NONZERO -pc_factor_shift_amount 1.e-5
408:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
409:       output_file: output/ex30_shiftnz.out

411:     test:
412:       suffix: shiftpd_2
413:       args: -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type cholesky -pc_factor_shift_type POSITIVE_DEFINITE
414:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)

416:     test:
417:       suffix: shift_cholesky_sbaij
418:       args: -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type cholesky -pc_factor_shift_type NONZERO -pc_factor_shift_amount 1.e-5 -mat_type sbaij
419:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
420:       output_file: output/ex30_shiftnz.out

422:     test:
423:       suffix: shiftpd_2_sbaij
424:       args: -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type cholesky -pc_factor_shift_type POSITIVE_DEFINITE -mat_type sbaij
425:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
426:       output_file: output/ex30_shiftpd_2.out

428:     test:
429:       suffix: shiftinblocks
430:       args:  -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type lu -pc_factor_shift_type INBLOCKS
431:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)

433:     test:
434:       suffix: shiftinblocks2
435:       args: -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type cholesky -pc_factor_shift_type INBLOCKS
436:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
437:       output_file: output/ex30_shiftinblocks.out

439:     test:
440:       suffix: shiftinblockssbaij
441:       args: -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type cholesky -pc_factor_shift_type INBLOCKS -mat_type sbaij
442:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
443:       output_file: output/ex30_shiftinblocks.out

445: TEST*/