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;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

229:   /*
230:      Create linear solver; set operators; set runtime options.
231:   */
232:   KSPCreate(PETSC_COMM_WORLD,&ksp);
233:   KSPSetInitialGuessNonzero(ksp,initialguess);
234:   num_numfac = 1;
235:   PetscOptionsGetInt(NULL,NULL,"-num_numfac",&num_numfac,NULL);
236:   while (num_numfac--) {

238:     KSPSetOperators(ksp,A,A);
239:     KSPSetFromOptions(ksp);

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

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

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

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

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

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

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

322:       /*
323:         Open a string viewer; then write info to it.
324:       */
325:       PetscViewerStringOpen(PETSC_COMM_WORLD,kspinfo,sizeof(kspinfo),&viewer);
326:       KSPView(ksp,viewer);
327:       PetscStrrchr(file[PetscPreLoadIt],'/',&matrixname);
328:       PetscPrintf(PETSC_COMM_WORLD,"%-8.8s %3D %2.0e %s \n", matrixname,its,rnorm,kspinfo);

330:       /*
331:         Destroy the viewer
332:       */
333:       PetscViewerDestroy(&viewer);
334:     }

336:     PetscOptionsGetString(NULL,NULL,"-solution",file[3],sizeof(file[3]),&flg);
337:     if (flg) {
338:       PetscViewer viewer;
339:       Vec         xstar;

341:       PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[3],FILE_MODE_READ,&viewer);
342:       VecCreate(PETSC_COMM_WORLD,&xstar);
343:       VecLoad(xstar,viewer);
344:       VecAXPY(xstar, -1.0, x);
345:       VecNorm(xstar, NORM_2, &enorm);
346:       PetscPrintf(PETSC_COMM_WORLD, "Error norm %g\n", (double)enorm);
347:       VecDestroy(&xstar);
348:       PetscViewerDestroy(&viewer);
349:     }
350:     if (outputSoln) {
351:       PetscViewer viewer;

353:       PetscViewerBinaryOpen(PETSC_COMM_WORLD,"solution.petsc",FILE_MODE_WRITE,&viewer);
354:       VecView(x, viewer);
355:       PetscViewerDestroy(&viewer);
356:     }

358:     flg  = PETSC_FALSE;
359:     PetscOptionsGetBool(NULL,NULL, "-ksp_reason", &flg,NULL);
360:     if (flg) {
361:       KSPConvergedReason reason;
362:       KSPGetConvergedReason(ksp,&reason);
363:       PetscPrintf(PETSC_COMM_WORLD,"KSPConvergedReason: %D\n", reason);
364:     }

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

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

382:   PetscFinalize();
383:   return 0;
384: }

386: /*TEST

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

393:     test:
394:       TODO: Matrix row/column sizes are not compatible with block size
395:       suffix: 2
396:       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
397:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)

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

404:     test:
405:       suffix: shiftpd
406:       args: -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type lu -pc_factor_shift_type POSITIVE_DEFINITE
407:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)

409:     test:
410:       suffix: shift_cholesky_aij
411:       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
412:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
413:       output_file: output/ex30_shiftnz.out

415:     test:
416:       suffix: shiftpd_2
417:       args: -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type cholesky -pc_factor_shift_type POSITIVE_DEFINITE
418:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)

420:     test:
421:       suffix: shift_cholesky_sbaij
422:       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
423:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
424:       output_file: output/ex30_shiftnz.out

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

432:     test:
433:       suffix: shiftinblocks
434:       args:  -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type lu -pc_factor_shift_type INBLOCKS
435:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)

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

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

449: TEST*/