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

      2: static char help[] = "Solves two linear systems in parallel with KSP.  The code\n\
      3: illustrates repeated solution of linear systems with the same preconditioner\n\
      4: method but different matrices (having the same nonzero structure).  The code\n\
      5: also uses multiple profiling stages.  Input arguments are\n\
      6:   -m <size> : problem size\n\
      7:   -mat_nonsym : use nonsymmetric matrix (default is symmetric)\n\n";

      9: /*T
     10:    Concepts: KSP^repeatedly solving linear systems;
     11:    Concepts: PetscLog^profiling multiple stages of code;
     12:    Processors: n
     13: T*/

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

     27: int main(int argc,char **args)
     28: {
     29:   KSP            ksp;              /* linear solver context */
     30:   Mat            C;                /* matrix */
     31:   Vec            x,u,b;            /* approx solution, RHS, exact solution */
     32:   PetscReal      norm;             /* norm of solution error */
     33:   PetscScalar    v,none = -1.0;
     34:   PetscInt       Ii,J,ldim,low,high,iglobal,Istart,Iend;
     36:   PetscInt       i,j,m = 3,n = 2,its;
     37:   PetscMPIInt    size,rank;
     38:   PetscBool      mat_nonsymmetric = PETSC_FALSE;
     39:   PetscBool      testnewC         = PETSC_FALSE;
     40: #if defined(PETSC_USE_LOG)
     41:   PetscLogStage stages[2];
     42: #endif

     44:   PetscInitialize(&argc,&args,(char*)0,help);
     45:   PetscOptionsGetInt(NULL,"-m",&m,NULL);
     46:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
     47:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
     48:   n    = 2*size;

     50:   /*
     51:      Set flag if we are doing a nonsymmetric problem; the default is symmetric.
     52:   */
     53:   PetscOptionsGetBool(NULL,"-mat_nonsym",&mat_nonsymmetric,NULL);

     55:   /*
     56:      Register two stages for separate profiling of the two linear solves.
     57:      Use the runtime option -log_summary for a printout of performance
     58:      statistics at the program's conlusion.
     59:   */
     60:   PetscLogStageRegister("Original Solve",&stages[0]);
     61:   PetscLogStageRegister("Second Solve",&stages[1]);

     63:   /* -------------- Stage 0: Solve Original System ---------------------- */
     64:   /*
     65:      Indicate to PETSc profiling that we're beginning the first stage
     66:   */
     67:   PetscLogStagePush(stages[0]);

     69:   /*
     70:      Create parallel matrix, specifying only its global dimensions.
     71:      When using MatCreate(), the matrix format can be specified at
     72:      runtime. Also, the parallel partitioning of the matrix is
     73:      determined by PETSc at runtime.
     74:   */
     75:   MatCreate(PETSC_COMM_WORLD,&C);
     76:   MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);
     77:   MatSetFromOptions(C);
     78:   MatSetUp(C);

     80:   /*
     81:      Currently, all PETSc parallel matrix formats are partitioned by
     82:      contiguous chunks of rows across the processors.  Determine which
     83:      rows of the matrix are locally owned.
     84:   */
     85:   MatGetOwnershipRange(C,&Istart,&Iend);

     87:   /*
     88:      Set matrix entries matrix in parallel.
     89:       - Each processor needs to insert only elements that it owns
     90:         locally (but any non-local elements will be sent to the
     91:         appropriate processor during matrix assembly).
     92:       - Always specify global row and columns of matrix entries.
     93:   */
     94:   for (Ii=Istart; Ii<Iend; Ii++) {
     95:     v = -1.0; i = Ii/n; j = Ii - i*n;
     96:     if (i>0)   {J = Ii - n; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
     97:     if (i<m-1) {J = Ii + n; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
     98:     if (j>0)   {J = Ii - 1; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
     99:     if (j<n-1) {J = Ii + 1; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
    100:     v = 4.0; MatSetValues(C,1,&Ii,1,&Ii,&v,ADD_VALUES);
    101:   }

    103:   /*
    104:      Make the matrix nonsymmetric if desired
    105:   */
    106:   if (mat_nonsymmetric) {
    107:     for (Ii=Istart; Ii<Iend; Ii++) {
    108:       v = -1.5; i = Ii/n;
    109:       if (i>1)   {J = Ii-n-1; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
    110:     }
    111:   } else {
    112:     MatSetOption(C,MAT_SYMMETRIC,PETSC_TRUE);
    113:     MatSetOption(C,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);
    114:   }

    116:   /*
    117:      Assemble matrix, using the 2-step process:
    118:        MatAssemblyBegin(), MatAssemblyEnd()
    119:      Computations can be done while messages are in transition
    120:      by placing code between these two statements.
    121:   */
    122:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
    123:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);

    125:   /*
    126:      Create parallel vectors.
    127:       - When using VecSetSizes(), we specify only the vector's global
    128:         dimension; the parallel partitioning is determined at runtime.
    129:       - Note: We form 1 vector from scratch and then duplicate as needed.
    130:   */
    131:   VecCreate(PETSC_COMM_WORLD,&u);
    132:   VecSetSizes(u,PETSC_DECIDE,m*n);
    133:   VecSetFromOptions(u);
    134:   VecDuplicate(u,&b);
    135:   VecDuplicate(b,&x);

    137:   /*
    138:      Currently, all parallel PETSc vectors are partitioned by
    139:      contiguous chunks across the processors.  Determine which
    140:      range of entries are locally owned.
    141:   */
    142:   VecGetOwnershipRange(x,&low,&high);

    144:   /*
    145:     Set elements within the exact solution vector in parallel.
    146:      - Each processor needs to insert only elements that it owns
    147:        locally (but any non-local entries will be sent to the
    148:        appropriate processor during vector assembly).
    149:      - Always specify global locations of vector entries.
    150:   */
    151:   VecGetLocalSize(x,&ldim);
    152:   for (i=0; i<ldim; i++) {
    153:     iglobal = i + low;
    154:     v       = (PetscScalar)(i + 100*rank);
    155:     VecSetValues(u,1,&iglobal,&v,INSERT_VALUES);
    156:   }

    158:   /*
    159:      Assemble vector, using the 2-step process:
    160:        VecAssemblyBegin(), VecAssemblyEnd()
    161:      Computations can be done while messages are in transition,
    162:      by placing code between these two statements.
    163:   */
    164:   VecAssemblyBegin(u);
    165:   VecAssemblyEnd(u);

    167:   /*
    168:      Compute right-hand-side vector
    169:   */
    170:   MatMult(C,u,b);

    172:   /*
    173:     Create linear solver context
    174:   */
    175:   KSPCreate(PETSC_COMM_WORLD,&ksp);

    177:   /*
    178:      Set operators. Here the matrix that defines the linear system
    179:      also serves as the preconditioning matrix.
    180:   */
    181:   KSPSetOperators(ksp,C,C,SAME_NONZERO_PATTERN);

    183:   /*
    184:      Set runtime options (e.g., -ksp_type <type> -pc_type <type>)
    185:   */
    186:   KSPSetFromOptions(ksp);

    188:   /*
    189:      Solve linear system.  Here we explicitly call KSPSetUp() for more
    190:      detailed performance monitoring of certain preconditioners, such
    191:      as ICC and ILU.  This call is optional, as KSPSetUp() will
    192:      automatically be called within KSPSolve() if it hasn't been
    193:      called already.
    194:   */
    195:   KSPSetUp(ksp);
    196:   KSPSolve(ksp,b,x);

    198:   /*
    199:      Check the error
    200:   */
    201:   VecAXPY(x,none,u);
    202:   VecNorm(x,NORM_2,&norm);
    203:   KSPGetIterationNumber(ksp,&its);
    204:   if (norm > 1.e-13) {
    205:     PetscPrintf(PETSC_COMM_WORLD,"Norm of error %G, Iterations %D\n",norm,its);
    206:   }

    208:   /* -------------- Stage 1: Solve Second System ---------------------- */
    209:   /*
    210:      Solve another linear system with the same method.  We reuse the KSP
    211:      context, matrix and vector data structures, and hence save the
    212:      overhead of creating new ones.

    214:      Indicate to PETSc profiling that we're concluding the first
    215:      stage with PetscLogStagePop(), and beginning the second stage with
    216:      PetscLogStagePush().
    217:   */
    218:   PetscLogStagePop();
    219:   PetscLogStagePush(stages[1]);

    221:   /*
    222:      Initialize all matrix entries to zero.  MatZeroEntries() retains the
    223:      nonzero structure of the matrix for sparse formats.
    224:   */
    225:   MatZeroEntries(C);

    227:   /*
    228:      Assemble matrix again.  Note that we retain the same matrix data
    229:      structure and the same nonzero pattern; we just change the values
    230:      of the matrix entries.
    231:   */
    232:   for (i=0; i<m; i++) {
    233:     for (j=2*rank; j<2*rank+2; j++) {
    234:       v = -1.0;  Ii = j + n*i;
    235:       if (i>0)   {J = Ii - n; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
    236:       if (i<m-1) {J = Ii + n; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
    237:       if (j>0)   {J = Ii - 1; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
    238:       if (j<n-1) {J = Ii + 1; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
    239:       v = 6.0; MatSetValues(C,1,&Ii,1,&Ii,&v,ADD_VALUES);
    240:     }
    241:   }
    242:   if (mat_nonsymmetric) {
    243:     for (Ii=Istart; Ii<Iend; Ii++) {
    244:       v = -1.5; i = Ii/n;
    245:       if (i>1)   {J = Ii-n-1; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
    246:     }
    247:   }
    248:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
    249:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);

    251:   PetscOptionsGetBool(NULL,"-test_newMat",&testnewC,NULL);
    252:   if (testnewC) {
    253:     /*
    254:      User may use a new matrix C with same nonzero pattern, e.g.
    255:       ./ex5 -ksp_monitor -mat_type sbaij -pc_type cholesky -pc_factor_mat_solver_package mumps -test_newMat
    256:     */
    257:     Mat Ctmp;
    258:     MatDuplicate(C,MAT_COPY_VALUES,&Ctmp);
    259:     MatDestroy(&C);
    260:     MatDuplicate(Ctmp,MAT_COPY_VALUES,&C);
    261:     MatDestroy(&Ctmp);
    262:   }
    263:   /*
    264:      Compute another right-hand-side vector
    265:   */
    266:   MatMult(C,u,b);

    268:   /*
    269:      Set operators. Here the matrix that defines the linear system
    270:      also serves as the preconditioning matrix.
    271:       - The flag SAME_NONZERO_PATTERN indicates that the
    272:         preconditioning matrix has identical nonzero structure
    273:         as during the last linear solve (although the values of
    274:         the entries have changed). Thus, we can save some
    275:         work in setting up the preconditioner (e.g., no need to
    276:         redo symbolic factorization for ILU/ICC preconditioners).
    277:       - If the nonzero structure of the matrix is different during
    278:         the second linear solve, then the flag DIFFERENT_NONZERO_PATTERN
    279:         must be used instead.  If you are unsure whether the
    280:         matrix structure has changed or not, use the flag
    281:         DIFFERENT_NONZERO_PATTERN.
    282:       - Caution:  If you specify SAME_NONZERO_PATTERN, PETSc
    283:         believes your assertion and does not check the structure
    284:         of the matrix.  If you erroneously claim that the structure
    285:         is the same when it actually is not, the new preconditioner
    286:         will not function correctly.  Thus, use this optimization
    287:         feature with caution!
    288:   */
    289:   KSPSetOperators(ksp,C,C,SAME_NONZERO_PATTERN);

    291:   /*
    292:      Solve linear system
    293:   */
    294:   KSPSetUp(ksp);
    295:   KSPSolve(ksp,b,x);

    297:   /*
    298:      Check the error
    299:   */
    300:   VecAXPY(x,none,u);
    301:   VecNorm(x,NORM_2,&norm);
    302:   KSPGetIterationNumber(ksp,&its);
    303:   if (norm > 1.e-4) {
    304:     PetscPrintf(PETSC_COMM_WORLD,"Norm of error %G, Iterations %D\n",norm,its);
    305:   }

    307:   /*
    308:      Free work space.  All PETSc objects should be destroyed when they
    309:      are no longer needed.
    310:   */
    311:   KSPDestroy(&ksp);
    312:   VecDestroy(&u);
    313:   VecDestroy(&x);
    314:   VecDestroy(&b);
    315:   MatDestroy(&C);

    317:   /*
    318:      Indicate to PETSc profiling that we're concluding the second stage
    319:   */
    320:   PetscLogStagePop();

    322:   PetscFinalize();
    323:   return 0;
    324: }