Actual source code: ex15.c

      1: static char help[] = "Solves a linear system in parallel with KSP.  Also\n\
      2: illustrates setting a user-defined shell preconditioner and using the\n\
      3: Input parameters include:\n\
      4:   -user_defined_pc : Activate a user-defined preconditioner\n\n";

      6: /*
      7:   Include "petscksp.h" so that we can use KSP solvers.  Note that this file
      8:   automatically includes:
      9:      petscsys.h       - base PETSc routines   petscvec.h - vectors
     10:      petscmat.h - matrices
     11:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
     12:      petscviewer.h - viewers               petscpc.h  - preconditioners
     13: */
     14: #include <petscksp.h>

     16: /* Define context for user-provided preconditioner */
     17: typedef struct {
     18:   Vec diag;
     19: } SampleShellPC;

     21: /* Declare routines for user-provided preconditioner */
     22: extern PetscErrorCode SampleShellPCCreate(SampleShellPC **);
     23: extern PetscErrorCode SampleShellPCSetUp(PC, Mat, Vec);
     24: extern PetscErrorCode SampleShellPCApply(PC, Vec x, Vec y);
     25: extern PetscErrorCode SampleShellPCDestroy(PC);

     27: /*
     28:    User-defined routines.  Note that immediately before each routine below,
     29:    If defined, this macro is used in the PETSc error handlers to provide a
     30:    complete traceback of routine names.  All PETSc library routines use this
     31:    macro, and users can optionally employ it as well in their application
     32:    codes.  Note that users can get a traceback of PETSc errors regardless of
     33:    provides the added traceback detail of the application routine names.
     34: */

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

     48:   PetscFunctionBeginUser;
     49:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
     50:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
     51:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));

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

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

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

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

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

    128:   /*
    129:      Set exact solution; then compute right-hand-side vector.
    130:   */
    131:   PetscCall(VecSet(u, one));
    132:   PetscCall(MatMult(A, u, b));

    134:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    135:                 Create the linear solver and set various options
    136:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

    149:   /*
    150:      Set linear solver defaults for this problem (optional).
    151:      - By extracting the KSP and PC contexts from the KSP context,
    152:        we can then directly call any KSP and PC routines
    153:        to set various options.
    154:   */
    155:   PetscCall(KSPGetPC(ksp, &pc));
    156:   PetscCall(KSPSetTolerances(ksp, 1.e-7, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT));

    158:   /*
    159:      Set a user-defined "shell" preconditioner if desired
    160:   */
    161:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-user_defined_pc", &user_defined_pc, NULL));
    162:   if (user_defined_pc) {
    163:     /* (Required) Indicate to PETSc that we're using a "shell" preconditioner */
    164:     PetscCall(PCSetType(pc, PCSHELL));

    166:     /* (Optional) Create a context for the user-defined preconditioner; this
    167:        context can be used to contain any application-specific data. */
    168:     PetscCall(SampleShellPCCreate(&shell));

    170:     /* (Required) Set the user-defined routine for applying the preconditioner */
    171:     PetscCall(PCShellSetApply(pc, SampleShellPCApply));
    172:     PetscCall(PCShellSetContext(pc, shell));

    174:     /* (Optional) Set user-defined function to free objects used by custom preconditioner */
    175:     PetscCall(PCShellSetDestroy(pc, SampleShellPCDestroy));

    177:     /* (Optional) Set a name for the preconditioner, used for PCView() */
    178:     PetscCall(PCShellSetName(pc, "MyPreconditioner"));

    180:     /* (Optional) Do any setup required for the preconditioner */
    181:     /* Note: This function could be set with PCShellSetSetUp and it would be called when necessary */
    182:     PetscCall(SampleShellPCSetUp(pc, A, x));

    184:   } else {
    185:     PetscCall(PCSetType(pc, PCJACOBI));
    186:   }

    188:   /*
    189:     Set runtime options, e.g.,
    190:         -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
    191:     These options will override those specified above as long as
    192:     KSPSetFromOptions() is called _after_ any other customization
    193:     routines.
    194:   */
    195:   PetscCall(KSPSetFromOptions(ksp));

    197:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    198:                       Solve the linear system
    199:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

    203:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    204:                       Check solution and clean up
    205:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    207:   /*
    208:      Check the error
    209:   */
    210:   PetscCall(VecAXPY(x, none, u));
    211:   PetscCall(VecNorm(x, NORM_2, &norm));
    212:   PetscCall(KSPGetIterationNumber(ksp, &its));
    213:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g iterations %" PetscInt_FMT "\n", (double)norm, its));

    215:   /*
    216:      Free work space.  All PETSc objects should be destroyed when they
    217:      are no longer needed.
    218:   */
    219:   PetscCall(KSPDestroy(&ksp));
    220:   PetscCall(VecDestroy(&u));
    221:   PetscCall(VecDestroy(&x));
    222:   PetscCall(VecDestroy(&b));
    223:   PetscCall(MatDestroy(&A));

    225:   PetscCall(PetscFinalize());
    226:   return 0;
    227: }

    229: /***********************************************************************/
    230: /*          Routines for a user-defined shell preconditioner           */
    231: /***********************************************************************/

    233: /*
    234:    SampleShellPCCreate - This routine creates a user-defined
    235:    preconditioner context.

    237:    Output Parameter:
    238: .  shell - user-defined preconditioner context
    239: */
    240: PetscErrorCode SampleShellPCCreate(SampleShellPC **shell)
    241: {
    242:   SampleShellPC *newctx;

    244:   PetscFunctionBeginUser;
    245:   PetscCall(PetscNew(&newctx));
    246:   newctx->diag = 0;
    247:   *shell       = newctx;
    248:   PetscFunctionReturn(PETSC_SUCCESS);
    249: }
    250: /* ------------------------------------------------------------------- */
    251: /*
    252:    SampleShellPCSetUp - This routine sets up a user-defined
    253:    preconditioner context.

    255:    Input Parameters:
    256: .  pc    - preconditioner object
    257: .  pmat  - preconditioner matrix
    258: .  x     - vector

    260:    Output Parameter:
    261: .  shell - fully set up user-defined preconditioner context

    263:    Notes:
    264:    In this example, we define the shell preconditioner to be Jacobi's
    265:    method.  Thus, here we create a work vector for storing the reciprocal
    266:    of the diagonal of the preconditioner matrix; this vector is then
    267:    used within the routine SampleShellPCApply().
    268: */
    269: PetscErrorCode SampleShellPCSetUp(PC pc, Mat pmat, Vec x)
    270: {
    271:   SampleShellPC *shell;
    272:   Vec            diag;

    274:   PetscFunctionBeginUser;
    275:   PetscCall(PCShellGetContext(pc, &shell));
    276:   PetscCall(VecDuplicate(x, &diag));
    277:   PetscCall(MatGetDiagonal(pmat, diag));
    278:   PetscCall(VecReciprocal(diag));

    280:   shell->diag = diag;
    281:   PetscFunctionReturn(PETSC_SUCCESS);
    282: }
    283: /* ------------------------------------------------------------------- */
    284: /*
    285:    SampleShellPCApply - This routine demonstrates the use of a
    286:    user-provided preconditioner.

    288:    Input Parameters:
    289: +  pc - preconditioner object
    290: -  x - input vector

    292:    Output Parameter:
    293: .  y - preconditioned vector

    295:    Notes:
    296:    This code implements the Jacobi preconditioner, merely as an
    297:    example of working with a PCSHELL.  Note that the Jacobi method
    298:    is already provided within PETSc.
    299: */
    300: PetscErrorCode SampleShellPCApply(PC pc, Vec x, Vec y)
    301: {
    302:   SampleShellPC *shell;

    304:   PetscFunctionBeginUser;
    305:   PetscCall(PCShellGetContext(pc, &shell));
    306:   PetscCall(VecPointwiseMult(y, x, shell->diag));
    307:   PetscFunctionReturn(PETSC_SUCCESS);
    308: }
    309: /* ------------------------------------------------------------------- */
    310: /*
    311:    SampleShellPCDestroy - This routine destroys a user-defined
    312:    preconditioner context.

    314:    Input Parameter:
    315: .  shell - user-defined preconditioner context
    316: */
    317: PetscErrorCode SampleShellPCDestroy(PC pc)
    318: {
    319:   SampleShellPC *shell;

    321:   PetscFunctionBeginUser;
    322:   PetscCall(PCShellGetContext(pc, &shell));
    323:   PetscCall(VecDestroy(&shell->diag));
    324:   PetscCall(PetscFree(shell));
    325:   PetscFunctionReturn(PETSC_SUCCESS);
    326: }

    328: /*TEST

    330:    build:
    331:       requires: !complex !single

    333:    test:
    334:       nsize: 2
    335:       args: -ksp_view -user_defined_pc -ksp_gmres_cgs_refinement_type refine_always

    337:    test:
    338:       suffix: tsirm
    339:       args: -m 60 -n 60 -ksp_type tsirm -pc_type ksp -ksp_monitor_short -ksp_ksp_type fgmres -ksp_ksp_rtol 1e-10 -ksp_pc_type mg -ksp_ksp_max_it 30
    340:       timeoutfactor: 4

    342: TEST*/