Actual source code: ex10.c

      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 -mat_type aijspooles/superlu/superlu_dist/aijmumps \n\
     19:         -ksp_type preonly -pc_type cholesky -mat_type sbaijspooles/dscpack/sbaijmumps \n\
     20:         -f0 <A> -fB <B> -mat_type sbaijmumps -ksp_type preonly -pc_type cholesky -test_inertia -mat_sigma <sigma> \n\
     21:    mpirun -np <np> ex10 -f0 <datafile> -ksp_type cg -pc_type asm -pc_asm_type basic -sub_pc_type icc -mat_type sbaij
     22:  \n\n";
     23: */
     24: /*T
     25:    Concepts: KSP^solving a linear system
     26:    Processors: n
     27: T*/

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

     41: int main(int argc,char **args)
     42: {
     43:   KSP            ksp;             /* linear solver context */
     44:   Mat            A,B;            /* matrix */
     45:   Vec            x,b,u;          /* approx solution, RHS, exact solution */
     46:   PetscViewer    fd;               /* viewer */
     47:   char           file[3][PETSC_MAX_PATH_LEN];     /* input file name */
     48:   PetscTruth     table,flg,flgB=PETSC_FALSE,trans=PETSC_FALSE,partition=PETSC_FALSE;
     50:   PetscInt       its;
     51:   PetscReal      norm;
     52:   PetscLogDouble tsetup,tsetup1,tsetup2,tsolve,tsolve1,tsolve2;
     53:   PetscScalar    zero = 0.0,none = -1.0;
     54:   PetscTruth     preload=PETSC_TRUE,diagonalscale,isSymmetric,cknorm=PETSC_FALSE,Test_MatDuplicate=PETSC_FALSE;
     55:   PetscInt       num_numfac;
     56:   PetscScalar    sigma;

     58:   PetscInitialize(&argc,&args,(char *)0,help);

     60:   PetscOptionsHasName(PETSC_NULL,"-table",&table);
     61:   PetscOptionsHasName(PETSC_NULL,"-trans",&trans);
     62:   PetscOptionsHasName(PETSC_NULL,"-partition",&partition);

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

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

     93:     /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
     94:                            Load system
     95:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     97:     /*
     98:        Open binary file.  Note that we use PETSC_FILE_RDONLY to indicate
     99:        reading from this file.
    100:     */
    101:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[PreLoadIt],PETSC_FILE_RDONLY,&fd);

    103:     /*
    104:        Load the matrix and vector; then destroy the viewer.
    105:     */
    106:     MatLoad(fd,MATAIJ,&A);

    108:     PetscExceptionTry1(VecLoad(fd,PETSC_NULL,&b),PETSC_ERR_FILE_READ);
    109:     if (PetscExceptionCaught(ierr,PETSC_ERR_FILE_READ)) { /* if file contains no RHS, then use a vector of all ones */
    110:       PetscInt    m;
    111:       PetscScalar one = 1.0;
    112:       PetscLogInfo((0,"Using vector of ones for RHS\n"));
    113:       MatGetLocalSize(A,&m,PETSC_NULL);
    114:       VecCreate(PETSC_COMM_WORLD,&b);
    115:       VecSetSizes(b,m,PETSC_DECIDE);
    116:       VecSetFromOptions(b);
    117:       VecSet(b,one);
    118:     } else

    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(PETSC_NULL,"-mat_sigma",&sigma,&flg);
    134:     if(flg) {
    135:       PetscOptionsGetString(PETSC_NULL,"-fB",file[2],PETSC_MAX_PATH_LEN-1,&flgB);
    136:       if (flgB){
    137:         /* load B to get A = A + sigma*B */
    138:         PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[2],PETSC_FILE_RDONLY,&fd);
    139:         MatLoad(fd,MATAIJ,&B);
    140:         PetscViewerDestroy(fd);
    141:         MatAXPY(A,sigma,B,DIFFERENT_NONZERO_PATTERN); /* A <- sigma*B + A */
    142:       } else {
    143:         MatShift(A,sigma);
    144:       }
    145:     }

    147:     /* Check whether A is symmetric */
    148:     PetscOptionsHasName(PETSC_NULL, "-check_symmetry", &flg);
    149:     if (flg) {
    150:       Mat Atrans;
    151:       MatTranspose(A, &Atrans);
    152:       MatEqual(A, Atrans, &isSymmetric);
    153:       if (isSymmetric) {
    154:         PetscPrintf(PETSC_COMM_WORLD,"A is symmetric \n");
    155:       } else {
    156:         PetscPrintf(PETSC_COMM_WORLD,"A is non-symmetric \n");
    157:       }
    158:       MatDestroy(Atrans);
    159:     }

    161:     /*
    162:        If the loaded matrix is larger than the vector (due to being padded
    163:        to match the block size of the system), then create a new padded vector.
    164:     */
    165:     {
    166:       PetscInt    m,n,j,mvec,start,end,indx;
    167:       Vec         tmp;
    168:       PetscScalar *bold;

    170:       /* Create a new vector b by padding the old one */
    171:       MatGetLocalSize(A,&m,&n);
    172:       VecCreate(PETSC_COMM_WORLD,&tmp);
    173:       VecSetSizes(tmp,m,PETSC_DECIDE);
    174:       VecSetFromOptions(tmp);
    175:       VecGetOwnershipRange(b,&start,&end);
    176:       VecGetLocalSize(b,&mvec);
    177:       VecGetArray(b,&bold);
    178:       for (j=0; j<mvec; j++) {
    179:         indx = start+j;
    180:         VecSetValues(tmp,1,&indx,bold+j,INSERT_VALUES);
    181:       }
    182:       VecRestoreArray(b,&bold);
    183:       VecDestroy(b);
    184:       VecAssemblyBegin(tmp);
    185:       VecAssemblyEnd(tmp);
    186:       b = tmp;
    187:     }
    188:     VecDuplicate(b,&x);
    189:     VecDuplicate(b,&u);
    190:     VecSet(x,zero);

    192:     /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
    193:                       Setup solve for system
    194:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

    222:       /* need to move the vector also */
    223:       ISDestroy(is);
    224:       ISDestroy(isn);
    225:       MatDestroy(A);
    226:       A    = BB;
    227:     }
    228:
    229:     /*
    230:        Conclude profiling last stage; begin profiling next stage.
    231:     */
    232:     PreLoadStage("KSPSetUp");

    234:     /*
    235:        We also explicitly time this stage via PetscGetTime()
    236:     */
    237:     PetscGetTime(&tsetup1);

    239:     /*
    240:        Create linear solver; set operators; set runtime options.
    241:     */
    242:     KSPCreate(PETSC_COMM_WORLD,&ksp);

    244:     num_numfac = 1;
    245:     PetscOptionsGetInt(PETSC_NULL,"-num_numfac",&num_numfac,PETSC_NULL);
    246:     while ( num_numfac-- ){
    247:       /* KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN); */
    248:     KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);
    249:     KSPSetFromOptions(ksp);

    251:     /*
    252:        Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
    253:        enable more precise profiling of setting up the preconditioner.
    254:        These calls are optional, since both will be called within
    255:        KSPSolve() if they haven't been called already.
    256:     */
    257:     KSPSetUp(ksp);
    258:     KSPSetUpOnBlocks(ksp);
    259:     PetscGetTime(&tsetup2);
    260:     tsetup = tsetup2 - tsetup1;

    262:     /*
    263:       Test MatGetInertia()
    264:       Usage:
    265:       ex10 -f0 <mat_binaryfile> -ksp_type preonly -pc_type cholesky -mat_type seqsbaij -test_inertia -mat_sigma <sigma>
    266:      */
    267:     PetscOptionsHasName(PETSC_NULL,"-test_inertia",&flg);
    268:     if (flg){
    269:       PC        pc;
    270:       PetscInt  nneg, nzero, npos;
    271:       Mat       F;
    272:
    273:       KSPGetPC(ksp,&pc);
    274:       PCGetFactoredMatrix(pc,&F);
    275:       MatGetInertia(F,&nneg,&nzero,&npos);
    276:       PetscPrintf(PETSC_COMM_SELF," MatInertia: nneg: %D, nzero: %D, npos: %D\n",nneg,nzero,npos);
    277:     }

    279:     /*
    280:        Tests "diagonal-scaling of preconditioned residual norm" as used
    281:        by many ODE integrator codes including PVODE. Note this is different
    282:        than diagonally scaling the matrix before computing the preconditioner
    283:     */
    284:     PetscOptionsHasName(PETSC_NULL,"-diagonal_scale",&diagonalscale);
    285:     if (diagonalscale) {
    286:       PC       pc;
    287:       PetscInt j,start,end,n;
    288:       Vec      scale;
    289:
    290:       KSPGetPC(ksp,&pc);
    291:       VecGetSize(x,&n);
    292:       VecDuplicate(x,&scale);
    293:       VecGetOwnershipRange(scale,&start,&end);
    294:       for (j=start; j<end; j++) {
    295:         VecSetValue(scale,j,((PetscReal)(j+1))/((PetscReal)n),INSERT_VALUES);
    296:       }
    297:       VecAssemblyBegin(scale);
    298:       VecAssemblyEnd(scale);
    299:       PCDiagonalScaleSet(pc,scale);
    300:       VecDestroy(scale);

    302:     }

    304:     /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
    305:                            Solve system
    306:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    308:     /*
    309:        Begin profiling next stage
    310:     */
    311:     PreLoadStage("KSPSolve");

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

    343:    /*
    344:        Conclude profiling this stage
    345:     */
    346:     PreLoadStage("Cleanup");

    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,none,b);
    361:     VecNorm(u,NORM_2,&norm);

    363:     /*
    364:        Write output (optinally using table for solver details).
    365:         - PetscPrintf() handles output for multiprocessor jobs
    366:           by printing from only one processor in the communicator.
    367:         - KSPView() prints information about the linear solver.
    368:     */
    369:     if (table) {
    370:       char        *matrixname,kspinfo[120];
    371:       PetscViewer viewer;

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

    382:       /*
    383:          Destroy the viewer
    384:       */
    385:       PetscViewerDestroy(viewer);
    386:     } else {
    387:       PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
    388:       PetscPrintf(PETSC_COMM_WORLD,"Residual norm %A\n",norm);
    389:     }

    391:     PetscOptionsHasName(PETSC_NULL, "-ksp_reason", &flg);
    392:     if (flg){
    393:       KSPConvergedReason reason;
    394:       KSPGetConvergedReason(ksp,&reason);
    395:       PetscPrintf(PETSC_COMM_WORLD,"KSPConvergedReason: %D\n", reason);
    396:     }
    397:
    398:     } /* while ( num_numfac-- ) */

    400:     /*
    401:        Free work space.  All PETSc objects should be destroyed when they
    402:        are no longer needed.
    403:     */
    404:     MatDestroy(A); VecDestroy(b);
    405:     VecDestroy(u); VecDestroy(x);
    406:     KSPDestroy(ksp);
    407:     if (flgB) { MatDestroy(B); }
    408:   PreLoadEnd();
    409:   /* -----------------------------------------------------------
    410:                       End of linear solver loop
    411:      ----------------------------------------------------------- */

    413:   PetscFinalize();
    414:   return 0;
    415: }