Actual source code: ex5.c

      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:      petsc.h       - base PETSc routines   petscvec.h - vectors
     19:      petscsys.h    - system routines       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       I,J,ldim,low,high,iglobal,Istart,Iend;
     36:   PetscInt       i,j,m = 3,n = 2,its;
     37:   PetscMPIInt    size,rank;
     38:   PetscTruth     mat_nonsymmetric;
     39:   int            stages[2];

     41:   PetscInitialize(&argc,&args,(char *)0,help);
     42:   PetscOptionsGetInt(PETSC_NULL,"-m",&m,PETSC_NULL);
     43:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
     44:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
     45:   n = 2*size;

     47:   /*
     48:      Set flag if we are doing a nonsymmetric problem; the default is symmetric.
     49:   */
     50:   PetscOptionsHasName(PETSC_NULL,"-mat_nonsym",&mat_nonsymmetric);

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

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

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

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

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

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

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

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

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

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

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

    163:   /*
    164:      Compute right-hand-side vector
    165:   */
    166:   MatMult(C,u,b);
    167:
    168:   /*
    169:     Create linear solver context
    170:   */
    171:   KSPCreate(PETSC_COMM_WORLD,&ksp);

    173:   /*
    174:      Set operators. Here the matrix that defines the linear system
    175:      also serves as the preconditioning matrix.
    176:   */
    177:   KSPSetOperators(ksp,C,C,DIFFERENT_NONZERO_PATTERN);

    179:   /*
    180:      Set runtime options (e.g., -ksp_type <type> -pc_type <type>)
    181:   */

    183:   KSPSetFromOptions(ksp);

    185:   /*
    186:      Solve linear system.  Here we explicitly call KSPSetUp() for more
    187:      detailed performance monitoring of certain preconditioners, such
    188:      as ICC and ILU.  This call is optional, as KSPSetUp() will
    189:      automatically be called within KSPSolve() if it hasn't been
    190:      called already.
    191:   */
    192:   KSPSetUp(ksp);
    193:   KSPSolve(ksp,b,x);
    194:
    195:   /*
    196:      Check the error
    197:   */
    198:   VecAXPY(x,none,u);
    199:   VecNorm(x,NORM_2,&norm);
    200:   KSPGetIterationNumber(ksp,&its);
    201:   PetscPrintf(PETSC_COMM_WORLD,"Norm of error %A, Iterations %D\n",norm,its);

    203:   /* -------------- Stage 1: Solve Second System ---------------------- */
    204:   /*
    205:      Solve another linear system with the same method.  We reuse the KSP
    206:      context, matrix and vector data structures, and hence save the
    207:      overhead of creating new ones.

    209:      Indicate to PETSc profiling that we're concluding the first
    210:      stage with PetscLogStagePop(), and beginning the second stage with
    211:      PetscLogStagePush().
    212:   */
    213:   PetscLogStagePop();
    214:   PetscLogStagePush(stages[1]);

    216:   /*
    217:      Initialize all matrix entries to zero.  MatZeroEntries() retains the
    218:      nonzero structure of the matrix for sparse formats.
    219:   */
    220:   MatZeroEntries(C);

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

    246:   /*
    247:      Compute another right-hand-side vector
    248:   */
    249:   MatMult(C,u,b);

    251:   /*
    252:      Set operators. Here the matrix that defines the linear system
    253:      also serves as the preconditioning matrix.
    254:       - The flag SAME_NONZERO_PATTERN indicates that the
    255:         preconditioning matrix has identical nonzero structure
    256:         as during the last linear solve (although the values of
    257:         the entries have changed). Thus, we can save some
    258:         work in setting up the preconditioner (e.g., no need to
    259:         redo symbolic factorization for ILU/ICC preconditioners).
    260:       - If the nonzero structure of the matrix is different during
    261:         the second linear solve, then the flag DIFFERENT_NONZERO_PATTERN
    262:         must be used instead.  If you are unsure whether the
    263:         matrix structure has changed or not, use the flag
    264:         DIFFERENT_NONZERO_PATTERN.
    265:       - Caution:  If you specify SAME_NONZERO_PATTERN, PETSc
    266:         believes your assertion and does not check the structure
    267:         of the matrix.  If you erroneously claim that the structure
    268:         is the same when it actually is not, the new preconditioner
    269:         will not function correctly.  Thus, use this optimization
    270:         feature with caution!
    271:   */
    272:   KSPSetOperators(ksp,C,C,SAME_NONZERO_PATTERN);

    274:   /*
    275:      Solve linear system
    276:   */
    277:   KSPSetUp(ksp);
    278:   KSPSolve(ksp,b,x);

    280:   /*
    281:      Check the error
    282:   */
    283:   VecAXPY(x,none,u);
    284:   VecNorm(x,NORM_2,&norm);
    285:   KSPGetIterationNumber(ksp,&its);
    286:   PetscPrintf(PETSC_COMM_WORLD,"Norm of error %A, Iterations %D\n",norm,its);

    288:   /*
    289:      Free work space.  All PETSc objects should be destroyed when they
    290:      are no longer needed.
    291:   */
    292:   KSPDestroy(ksp);
    293:   VecDestroy(u);
    294:   VecDestroy(x);
    295:   VecDestroy(b);
    296:   MatDestroy(C);

    298:   /*
    299:      Indicate to PETSc profiling that we're concluding the second stage
    300:   */
    301:   PetscLogStagePop();

    303:   PetscFinalize();
    304:   return 0;
    305: }