Actual source code: ex30.c

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

  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/examples/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: */
 23: /*T
 24:    Concepts: KSP solving a linear system
 25:    Processors: n
 26: T*/

 28:  #include <petscksp.h>

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

 47:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
 48:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 49:   PetscOptionsGetBool(NULL,NULL,"-table",&table,NULL);
 50:   PetscOptionsGetBool(NULL,NULL,"-trans",&trans,NULL);
 51:   PetscOptionsGetBool(NULL,NULL,"-partition",&partition,NULL);
 52:   PetscOptionsGetBool(NULL,NULL,"-initialguess",&initialguess,NULL);
 53:   PetscOptionsGetBool(NULL,NULL,"-output_solution",&outputSoln,NULL);
 54:   PetscOptionsGetBool(NULL,NULL,"-ckrnorm",&ckrnorm,NULL);
 55:   PetscOptionsGetBool(NULL,NULL,"-ckerror",&ckerror,NULL);

 57:   /*
 58:      Determine files from which we read the two linear systems
 59:      (matrix and right-hand-side vector).
 60:   */
 61:   PetscOptionsGetString(NULL,NULL,"-f",file[0],PETSC_MAX_PATH_LEN,&flg);
 62:   if (flg) {
 63:     PetscStrcpy(file[1],file[0]);
 64:     preload = PETSC_FALSE;
 65:   } else {
 66:     PetscOptionsGetString(NULL,NULL,"-f0",file[0],PETSC_MAX_PATH_LEN,&flg);
 67:     if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate binary file with the -f0 or -f option");
 68:     PetscOptionsGetString(NULL,NULL,"-f1",file[1],PETSC_MAX_PATH_LEN,&flg);
 69:     if (!flg) preload = PETSC_FALSE;   /* don't bother with second system */
 70:   }

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

 86:   /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
 87:                          Load system
 88:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 90:   /*
 91:      Open binary file.  Note that we use FILE_MODE_READ to indicate
 92:      reading from this file.
 93:   */
 94:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[PetscPreLoadIt],FILE_MODE_READ,&fd);

 96:   /*
 97:      Load the matrix and vector; then destroy the viewer.
 98:   */
 99:   MatCreate(PETSC_COMM_WORLD,&A);
100:   MatSetFromOptions(A);
101:   MatLoad(A,fd);

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

122:   /* Test MatDuplicate() */
123:   if (Test_MatDuplicate) {
124:     MatDuplicate(A,MAT_COPY_VALUES,&B);
125:     MatEqual(A,B,&flg);
126:     if (!flg) {
127:       PetscPrintf(PETSC_COMM_WORLD,"  A != B \n");
128:     }
129:     MatDestroy(&B);
130:   }

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

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

173:   /* Check whether A is symmetric */
174:   flg  = PETSC_FALSE;
175:   PetscOptionsGetBool(NULL,NULL, "-check_symmetry", &flg,NULL);
176:   if (flg) {
177:     Mat Atrans;
178:     MatTranspose(A, MAT_INITIAL_MATRIX,&Atrans);
179:     MatEqual(A, Atrans, &isSymmetric);
180:     if (isSymmetric) {
181:       PetscPrintf(PETSC_COMM_WORLD,"A is symmetric \n");
182:     } else {
183:       PetscPrintf(PETSC_COMM_WORLD,"A is non-symmetric \n");
184:     }
185:     MatDestroy(&Atrans);
186:   }

188:   VecDuplicate(b,&b2);
189:   VecDuplicate(b,&x);
190:   PetscObjectSetName((PetscObject)x, "Solution vector");
191:   VecDuplicate(b,&u);
192:   PetscObjectSetName((PetscObject)u, "True Solution vector");
193:   VecSet(x,0.0);

195:   if (ckerror) {   /* Set true solution */
196:     VecSet(u,1.0);
197:     MatMult(A,u,b);
198:   }

200:   /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
201:                     Setup solve for system
202:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

228:     /* need to move the vector also */
229:     ISDestroy(&is);
230:     MatDestroy(&A);
231:     A    = BB;
232:   }

234:   /*
235:      Create linear solver; set operators; set runtime options.
236:   */
237:   KSPCreate(PETSC_COMM_WORLD,&ksp);
238:   KSPSetInitialGuessNonzero(ksp,initialguess);
239:   num_numfac = 1;
240:   PetscOptionsGetInt(NULL,NULL,"-num_numfac",&num_numfac,NULL);
241:   while (num_numfac--) {

243:     KSPSetOperators(ksp,A,A);
244:     KSPSetFromOptions(ksp);

246:     /*
247:      Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
248:      enable more precise profiling of setting up the preconditioner.
249:      These calls are optional, since both will be called within
250:      KSPSolve() if they haven't been called already.
251:     */
252:     KSPSetUp(ksp);
253:     KSPSetUpOnBlocks(ksp);

255:     /*
256:      Tests "diagonal-scaling of preconditioned residual norm" as used
257:      by many ODE integrator codes including SUNDIALS. Note this is different
258:      than diagonally scaling the matrix before computing the preconditioner
259:     */
260:     diagonalscale = PETSC_FALSE;
261:     PetscOptionsGetBool(NULL,NULL,"-diagonal_scale",&diagonalscale,NULL);
262:     if (diagonalscale) {
263:       PC       pc;
264:       PetscInt j,start,end,n;
265:       Vec      scale;

267:       KSPGetPC(ksp,&pc);
268:       VecGetSize(x,&n);
269:       VecDuplicate(x,&scale);
270:       VecGetOwnershipRange(scale,&start,&end);
271:       for (j=start; j<end; j++) {
272:         VecSetValue(scale,j,((PetscReal)(j+1))/((PetscReal)n),INSERT_VALUES);
273:       }
274:       VecAssemblyBegin(scale);
275:       VecAssemblyEnd(scale);
276:       PCSetDiagonalScale(pc,scale);
277:       VecDestroy(&scale);
278:     }

280:     /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
281:                          Solve system
282:       - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
283:     /*
284:      Solve linear system; 
285:     */
286:     if (trans) {
287:       KSPSolveTranspose(ksp,b,x);
288:       KSPGetIterationNumber(ksp,&its);
289:     } else {
290:       PetscInt num_rhs=1;
291:       PetscOptionsGetInt(NULL,NULL,"-num_rhs",&num_rhs,NULL);

293:       while (num_rhs--) {
294:         KSPSolve(ksp,b,x);
295:       }
296:       KSPGetIterationNumber(ksp,&its);
297:       if (ckrnorm) {     /* Check residual for each rhs */
298:         if (trans) {
299:           MatMultTranspose(A,x,b2);
300:         } else {
301:           MatMult(A,x,b2);
302:         }
303:         VecAXPY(b2,-1.0,b);
304:         VecNorm(b2,NORM_2,&rnorm);
305:         PetscPrintf(PETSC_COMM_WORLD,"  Number of iterations = %3D\n",its);
306:         PetscPrintf(PETSC_COMM_WORLD,"  Residual norm %g\n",(double)rnorm);
307:       }
308:       if (ckerror && !trans) {    /* Check error for each rhs */
309:         /* VecView(x,PETSC_VIEWER_STDOUT_WORLD); */
310:         VecAXPY(u,-1.0,x);
311:         VecNorm(u,NORM_2,&enorm);
312:         PetscPrintf(PETSC_COMM_WORLD,"  Error norm %g\n",(double)enorm);
313:       }

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


318:     /*
319:      Write output (optinally using table for solver details).
320:       - PetscPrintf() handles output for multiprocessor jobs
321:         by printing from only one processor in the communicator.
322:       - KSPView() prints information about the linear solver.
323:     */
324:     if (table && ckrnorm) {
325:       char        *matrixname,kspinfo[120];
326:       PetscViewer viewer;

328:       /*
329:         Open a string viewer; then write info to it.
330:       */
331:       PetscViewerStringOpen(PETSC_COMM_WORLD,kspinfo,sizeof(kspinfo),&viewer);
332:       KSPView(ksp,viewer);
333:       PetscStrrchr(file[PetscPreLoadIt],'/',&matrixname);
334:       PetscPrintf(PETSC_COMM_WORLD,"%-8.8s %3D %2.0e %s \n", matrixname,its,rnorm,kspinfo);

336:       /*
337:         Destroy the viewer
338:       */
339:       PetscViewerDestroy(&viewer);
340:     }

342:     PetscOptionsGetString(NULL,NULL,"-solution",file[3],PETSC_MAX_PATH_LEN,&flg);
343:     if (flg) {
344:       PetscViewer viewer;
345:       Vec         xstar;

347:       PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[3],FILE_MODE_READ,&viewer);
348:       VecCreate(PETSC_COMM_WORLD,&xstar);
349:       VecLoad(xstar,viewer);
350:       VecAXPY(xstar, -1.0, x);
351:       VecNorm(xstar, NORM_2, &enorm);
352:       PetscPrintf(PETSC_COMM_WORLD, "Error norm %g\n", (double)enorm);
353:       VecDestroy(&xstar);
354:       PetscViewerDestroy(&viewer);
355:     }
356:     if (outputSoln) {
357:       PetscViewer viewer;

359:       PetscViewerBinaryOpen(PETSC_COMM_WORLD,"solution.petsc",FILE_MODE_WRITE,&viewer);
360:       VecView(x, viewer);
361:       PetscViewerDestroy(&viewer);
362:     }

364:     flg  = PETSC_FALSE;
365:     PetscOptionsGetBool(NULL,NULL, "-ksp_reason", &flg,NULL);
366:     if (flg) {
367:       KSPConvergedReason reason;
368:       KSPGetConvergedReason(ksp,&reason);
369:       PetscPrintf(PETSC_COMM_WORLD,"KSPConvergedReason: %D\n", reason);
370:     }

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

374:   /*
375:      Free work space.  All PETSc objects should be destroyed when they
376:      are no longer needed.
377:   */
378:   MatDestroy(&A); VecDestroy(&b);
379:   VecDestroy(&u); VecDestroy(&x);
380:   VecDestroy(&b2);
381:   KSPDestroy(&ksp);
382:   if (flgB) { MatDestroy(&B); }
383:   PetscPreLoadEnd();
384:   /* -----------------------------------------------------------
385:                       End of linear solver loop
386:      ----------------------------------------------------------- */

388:   PetscFinalize();
389:   return ierr;
390: }

392: /*TEST

394:     test:
395:       args: -f ${DATAFILESPATH}/matrices/small -ksp_type preonly -pc_type ilu -pc_factor_mat_ordering_type natural -num_numfac 2 -pc_factor_reuse_fill
396:       requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
397:       output_file: output/ex30.out

399:     test:
400:       suffix: 2
401:       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
402:       requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)

404:     test:
405:       suffix: shiftnz
406:       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
407:       requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)

409:     test:
410:       suffix: shiftpd
411:       args: -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type lu -pc_factor_shift_type POSITIVE_DEFINITE
412:       requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)

414:     test:
415:       suffix: shift_cholesky_aij
416:       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 
417:       requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
418:       output_file: output/ex30_shiftnz.out

420:     test:
421:       suffix: shiftpd_2
422:       args: -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type cholesky -pc_factor_shift_type POSITIVE_DEFINITE
423:       requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)

425:     test:
426:       suffix: shift_cholesky_sbaij
427:       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
428:       requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
429:       output_file: output/ex30_shiftnz.out

431:     test:
432:       suffix: shiftpd_2_sbaij
433:       args: -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type cholesky -pc_factor_shift_type POSITIVE_DEFINITE -mat_type sbaij
434:       requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
435:       output_file: output/ex30_shiftpd_2.out

437:     test:
438:       suffix: shiftinblocks
439:       args:  -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type lu -pc_factor_shift_type INBLOCKS
440:       requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)

442:     test:
443:       suffix: shiftinblocks2
444:       args: -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type cholesky -pc_factor_shift_type INBLOCKS
445:       requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
446:       output_file: output/ex30_shiftinblocks.out

448:     test:
449:       suffix: shiftinblockssbaij
450:       args: -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type cholesky -pc_factor_shift_type INBLOCKS -mat_type sbaij
451:       requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
452:       output_file: output/ex30_shiftinblocks.out

454: TEST*/