Actual source code: ex12.c

      1: static char help[] = "Solves a linear system in parallel with KSP,\n\
      2: demonstrating how to register a new preconditioner (PC) type.\n\
      3: Input parameters include:\n\
      4:   -m <mesh_x>       : number of mesh points in x-direction\n\
      5:   -n <mesh_y>       : number of mesh points in y-direction\n\n";

      7: /*
      8:    Demonstrates registering a new preconditioner (PC) type.

     10:    To register a PC type whose code is linked into the executable,
     11:    use PCRegister(). To register a PC type in a dynamic library use PCRegister()

     13:    Also provide the prototype for your PCCreate_XXX() function. In
     14:    this example we use the PETSc implementation of the Jacobi method,
     15:    PCCreate_Jacobi() just as an example.

     17:    See the file src/ksp/pc/impls/jacobi/jacobi.c for details on how to
     18:    write a new PC component.

     20:    See the manual page PCRegister() for details on how to register a method.
     21: */

     23: /*
     24:   Include "petscksp.h" so that we can use KSP solvers.  Note that this file
     25:   automatically includes:
     26:      petscsys.h       - base PETSc routines   petscvec.h - vectors
     27:      petscmat.h - matrices
     28:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
     29:      petscviewer.h - viewers               petscpc.h  - preconditioners
     30: */
     31: #include <petscksp.h>

     33: PETSC_EXTERN PetscErrorCode PCCreate_Jacobi(PC);

     35: int main(int argc, char **args)
     36: {
     37:   Vec         x, b, u; /* approx solution, RHS, exact solution */
     38:   Mat         A;       /* linear system matrix */
     39:   KSP         ksp;     /* linear solver context */
     40:   PetscReal   norm;    /* norm of solution error */
     41:   PetscInt    i, j, Ii, J, Istart, Iend, m = 8, n = 7, its;
     42:   PetscScalar v, one = 1.0;
     43:   PC          pc; /* preconditioner context */

     45:   PetscFunctionBeginUser;
     46:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
     47:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
     48:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));

     50:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     51:          Compute the matrix and right-hand-side vector that define
     52:          the linear system, Ax = b.
     53:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     54:   /*
     55:      Create parallel matrix, specifying only its global dimensions.
     56:      When using MatCreate(), the matrix format can be specified at
     57:      runtime. Also, the parallel partitioning of the matrix can be
     58:      determined by PETSc at runtime.
     59:   */
     60:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
     61:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m * n, m * n));
     62:   PetscCall(MatSetFromOptions(A));
     63:   PetscCall(MatSetUp(A));

     65:   /*
     66:      Currently, all PETSc parallel matrix formats are partitioned by
     67:      contiguous chunks of rows across the processors.  Determine which
     68:      rows of the matrix are locally owned.
     69:   */
     70:   PetscCall(MatGetOwnershipRange(A, &Istart, &Iend));

     72:   /*
     73:      Set matrix elements for the 2-D, five-point stencil in parallel.
     74:       - Each processor needs to insert only elements that it owns
     75:         locally (but any non-local elements will be sent to the
     76:         appropriate processor during matrix assembly).
     77:       - Always specify global rows and columns of matrix entries.
     78:    */
     79:   for (Ii = Istart; Ii < Iend; Ii++) {
     80:     v = -1.0;
     81:     i = Ii / n;
     82:     j = Ii - i * n;
     83:     if (i > 0) {
     84:       J = Ii - n;
     85:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
     86:     }
     87:     if (i < m - 1) {
     88:       J = Ii + n;
     89:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
     90:     }
     91:     if (j > 0) {
     92:       J = Ii - 1;
     93:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
     94:     }
     95:     if (j < n - 1) {
     96:       J = Ii + 1;
     97:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
     98:     }
     99:     v = 4.0;
    100:     PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &v, INSERT_VALUES));
    101:   }

    103:   /*
    104:      Assemble matrix, using the 2-step process:
    105:        MatAssemblyBegin(), MatAssemblyEnd()
    106:      Computations can be done while messages are in transition
    107:      by placing code between these two statements.
    108:   */
    109:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    110:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

    112:   /*
    113:      Create parallel vectors.
    114:       - When using VecCreate(), VecSetSizes() and VecSetFromOptions(),
    115:       we specify only the vector's global
    116:         dimension; the parallel partitioning is determined at runtime.
    117:       - When solving a linear system, the vectors and matrices MUST
    118:         be partitioned accordingly.  PETSc automatically generates
    119:         appropriately partitioned matrices and vectors when MatCreate()
    120:         and VecCreate() are used with the same communicator.
    121:       - Note: We form 1 vector from scratch and then duplicate as needed.
    122:   */
    123:   PetscCall(VecCreate(PETSC_COMM_WORLD, &u));
    124:   PetscCall(VecSetSizes(u, PETSC_DECIDE, m * n));
    125:   PetscCall(VecSetFromOptions(u));
    126:   PetscCall(VecDuplicate(u, &b));
    127:   PetscCall(VecDuplicate(b, &x));

    129:   /*
    130:      Set exact solution; then compute right-hand-side vector.
    131:      Use an exact solution of a vector with all elements of 1.0;
    132:   */
    133:   PetscCall(VecSet(u, one));
    134:   PetscCall(MatMult(A, u, b));

    136:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    137:                 Create the linear solver and set various options
    138:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    140:   /*
    141:      Create linear solver context
    142:   */
    143:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));

    145:   /*
    146:      Set operators. Here the matrix that defines the linear system
    147:      also serves as the preconditioning matrix.
    148:   */
    149:   PetscCall(KSPSetOperators(ksp, A, A));

    151:   /*
    152:        First register a new PC type with the command PCRegister()
    153:   */
    154:   PetscCall(PCRegister("ourjacobi", PCCreate_Jacobi));

    156:   /*
    157:      Set the PC type to be the new method
    158:   */
    159:   PetscCall(KSPGetPC(ksp, &pc));
    160:   PetscCall(PCSetType(pc, "ourjacobi"));

    162:   /*
    163:     Set runtime options, e.g.,
    164:         -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
    165:     These options will override those specified above as long as
    166:     KSPSetFromOptions() is called _after_ any other customization
    167:     routines.
    168:   */
    169:   PetscCall(KSPSetFromOptions(ksp));

    171:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    172:                       Solve the linear system
    173:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    175:   PetscCall(KSPSolve(ksp, b, x));

    177:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    178:                       Check solution and clean up
    179:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    181:   /*
    182:      Check the error
    183:   */
    184:   PetscCall(VecAXPY(x, -1.0, u));
    185:   PetscCall(VecNorm(x, NORM_2, &norm));
    186:   PetscCall(KSPGetIterationNumber(ksp, &its));

    188:   /*
    189:      Print convergence information.  PetscPrintf() produces a single
    190:      print statement from all processes that share a communicator.
    191:   */
    192:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g iterations %" PetscInt_FMT "\n", (double)norm, its));

    194:   /*
    195:      Free work space.  All PETSc objects should be destroyed when they
    196:      are no longer needed.
    197:   */
    198:   PetscCall(KSPDestroy(&ksp));
    199:   PetscCall(VecDestroy(&u));
    200:   PetscCall(VecDestroy(&x));
    201:   PetscCall(VecDestroy(&b));
    202:   PetscCall(MatDestroy(&A));

    204:   /*
    205:      Always call PetscFinalize() before exiting a program.  This routine
    206:        - finalizes the PETSc libraries as well as MPI
    207:        - provides summary and diagnostic information if certain runtime
    208:          options are chosen (e.g., -log_view).
    209:   */
    210:   PetscCall(PetscFinalize());
    211:   return 0;
    212: }

    214: /*TEST

    216:    test:
    217:       args: -ksp_gmres_cgs_refinement_type refine_always

    219: TEST*/