Actual source code: ex9.c

      2: static char help[] = "The solution of 2 different linear systems with different linear solvers.\n\
      3: Also, this example illustrates the repeated\n\
      4: solution of linear systems, while reusing matrix, vector, and solver data\n\
      5: structures throughout the process.  Note the various stages of event logging.\n\n";

      7: /*T
      8:    Concepts: KSP^repeatedly solving linear systems;
      9:    Concepts: PetscLog^profiling multiple stages of code;
     10:    Concepts: PetscLog^user-defined event profiling;
     11:    Processors: n
     12: T*/

     14: /*
     15:   Include "petscksp.h" so that we can use KSP solvers.  Note that this file
     16:   automatically includes:
     17:      petsc.h       - base PETSc routines   petscvec.h - vectors
     18:      petscsys.h    - system routines       petscmat.h - matrices
     19:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
     20:      petscviewer.h - viewers               petscpc.h  - preconditioners
     21: */
     22:  #include petscksp.h

     24: /*
     25:    Declare user-defined routines
     26: */

     32: int main(int argc,char **args)
     33: {
     34:   Vec            x1,b1,x2,b2; /* solution and RHS vectors for systems #1 and #2 */
     35:   Vec            u;              /* exact solution vector */
     36:   Mat            C1,C2;         /* matrices for systems #1 and #2 */
     37:   KSP            ksp1,ksp2;   /* KSP contexts for systems #1 and #2 */
     38:   PetscInt       ntimes = 3;     /* number of times to solve the linear systems */
     39:   PetscEvent     CHECK_ERROR;    /* event number for error checking */
     40:   PetscInt       ldim,low,high,iglobal,Istart,Iend,Istart2,Iend2;
     41:   PetscInt       I,J,i,j,m = 3,n = 2,its,t;
     43:   int            stages[3];
     44:   PetscTruth     flg;
     45:   PetscScalar    v;
     46:   PetscMPIInt    rank,size;

     48:   PetscInitialize(&argc,&args,(char *)0,help);
     49:   PetscOptionsGetInt(PETSC_NULL,"-m",&m,PETSC_NULL);
     50:   PetscOptionsGetInt(PETSC_NULL,"-t",&ntimes,PETSC_NULL);
     51:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
     52:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
     53:   n = 2*size;

     55:   /*
     56:      Register various stages for profiling
     57:   */
     58:   PetscLogStageRegister(&stages[0],"Prelim setup");
     59:   PetscLogStageRegister(&stages[1],"Linear System 1");
     60:   PetscLogStageRegister(&stages[2],"Linear System 2");

     62:   /*
     63:      Register a user-defined event for profiling (error checking).
     64:   */
     65:   CHECK_ERROR = 0;
     66:   PetscLogEventRegister(&CHECK_ERROR,"Check Error",KSP_COOKIE);

     68:   /* - - - - - - - - - - - - Stage 0: - - - - - - - - - - - - - -
     69:                         Preliminary Setup
     70:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     72:   PetscLogStagePush(stages[0]);

     74:   /*
     75:      Create data structures for first linear system.
     76:       - Create parallel matrix, specifying only its global dimensions.
     77:         When using MatCreate(), the matrix format can be specified at
     78:         runtime. Also, the parallel partitioning of the matrix is
     79:         determined by PETSc at runtime.
     80:       - Create parallel vectors.
     81:         - When using VecSetSizes(), we specify only the vector's global
     82:           dimension; the parallel partitioning is determined at runtime.
     83:         - Note: We form 1 vector from scratch and then duplicate as needed.
     84:   */
     85:   MatCreate(PETSC_COMM_WORLD,&C1);
     86:   MatSetSizes(C1,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);
     87:   MatSetFromOptions(C1);
     88:   MatGetOwnershipRange(C1,&Istart,&Iend);
     89:   VecCreate(PETSC_COMM_WORLD,&u);
     90:   VecSetSizes(u,PETSC_DECIDE,m*n);
     91:   VecSetFromOptions(u);
     92:   VecDuplicate(u,&b1);
     93:   VecDuplicate(u,&x1);

     95:   /*
     96:      Create first linear solver context.
     97:      Set runtime options (e.g., -pc_type <type>).
     98:      Note that the first linear system uses the default option
     99:      names, while the second linear systme uses a different
    100:      options prefix.
    101:   */
    102:   KSPCreate(PETSC_COMM_WORLD,&ksp1);
    103:   KSPSetFromOptions(ksp1);

    105:   /*
    106:      Set user-defined monitoring routine for first linear system.
    107:   */
    108:   PetscOptionsHasName(PETSC_NULL,"-my_ksp_monitor",&flg);
    109:   if (flg) {KSPSetMonitor(ksp1,MyKSPMonitor,PETSC_NULL,0);}

    111:   /*
    112:      Create data structures for second linear system.
    113:   */
    114:   MatCreate(PETSC_COMM_WORLD,&C2);
    115:   MatSetSizes(C2,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);
    116:   MatSetFromOptions(C2);
    117:   MatGetOwnershipRange(C2,&Istart2,&Iend2);
    118:   VecDuplicate(u,&b2);
    119:   VecDuplicate(u,&x2);

    121:   /*
    122:      Create second linear solver context
    123:   */
    124:   KSPCreate(PETSC_COMM_WORLD,&ksp2);

    126:   /*
    127:      Set different options prefix for second linear system.
    128:      Set runtime options (e.g., -s2_pc_type <type>)
    129:   */
    130:   KSPAppendOptionsPrefix(ksp2,"s2_");
    131:   KSPSetFromOptions(ksp2);

    133:   /*
    134:      Assemble exact solution vector in parallel.  Note that each
    135:      processor needs to set only its local part of the vector.
    136:   */
    137:   VecGetLocalSize(u,&ldim);
    138:   VecGetOwnershipRange(u,&low,&high);
    139:   for (i=0; i<ldim; i++) {
    140:     iglobal = i + low;
    141:     v = (PetscScalar)(i + 100*rank);
    142:     VecSetValues(u,1,&iglobal,&v,ADD_VALUES);
    143:   }
    144:   VecAssemblyBegin(u);
    145:   VecAssemblyEnd(u);

    147:   /*
    148:      Log the number of flops for computing vector entries
    149:   */
    150:   PetscLogFlops(2*ldim);

    152:   /*
    153:      End curent profiling stage
    154:   */
    155:   PetscLogStagePop();

    157:   /* --------------------------------------------------------------
    158:                         Linear solver loop:
    159:       Solve 2 different linear systems several times in succession
    160:      -------------------------------------------------------------- */

    162:   for (t=0; t<ntimes; t++) {

    164:     /* - - - - - - - - - - - - Stage 1: - - - - - - - - - - - - - -
    165:                  Assemble and solve first linear system
    166:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    168:     /*
    169:        Begin profiling stage #1
    170:     */
    171:     PetscLogStagePush(stages[1]);

    173:     /*
    174:        Initialize all matrix entries to zero.  MatZeroEntries() retains
    175:        the nonzero structure of the matrix for sparse formats.
    176:     */
    177:     if (t > 0) {MatZeroEntries(C1);}

    179:     /*
    180:        Set matrix entries in parallel.  Also, log the number of flops
    181:        for computing matrix entries.
    182:         - Each processor needs to insert only elements that it owns
    183:           locally (but any non-local elements will be sent to the
    184:           appropriate processor during matrix assembly).
    185:         - Always specify global row and columns of matrix entries.
    186:     */
    187:     for (I=Istart; I<Iend; I++) {
    188:       v = -1.0; i = I/n; j = I - i*n;
    189:       if (i>0)   {J = I - n; MatSetValues(C1,1,&I,1,&J,&v,ADD_VALUES);}
    190:       if (i<m-1) {J = I + n; MatSetValues(C1,1,&I,1,&J,&v,ADD_VALUES);}
    191:       if (j>0)   {J = I - 1; MatSetValues(C1,1,&I,1,&J,&v,ADD_VALUES);}
    192:       if (j<n-1) {J = I + 1; MatSetValues(C1,1,&I,1,&J,&v,ADD_VALUES);}
    193:       v = 4.0; MatSetValues(C1,1,&I,1,&I,&v,ADD_VALUES);
    194:     }
    195:     for (I=Istart; I<Iend; I++) { /* Make matrix nonsymmetric */
    196:       v = -1.0*(t+0.5); i = I/n;
    197:       if (i>0)   {J = I - n; MatSetValues(C1,1,&I,1,&J,&v,ADD_VALUES);}
    198:     }
    199:     PetscLogFlops(2*(Istart-Iend));

    201:     /*
    202:        Assemble matrix, using the 2-step process:
    203:          MatAssemblyBegin(), MatAssemblyEnd()
    204:        Computations can be done while messages are in transition
    205:        by placing code between these two statements.
    206:     */
    207:     MatAssemblyBegin(C1,MAT_FINAL_ASSEMBLY);
    208:     MatAssemblyEnd(C1,MAT_FINAL_ASSEMBLY);

    210:     /*
    211:        Indicate same nonzero structure of successive linear system matrices
    212:     */
    213:     MatSetOption(C1,MAT_NO_NEW_NONZERO_LOCATIONS);

    215:     /*
    216:        Compute right-hand-side vector
    217:     */
    218:     MatMult(C1,u,b1);

    220:     /*
    221:        Set operators. Here the matrix that defines the linear system
    222:        also serves as the preconditioning matrix.
    223:         - The flag SAME_NONZERO_PATTERN indicates that the
    224:           preconditioning matrix has identical nonzero structure
    225:           as during the last linear solve (although the values of
    226:           the entries have changed). Thus, we can save some
    227:           work in setting up the preconditioner (e.g., no need to
    228:           redo symbolic factorization for ILU/ICC preconditioners).
    229:         - If the nonzero structure of the matrix is different during
    230:           the second linear solve, then the flag DIFFERENT_NONZERO_PATTERN
    231:           must be used instead.  If you are unsure whether the
    232:           matrix structure has changed or not, use the flag
    233:           DIFFERENT_NONZERO_PATTERN.
    234:         - Caution:  If you specify SAME_NONZERO_PATTERN, PETSc
    235:           believes your assertion and does not check the structure
    236:           of the matrix.  If you erroneously claim that the structure
    237:           is the same when it actually is not, the new preconditioner
    238:           will not function correctly.  Thus, use this optimization
    239:           feature with caution!
    240:     */
    241:     KSPSetOperators(ksp1,C1,C1,SAME_NONZERO_PATTERN);

    243:     /*
    244:        Use the previous solution of linear system #1 as the initial
    245:        guess for the next solve of linear system #1.  The user MUST
    246:        call KSPSetInitialGuessNonzero() in indicate use of an initial
    247:        guess vector; otherwise, an initial guess of zero is used.
    248:     */
    249:     if (t>0) {
    250:       KSPSetInitialGuessNonzero(ksp1,PETSC_TRUE);
    251:     }

    253:     /*
    254:        Solve the first linear system.  Here we explicitly call
    255:        KSPSetUp() for more detailed performance monitoring of
    256:        certain preconditioners, such as ICC and ILU.  This call
    257:        is optional, ase KSPSetUp() will automatically be called
    258:        within KSPSolve() if it hasn't been called already.
    259:     */
    260:     KSPSetUp(ksp1);
    261:     KSPSolve(ksp1,b1,x1);
    262:     KSPGetIterationNumber(ksp1,&its);

    264:     /*
    265:        Check error of solution to first linear system
    266:     */
    267:     CheckError(u,x1,b1,its,CHECK_ERROR);

    269:     /* - - - - - - - - - - - - Stage 2: - - - - - - - - - - - - - -
    270:                  Assemble and solve second linear system
    271:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    273:     /*
    274:        Conclude profiling stage #1; begin profiling stage #2
    275:     */
    276:     PetscLogStagePop();
    277:     PetscLogStagePush(stages[2]);

    279:     /*
    280:        Initialize all matrix entries to zero
    281:     */
    282:     if (t > 0) {MatZeroEntries(C2);}

    284:    /*
    285:       Assemble matrix in parallel. Also, log the number of flops
    286:       for computing matrix entries.
    287:        - To illustrate the features of parallel matrix assembly, we
    288:          intentionally set the values differently from the way in
    289:          which the matrix is distributed across the processors.  Each
    290:          entry that is not owned locally will be sent to the appropriate
    291:          processor during MatAssemblyBegin() and MatAssemblyEnd().
    292:        - For best efficiency the user should strive to set as many
    293:          entries locally as possible.
    294:     */
    295:     for (i=0; i<m; i++) {
    296:       for (j=2*rank; j<2*rank+2; j++) {
    297:         v = -1.0;  I = j + n*i;
    298:         if (i>0)   {J = I - n; MatSetValues(C2,1,&I,1,&J,&v,ADD_VALUES);}
    299:         if (i<m-1) {J = I + n; MatSetValues(C2,1,&I,1,&J,&v,ADD_VALUES);}
    300:         if (j>0)   {J = I - 1; MatSetValues(C2,1,&I,1,&J,&v,ADD_VALUES);}
    301:         if (j<n-1) {J = I + 1; MatSetValues(C2,1,&I,1,&J,&v,ADD_VALUES);}
    302:         v = 6.0 + t*0.5; MatSetValues(C2,1,&I,1,&I,&v,ADD_VALUES);
    303:       }
    304:     }
    305:     for (I=Istart2; I<Iend2; I++) { /* Make matrix nonsymmetric */
    306:       v = -1.0*(t+0.5); i = I/n;
    307:       if (i>0)   {J = I - n; MatSetValues(C2,1,&I,1,&J,&v,ADD_VALUES);}
    308:     }
    309:     MatAssemblyBegin(C2,MAT_FINAL_ASSEMBLY);
    310:     MatAssemblyEnd(C2,MAT_FINAL_ASSEMBLY);
    311:     PetscLogFlops(2*(Istart-Iend));

    313:     /*
    314:        Indicate same nonzero structure of successive linear system matrices
    315:     */
    316:     MatSetOption(C2,MAT_NO_NEW_NONZERO_LOCATIONS);

    318:     /*
    319:        Compute right-hand-side vector
    320:     */
    321:     MatMult(C2,u,b2);

    323:     /*
    324:        Set operators. Here the matrix that defines the linear system
    325:        also serves as the preconditioning matrix.  Indicate same nonzero
    326:        structure of successive preconditioner matrices by setting flag
    327:        SAME_NONZERO_PATTERN.
    328:     */
    329:     KSPSetOperators(ksp2,C2,C2,SAME_NONZERO_PATTERN);

    331:     /*
    332:        Solve the second linear system
    333:     */
    334:     KSPSetUp(ksp2);
    335:     KSPSolve(ksp2,b2,x2);
    336:     KSPGetIterationNumber(ksp2,&its);

    338:     /*
    339:        Check error of solution to second linear system
    340:     */
    341:     CheckError(u,x2,b2,its,CHECK_ERROR);

    343:     /*
    344:        Conclude profiling stage #2
    345:     */
    346:     PetscLogStagePop();
    347:   }
    348:   /* --------------------------------------------------------------
    349:                        End of linear solver loop
    350:      -------------------------------------------------------------- */

    352:   /*
    353:      Free work space.  All PETSc objects should be destroyed when they
    354:      are no longer needed.
    355:   */
    356:   KSPDestroy(ksp1); KSPDestroy(ksp2);
    357:   VecDestroy(x1);   VecDestroy(x2);
    358:   VecDestroy(b1);   VecDestroy(b2);
    359:   MatDestroy(C1);   MatDestroy(C2);
    360:   VecDestroy(u);

    362:   PetscFinalize();
    363:   return 0;
    364: }
    367: /* ------------------------------------------------------------- */
    368: /*
    369:     CheckError - Checks the error of the solution.

    371:     Input Parameters:
    372:     u - exact solution
    373:     x - approximate solution
    374:     b - work vector
    375:     its - number of iterations for convergence
    376:     CHECK_ERROR - the event number for error checking
    377:                   (for use with profiling)

    379:     Notes:
    380:     In order to profile this section of code separately from the
    381:     rest of the program, we register it as an "event" with
    382:     PetscLogEventRegister() in the main program.  Then, we indicate
    383:     the start and end of this event by respectively calling
    384:         PetscLogEventBegin(CHECK_ERROR,u,x,b,0);
    385:         PetscLogEventEnd(CHECK_ERROR,u,x,b,0);
    386:     Here, we specify the objects most closely associated with
    387:     the event (the vectors u,x,b).  Such information is optional;
    388:     we could instead just use 0 instead for all objects.
    389: */
    390: PetscErrorCode CheckError(Vec u,Vec x,Vec b,PetscInt its,PetscEvent CHECK_ERROR)
    391: {
    392:   PetscScalar    none = -1.0;
    393:   PetscReal      norm;

    396:   PetscLogEventBegin(CHECK_ERROR,u,x,b,0);

    398:   /*
    399:      Compute error of the solution, using b as a work vector.
    400:   */
    401:   VecCopy(x,b);
    402:   VecAXPY(b,none,u);
    403:   VecNorm(b,NORM_2,&norm);
    404:   PetscPrintf(PETSC_COMM_WORLD,"Norm of error %A, Iterations %D\n",norm,its);
    405:   PetscLogEventEnd(CHECK_ERROR,u,x,b,0);
    406:   return 0;
    407: }
    408: /* ------------------------------------------------------------- */
    411: /*
    412:    MyKSPMonitor - This is a user-defined routine for monitoring
    413:    the KSP iterative solvers.

    415:    Input Parameters:
    416:      ksp   - iterative context
    417:      n     - iteration number
    418:      rnorm - 2-norm (preconditioned) residual value (may be estimated)
    419:      dummy - optional user-defined monitor context (unused here)
    420: */
    421: PetscErrorCode MyKSPMonitor(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
    422: {
    423:   Vec            x;

    426:   /*
    427:      Build the solution vector
    428:   */
    429:   KSPBuildSolution(ksp,PETSC_NULL,&x);

    431:   /*
    432:      Write the solution vector and residual norm to stdout.
    433:       - PetscPrintf() handles output for multiprocessor jobs
    434:         by printing from only one processor in the communicator.
    435:       - The parallel viewer PETSC_VIEWER_STDOUT_WORLD handles
    436:         data from multiple processors so that the output
    437:         is not jumbled.
    438:   */
    439:   PetscPrintf(PETSC_COMM_WORLD,"iteration %D solution vector:\n",n);
    440:   VecView(x,PETSC_VIEWER_STDOUT_WORLD);
    441:   PetscPrintf(PETSC_COMM_WORLD,"iteration %D KSP Residual norm %14.12e \n",n,rnorm);
    442:   return 0;
    443: }