Actual source code: ex10.c
petsc-3.4.2 2013-07-02

      2: static char help[] = "Reads a PETSc matrix and vector from a file and solves a linear system.\n\
      3: This version first preloads and solves a small system, then loads \n\
      4: another (larger) system and solves it as well.  This example illustrates\n\
      5: preloading of instructions with the smaller system so that more accurate\n\
      6: performance monitoring can be done with the larger one (that actually\n\
      7: is the system of interest).  See the 'Performance Hints' chapter of the\n\
      8: users manual for a discussion of preloading.  Input parameters include\n\
      9:   -f0 <input_file> : first file to load (small system)\n\
     10:   -f1 <input_file> : second file to load (larger system)\n\n\
     11:   -trans  : solve transpose system instead\n\n";
     12: /*
     13:   This code can be used to test PETSc interface to other packages.\n\
     14:   Examples of command line options:       \n\
     15:    ./ex10 -f0 <datafile> -ksp_type preonly  \n\
     16:         -help -ksp_view                  \n\
     17:         -num_numfac <num_numfac> -num_rhs <num_rhs> \n\
     18:         -ksp_type preonly -pc_type lu -pc_factor_mat_solver_package superlu or superlu_dist or mumps \n\
     19:         -ksp_type preonly -pc_type cholesky -pc_factor_mat_solver_package mumps \n\
     20:    mpiexec -n <np> ./ex10 -f0 <datafile> -ksp_type cg -pc_type asm -pc_asm_type basic -sub_pc_type icc -mat_type sbaij
     21:  \n\n";
     22: */
     23: /*T
     24:    Concepts: KSP^solving a linear system
     25:    Processors: n
     26: T*/

     28: /*
     29:   Include "petscksp.h" so that we can use KSP solvers.  Note that this file
     30:   automatically includes:
     31:      petscsys.h       - base PETSc routines   petscvec.h - vectors
     32:      petscmat.h - matrices
     33:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
     34:      petscviewer.h - viewers               petscpc.h  - preconditioners
     35: */
     36: #include <petscksp.h>
     37: #include <petsctime.h>

     41: int main(int argc,char **args)
     42: {
     43:   KSP            ksp;             /* linear solver context */
     44:   Mat            A;               /* matrix */
     45:   Vec            x,b,u;           /* approx solution, RHS, exact solution */
     46:   PetscViewer    fd;              /* viewer */
     47:   char           file[4][PETSC_MAX_PATH_LEN];     /* input file name */
     48:   PetscBool      table     =PETSC_FALSE,flg,trans=PETSC_FALSE,initialguess = PETSC_FALSE;
     49:   PetscBool      outputSoln=PETSC_FALSE;
     51:   PetscInt       its,num_numfac,m,n,M;
     52:   PetscReal      norm;
     53:   PetscLogDouble tsetup,tsetup1,tsetup2,tsolve,tsolve1,tsolve2;
     54:   PetscBool      preload=PETSC_TRUE,isSymmetric,cknorm=PETSC_FALSE,initialguessfile = PETSC_FALSE;
     55:   PetscMPIInt    rank;
     56:   char           initialguessfilename[PETSC_MAX_PATH_LEN];

     58:   PetscInitialize(&argc,&args,(char*)0,help);
     59:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
     60:   PetscOptionsGetBool(NULL,"-table",&table,NULL);
     61:   PetscOptionsGetBool(NULL,"-trans",&trans,NULL);
     62:   PetscOptionsGetBool(NULL,"-initialguess",&initialguess,NULL);
     63:   PetscOptionsGetBool(NULL,"-output_solution",&outputSoln,NULL);
     64:   PetscOptionsGetString(NULL,"-initialguessfilename",initialguessfilename,PETSC_MAX_PATH_LEN,&initialguessfile);

     66:   /*
     67:      Determine files from which we read the two linear systems
     68:      (matrix and right-hand-side vector).
     69:   */
     70:   PetscOptionsGetString(NULL,"-f",file[0],PETSC_MAX_PATH_LEN,&flg);
     71:   if (flg) {
     72:     PetscStrcpy(file[1],file[0]);
     73:     preload = PETSC_FALSE;
     74:   } else {
     75:     PetscOptionsGetString(NULL,"-f0",file[0],PETSC_MAX_PATH_LEN,&flg);
     76:     if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate binary file with the -f0 or -f option");
     77:     PetscOptionsGetString(NULL,"-f1",file[1],PETSC_MAX_PATH_LEN,&flg);
     78:     if (!flg) preload = PETSC_FALSE;   /* don't bother with second system */
     79:   }

     81:   /* -----------------------------------------------------------
     82:                   Beginning of linear solver loop
     83:      ----------------------------------------------------------- */
     84:   /*
     85:      Loop through the linear solve 2 times.
     86:       - The intention here is to preload and solve a small system;
     87:         then load another (larger) system and solve it as well.
     88:         This process preloads the instructions with the smaller
     89:         system so that more accurate performance monitoring (via
     90:         -log_summary) can be done with the larger one (that actually
     91:         is the system of interest).
     92:   */
     93:   PetscPreLoadBegin(preload,"Load system");

     95:   /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
     96:                          Load system
     97:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     99:   /*
    100:      Open binary file.  Note that we use FILE_MODE_READ to indicate
    101:      reading from this file.
    102:   */
    103:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[PetscPreLoadIt],FILE_MODE_READ,&fd);

    105:   /*
    106:      Load the matrix and vector; then destroy the viewer.
    107:   */
    108:   MatCreate(PETSC_COMM_WORLD,&A);
    109:   MatSetFromOptions(A);
    110:   MatLoad(A,fd);

    112:   flg  = PETSC_FALSE;
    113:   PetscOptionsGetString(NULL,"-rhs",file[2],PETSC_MAX_PATH_LEN,&flg);
    114:   VecCreate(PETSC_COMM_WORLD,&b);
    115:   if (flg) {   /* rhs is stored in a separate file */
    116:     if (file[2][0] == '0') {
    117:       PetscInt    m;
    118:       PetscScalar one = 1.0;
    119:       PetscInfo(0,"Using vector of ones for RHS\n");
    120:       MatGetLocalSize(A,&m,NULL);
    121:       VecSetSizes(b,m,PETSC_DECIDE);
    122:       VecSetFromOptions(b);
    123:       VecSet(b,one);
    124:     } else {
    125:       PetscViewerDestroy(&fd);
    126:       PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[2],FILE_MODE_READ,&fd);
    127:       VecSetFromOptions(b);
    128:       VecLoad(b,fd);
    129:     }
    130:   } else {   /* rhs is stored in the same file as matrix */
    131:     VecSetFromOptions(b);
    132:     VecLoad(b,fd);
    133:   }
    134:   PetscViewerDestroy(&fd);

    136:   /* Make A singular for testing zero-pivot of ilu factorization        */
    137:   /* Example: ./ex10 -f0 <datafile> -test_zeropivot -set_row_zero -pc_factor_shift_nonzero */
    138:   flg  = PETSC_FALSE;
    139:   PetscOptionsGetBool(NULL, "-test_zeropivot", &flg,NULL);
    140:   if (flg) {
    141:     PetscInt          row,ncols;
    142:     const PetscInt    *cols;
    143:     const PetscScalar *vals;
    144:     PetscBool         flg1=PETSC_FALSE;
    145:     PetscScalar       *zeros;
    146:     row  = 0;
    147:     MatGetRow(A,row,&ncols,&cols,&vals);
    148:     PetscMalloc(sizeof(PetscScalar)*(ncols+1),&zeros);
    149:     PetscMemzero(zeros,(ncols+1)*sizeof(PetscScalar));
    150:     PetscOptionsGetBool(NULL, "-set_row_zero", &flg1,NULL);
    151:     if (flg1) {   /* set entire row as zero */
    152:       MatSetValues(A,1,&row,ncols,cols,zeros,INSERT_VALUES);
    153:     } else {   /* only set (row,row) entry as zero */
    154:       MatSetValues(A,1,&row,1,&row,zeros,INSERT_VALUES);
    155:     }
    156:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
    157:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
    158:   }

    160:   /* Check whether A is symmetric */
    161:   flg  = PETSC_FALSE;
    162:   PetscOptionsGetBool(NULL, "-check_symmetry", &flg,NULL);
    163:   if (flg) {
    164:     Mat Atrans;
    165:     MatTranspose(A, MAT_INITIAL_MATRIX,&Atrans);
    166:     MatEqual(A, Atrans, &isSymmetric);
    167:     if (isSymmetric) {
    168:       MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
    169:     } else {
    170:       PetscPrintf(PETSC_COMM_WORLD,"Warning: A is non-symmetric \n");
    171:     }
    172:     MatDestroy(&Atrans);
    173:   }

    175:   /*
    176:      If the loaded matrix is larger than the vector (due to being padded
    177:      to match the block size of the system), then create a new padded vector.
    178:   */

    180:   MatGetLocalSize(A,&m,&n);
    181:   /*  if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%d, %d)", m, n);*/
    182:   MatGetSize(A,&M,NULL);
    183:   VecGetSize(b,&m);
    184:   if (M != m) {   /* Create a new vector b by padding the old one */
    185:     PetscInt    j,mvec,start,end,indx;
    186:     Vec         tmp;
    187:     PetscScalar *bold;

    189:     VecCreate(PETSC_COMM_WORLD,&tmp);
    190:     VecSetSizes(tmp,n,PETSC_DECIDE);
    191:     VecSetFromOptions(tmp);
    192:     VecGetOwnershipRange(b,&start,&end);
    193:     VecGetLocalSize(b,&mvec);
    194:     VecGetArray(b,&bold);
    195:     for (j=0; j<mvec; j++) {
    196:       indx = start+j;
    197:       VecSetValues(tmp,1,&indx,bold+j,INSERT_VALUES);
    198:     }
    199:     VecRestoreArray(b,&bold);
    200:     VecDestroy(&b);
    201:     VecAssemblyBegin(tmp);
    202:     VecAssemblyEnd(tmp);
    203:     b    = tmp;
    204:   }

    206:   MatGetVecs(A,&x,NULL);
    207:   VecDuplicate(b,&u);
    208:   if (initialguessfile) {
    209:     PetscViewer viewer2;
    210:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,initialguessfilename,FILE_MODE_READ,&viewer2);
    211:     VecLoad(x,viewer2);
    212:     PetscViewerDestroy(&viewer2);
    213:     initialguess = PETSC_TRUE;
    214:   } else if (initialguess) {
    215:     VecSet(x,1.0);
    216:   } else {
    217:     VecSet(x,0.0);
    218:   }

    221:   /* Check scaling in A */
    222:   flg  = PETSC_FALSE;
    223:   PetscOptionsGetBool(NULL, "-check_scaling", &flg,NULL);
    224:   if (flg) {
    225:     Vec       max, min;
    226:     PetscInt  idx;
    227:     PetscReal val;

    229:     VecDuplicate(x, &max);
    230:     VecDuplicate(x, &min);
    231:     MatGetRowMaxAbs(A, max, NULL);
    232:     MatGetRowMinAbs(A, min, NULL);
    233:     {
    234:       PetscViewer viewer;

    236:       PetscViewerASCIIOpen(PETSC_COMM_WORLD, "max.data", &viewer);
    237:       VecView(max, viewer);
    238:       PetscViewerDestroy(&viewer);
    239:       PetscViewerASCIIOpen(PETSC_COMM_WORLD, "min.data", &viewer);
    240:       VecView(min, viewer);
    241:       PetscViewerDestroy(&viewer);
    242:     }
    243:     VecView(max, PETSC_VIEWER_DRAW_WORLD);
    244:     VecMax(max, &idx, &val);
    245:     PetscPrintf(PETSC_COMM_WORLD, "Largest max row element %G at row %d\n", val, idx);
    246:     VecView(min, PETSC_VIEWER_DRAW_WORLD);
    247:     VecMin(min, &idx, &val);
    248:     PetscPrintf(PETSC_COMM_WORLD, "Smallest min row element %G at row %d\n", val, idx);
    249:     VecMin(max, &idx, &val);
    250:     PetscPrintf(PETSC_COMM_WORLD, "Smallest max row element %G at row %d\n", val, idx);
    251:     VecPointwiseDivide(max, max, min);
    252:     VecMax(max, &idx, &val);
    253:     PetscPrintf(PETSC_COMM_WORLD, "Largest row ratio %G at row %d\n", val, idx);
    254:     VecView(max, PETSC_VIEWER_DRAW_WORLD);
    255:     VecDestroy(&max);
    256:     VecDestroy(&min);
    257:   }

    259:   /*  MatView(A,PETSC_VIEWER_STDOUT_WORLD); */
    260:   /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
    261:                     Setup solve for system
    262:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    263:   /*
    264:      Conclude profiling last stage; begin profiling next stage.
    265:   */
    266:   PetscPreLoadStage("KSPSetUpSolve");

    268:   /*
    269:      We also explicitly time this stage via PetscTime()
    270:   */
    271:   PetscTime(&tsetup1);

    273:   /*
    274:      Create linear solver; set operators; set runtime options.
    275:   */
    276:   KSPCreate(PETSC_COMM_WORLD,&ksp);
    277:   KSPSetInitialGuessNonzero(ksp,initialguess);
    278:   num_numfac = 1;
    279:   PetscOptionsGetInt(NULL,"-num_numfac",&num_numfac,NULL);
    280:   while (num_numfac--) {
    281:     PetscBool lsqr;
    282:     char      str[32];
    283:     PetscOptionsGetString(NULL,"-ksp_type",str,32,&lsqr);
    284:     if (lsqr) {
    285:       PetscStrcmp("lsqr",str,&lsqr);
    286:     }
    287:     if (lsqr) {
    288:       Mat BtB;
    289:       MatTransposeMatMult(A,A,MAT_INITIAL_MATRIX,4,&BtB);
    290:       KSPSetOperators(ksp,A,BtB,SAME_NONZERO_PATTERN);
    291:       MatDestroy(&BtB);
    292:     } else {
    293:       KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);
    294:     }
    295:     KSPSetFromOptions(ksp);

    297:     /*
    298:      Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
    299:      enable more precise profiling of setting up the preconditioner.
    300:      These calls are optional, since both will be called within
    301:      KSPSolve() if they haven't been called already.
    302:     */
    303:     KSPSetUp(ksp);
    304:     KSPSetUpOnBlocks(ksp);
    305:     PetscTime(&tsetup2);
    306:     tsetup = tsetup2 - tsetup1;

    308:     /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
    309:                          Solve system
    310:       - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    312:     /*
    313:      Solve linear system; we also explicitly time this stage.
    314:     */
    315:     PetscTime(&tsolve1);
    316:     if (trans) {
    317:       KSPSolveTranspose(ksp,b,x);
    318:       KSPGetIterationNumber(ksp,&its);
    319:     } else {
    320:       PetscInt num_rhs=1;
    321:       PetscOptionsGetInt(NULL,"-num_rhs",&num_rhs,NULL);
    322:       cknorm = PETSC_FALSE;
    323:       PetscOptionsGetBool(NULL,"-cknorm",&cknorm,NULL);
    324:       while (num_rhs--) {
    325:         if (num_rhs == 1) VecSet(x,0.0);
    326:         KSPSolve(ksp,b,x);
    327:       }
    328:       KSPGetIterationNumber(ksp,&its);
    329:       if (cknorm) {     /* Check error for each rhs */
    330:         if (trans) {
    331:           MatMultTranspose(A,x,u);
    332:         } else {
    333:           MatMult(A,x,u);
    334:         }
    335:         VecAXPY(u,-1.0,b);
    336:         VecNorm(u,NORM_2,&norm);
    337:         PetscPrintf(PETSC_COMM_WORLD,"  Number of iterations = %3D\n",its);
    338:         if (norm < 1.e-12) {
    339:           PetscPrintf(PETSC_COMM_WORLD,"  Residual norm < 1.e-12\n");
    340:         } else {
    341:           PetscPrintf(PETSC_COMM_WORLD,"  Residual norm %G\n",norm);
    342:         }
    343:       }
    344:     }   /* while (num_rhs--) */
    345:     PetscTime(&tsolve2);
    346:     tsolve = tsolve2 - tsolve1;

    348:     /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
    349:           Check error, print output, free data structures.
    350:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    352:     /*
    353:        Check error
    354:     */
    355:     if (trans) {
    356:       MatMultTranspose(A,x,u);
    357:     } else {
    358:       MatMult(A,x,u);
    359:     }
    360:     VecAXPY(u,-1.0,b);
    361:     VecNorm(u,NORM_2,&norm);
    362:     /*
    363:      Write output (optinally using table for solver details).
    364:       - PetscPrintf() handles output for multiprocessor jobs
    365:         by printing from only one processor in the communicator.
    366:       - KSPView() prints information about the linear solver.
    367:     */
    368:     if (table) {
    369:       char        *matrixname,kspinfo[120];
    370:       PetscViewer viewer;

    372:       /*
    373:        Open a string viewer; then write info to it.
    374:       */
    375:       PetscViewerStringOpen(PETSC_COMM_WORLD,kspinfo,120,&viewer);
    376:       KSPView(ksp,viewer);
    377:       PetscStrrchr(file[PetscPreLoadIt],'/',&matrixname);
    378:       PetscPrintf(PETSC_COMM_WORLD,"%-8.8s %3D %2.0e %2.1e %2.1e %2.1e %s \n",
    379:                          matrixname,its,norm,tsetup+tsolve,tsetup,tsolve,kspinfo);

    381:       /*
    382:         Destroy the viewer
    383:       */
    384:       PetscViewerDestroy(&viewer);
    385:     } else {
    386:       PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
    387:       if (norm < 1.e-12) {
    388:         PetscPrintf(PETSC_COMM_WORLD,"  Residual norm < 1.e-12\n");
    389:       } else {
    390:         PetscPrintf(PETSC_COMM_WORLD,"Residual norm %G\n",norm);
    391:       }
    392:     }
    393:     PetscOptionsGetString(NULL,"-solution",file[3],PETSC_MAX_PATH_LEN,&flg);
    394:     if (flg) {
    395:       PetscViewer viewer;
    396:       Vec         xstar;
    397:       PetscReal   norm;

    399:       PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[3],FILE_MODE_READ,&viewer);
    400:       VecCreate(PETSC_COMM_WORLD,&xstar);
    401:       VecLoad(xstar,viewer);
    402:       VecAXPY(xstar, -1.0, x);
    403:       VecNorm(xstar, NORM_2, &norm);
    404:       PetscPrintf(PETSC_COMM_WORLD, "Error norm %G\n", norm);
    405:       VecDestroy(&xstar);
    406:       PetscViewerDestroy(&viewer);
    407:     }
    408:     if (outputSoln) {
    409:       PetscViewer viewer;

    411:       PetscViewerBinaryOpen(PETSC_COMM_WORLD,"solution.petsc",FILE_MODE_WRITE,&viewer);
    412:       VecView(x, viewer);
    413:       PetscViewerDestroy(&viewer);
    414:     }

    416:     flg  = PETSC_FALSE;
    417:     PetscOptionsGetBool(NULL, "-ksp_reason", &flg,NULL);
    418:     if (flg) {
    419:       KSPConvergedReason reason;
    420:       KSPGetConvergedReason(ksp,&reason);
    421:       PetscPrintf(PETSC_COMM_WORLD,"KSPConvergedReason: %D\n", reason);
    422:     }

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

    426:   /*
    427:      Free work space.  All PETSc objects should be destroyed when they
    428:      are no longer needed.
    429:   */
    430:   MatDestroy(&A); VecDestroy(&b);
    431:   VecDestroy(&u); VecDestroy(&x);
    432:   KSPDestroy(&ksp);
    433:   PetscPreLoadEnd();
    434:   /* -----------------------------------------------------------
    435:                       End of linear solver loop
    436:      ----------------------------------------------------------- */

    438:   PetscFinalize();
    439:   return 0;
    440: }