Actual source code: ex23.c

      1: static char help[] = "Solves a tridiagonal linear system.\n\n";

      3: /*
      4:   Include "petscksp.h" so that we can use KSP solvers.  Note that this file
      5:   automatically includes:
      6:      petscsys.h       - base PETSc routines   petscvec.h - vectors
      7:      petscmat.h - matrices
      8:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
      9:      petscviewer.h - viewers               petscpc.h  - preconditioners

     11:   Note:  The corresponding uniprocessor example is ex1.c
     12: */
     13: #include <petscksp.h>

     15: int main(int argc, char **args)
     16: {
     17:   Vec         x, b, u;                                   /* approx solution, RHS, exact solution */
     18:   Mat         A;                                         /* linear system matrix */
     19:   KSP         ksp;                                       /* linear solver context */
     20:   PC          pc;                                        /* preconditioner context */
     21:   PetscReal   norm, tol = 1000. * PETSC_MACHINE_EPSILON; /* norm of solution error */
     22:   PetscInt    i, n      = 10, col[3], its, rstart, rend, nlocal;
     23:   PetscScalar one = 1.0, value[3];

     25:   PetscFunctionBeginUser;
     26:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
     27:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));

     29:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     30:          Compute the matrix and right-hand-side vector that define
     31:          the linear system, Ax = b.
     32:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     34:   /*
     35:      Create vectors.  Note that we form 1 vector from scratch and
     36:      then duplicate as needed. For this simple case let PETSc decide how
     37:      many elements of the vector are stored on each processor. The second
     38:      argument to VecSetSizes() below causes PETSc to decide.
     39:   */
     40:   PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
     41:   PetscCall(VecSetSizes(x, PETSC_DECIDE, n));
     42:   PetscCall(VecSetFromOptions(x));
     43:   PetscCall(VecDuplicate(x, &b));
     44:   PetscCall(VecDuplicate(x, &u));

     46:   /* Identify the starting and ending mesh points on each
     47:      processor for the interior part of the mesh. We let PETSc decide
     48:      above. */

     50:   PetscCall(VecGetOwnershipRange(x, &rstart, &rend));
     51:   PetscCall(VecGetLocalSize(x, &nlocal));

     53:   /*
     54:      Create matrix.  When using MatCreate(), the matrix format can
     55:      be specified at runtime.

     57:      Performance tuning note:  For problems of substantial size,
     58:      preallocation of matrix memory is crucial for attaining good
     59:      performance. See the matrix chapter of the users manual for details.

     61:      We pass in nlocal as the "local" size of the matrix to force it
     62:      to have the same parallel layout as the vector created above.
     63:   */
     64:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
     65:   PetscCall(MatSetSizes(A, nlocal, nlocal, n, n));
     66:   PetscCall(MatSetFromOptions(A));
     67:   PetscCall(MatSetUp(A));

     69:   /*
     70:      Assemble matrix.

     72:      The linear system is distributed across the processors by
     73:      chunks of contiguous rows, which correspond to contiguous
     74:      sections of the mesh on which the problem is discretized.
     75:      For matrix assembly, each processor contributes entries for
     76:      the part that it owns locally.
     77:   */

     79:   if (!rstart) {
     80:     rstart   = 1;
     81:     i        = 0;
     82:     col[0]   = 0;
     83:     col[1]   = 1;
     84:     value[0] = 2.0;
     85:     value[1] = -1.0;
     86:     PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
     87:   }
     88:   if (rend == n) {
     89:     rend     = n - 1;
     90:     i        = n - 1;
     91:     col[0]   = n - 2;
     92:     col[1]   = n - 1;
     93:     value[0] = -1.0;
     94:     value[1] = 2.0;
     95:     PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
     96:   }

     98:   /* Set entries corresponding to the mesh interior */
     99:   value[0] = -1.0;
    100:   value[1] = 2.0;
    101:   value[2] = -1.0;
    102:   for (i = rstart; i < rend; i++) {
    103:     col[0] = i - 1;
    104:     col[1] = i;
    105:     col[2] = i + 1;
    106:     PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
    107:   }

    109:   /* Assemble the matrix */
    110:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    111:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

    113:   /*
    114:      Set exact solution; then compute right-hand-side vector.
    115:   */
    116:   PetscCall(VecSet(u, one));
    117:   PetscCall(MatMult(A, u, b));

    119:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    120:                 Create the linear solver and set various options
    121:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    122:   /*
    123:      Create linear solver context
    124:   */
    125:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));

    127:   /*
    128:      Set operators. Here the matrix that defines the linear system
    129:      also serves as the preconditioning matrix.
    130:   */
    131:   PetscCall(KSPSetOperators(ksp, A, A));

    133:   /*
    134:      Set linear solver defaults for this problem (optional).
    135:      - By extracting the KSP and PC contexts from the KSP context,
    136:        we can then directly call any KSP and PC routines to set
    137:        various options.
    138:      - The following four statements are optional; all of these
    139:        parameters could alternatively be specified at runtime via
    140:        KSPSetFromOptions();
    141:   */
    142:   PetscCall(KSPGetPC(ksp, &pc));
    143:   PetscCall(PCSetType(pc, PCJACOBI));
    144:   PetscCall(KSPSetTolerances(ksp, 1.e-7, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT));

    146:   /*
    147:     Set runtime options, e.g.,
    148:         -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
    149:     These options will override those specified above as long as
    150:     KSPSetFromOptions() is called _after_ any other customization
    151:     routines.
    152:   */
    153:   PetscCall(KSPSetFromOptions(ksp));

    155:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    156:                       Solve the linear system
    157:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    158:   /*
    159:      Solve linear system
    160:   */
    161:   PetscCall(KSPSolve(ksp, b, x));

    163:   /*
    164:      View solver info; we could instead use the option -ksp_view to
    165:      print this info to the screen at the conclusion of KSPSolve().
    166:   */
    167:   PetscCall(KSPView(ksp, PETSC_VIEWER_STDOUT_WORLD));

    169:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    170:                       Check solution and clean up
    171:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    172:   /*
    173:      Check the error
    174:   */
    175:   PetscCall(VecAXPY(x, -1.0, u));
    176:   PetscCall(VecNorm(x, NORM_2, &norm));
    177:   PetscCall(KSPGetIterationNumber(ksp, &its));
    178:   if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its));

    180:   /*
    181:      Free work space.  All PETSc objects should be destroyed when they
    182:      are no longer needed.
    183:   */
    184:   PetscCall(VecDestroy(&x));
    185:   PetscCall(VecDestroy(&u));
    186:   PetscCall(VecDestroy(&b));
    187:   PetscCall(MatDestroy(&A));
    188:   PetscCall(KSPDestroy(&ksp));

    190:   /*
    191:      Always call PetscFinalize() before exiting a program.  This routine
    192:        - finalizes the PETSc libraries as well as MPI
    193:        - provides summary and diagnostic information if certain runtime
    194:          options are chosen (e.g., -log_view).
    195:   */
    196:   PetscCall(PetscFinalize());
    197:   return 0;
    198: }

    200: /*TEST

    202:    build:
    203:       requires: !complex !single

    205:    test:
    206:       args: -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always

    208:    test:
    209:       suffix: 2
    210:       nsize: 3
    211:       args: -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always

    213:    test:
    214:       suffix: 3
    215:       nsize: 2
    216:       args: -ksp_monitor_short -ksp_rtol 1e-6 -ksp_type pipefgmres

    218: TEST*/