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: */
 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;
 39:   PetscInt       its,num_numfac;
 40:   PetscReal      rnorm,enorm;
 41:   PetscBool      preload=PETSC_TRUE,diagonalscale,isSymmetric,ckrnorm=PETSC_TRUE,Test_MatDuplicate=PETSC_FALSE,ckerror=PETSC_FALSE;
 42:   PetscMPIInt    rank;
 43:   PetscScalar    sigma;
 44:   PetscInt       m;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

334:       /*
335:         Destroy the viewer
336:       */
337:       PetscViewerDestroy(&viewer);
338:     }

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

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

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

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

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

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

386:   PetscFinalize();
387:   return 0;
388: }

390: /*TEST

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

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

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

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

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

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

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

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

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

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

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

453: TEST*/