Actual source code: ex30.c
petsc-3.12.4 2020-02-04
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*/