Actual source code: ex3.c

      1: static char help[] = "Newton methods to solve u'' + u^{2} = f in parallel.\n\
      2: This example employs a user-defined monitoring routine and optionally a user-defined\n\
      3: routine to check candidate iterates produced by line search routines.\n\
      4: The command line options include:\n\
      5:   -pre_check_iterates : activate checking of iterates\n\
      6:   -post_check_iterates : activate checking of iterates\n\
      7:   -check_tol <tol>: set tolerance for iterate checking\n\
      8:   -user_precond : activate a (trivial) user-defined preconditioner\n\n";

     10: /*
     11:    Include "petscdm.h" so that we can use data management objects (DMs)
     12:    Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
     13:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
     14:    file automatically includes:
     15:      petscsys.h    - base PETSc routines
     16:      petscvec.h    - vectors
     17:      petscmat.h    - matrices
     18:      petscis.h     - index sets
     19:      petscksp.h    - Krylov subspace methods
     20:      petscviewer.h - viewers
     21:      petscpc.h     - preconditioners
     22:      petscksp.h    - linear solvers
     23: */

     25: #include <petscdm.h>
     26: #include <petscdmda.h>
     27: #include <petscsnes.h>

     29: /*
     30:    User-defined routines.
     31: */
     32: PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *);
     33: PetscErrorCode FormFunction(SNES, Vec, Vec, void *);
     34: PetscErrorCode FormInitialGuess(Vec);
     35: PetscErrorCode Monitor(SNES, PetscInt, PetscReal, void *);
     36: PetscErrorCode PreCheck(SNESLineSearch, Vec, Vec, PetscBool *, void *);
     37: PetscErrorCode PostCheck(SNESLineSearch, Vec, Vec, Vec, PetscBool *, PetscBool *, void *);
     38: PetscErrorCode PostSetSubKSP(SNESLineSearch, Vec, Vec, Vec, PetscBool *, PetscBool *, void *);
     39: PetscErrorCode MatrixFreePreconditioner(PC, Vec, Vec);

     41: /*
     42:    User-defined application context
     43: */
     44: typedef struct {
     45:   DM          da;    /* distributed array */
     46:   Vec         F;     /* right-hand side of PDE */
     47:   PetscMPIInt rank;  /* rank of processor */
     48:   PetscMPIInt size;  /* size of communicator */
     49:   PetscReal   h;     /* mesh spacing */
     50:   PetscBool   sjerr; /* if or not to test jacobian domain error */
     51: } ApplicationCtx;

     53: /*
     54:    User-defined context for monitoring
     55: */
     56: typedef struct {
     57:   PetscViewer viewer;
     58: } MonitorCtx;

     60: /*
     61:    User-defined context for checking candidate iterates that are
     62:    determined by line search methods
     63: */
     64: typedef struct {
     65:   Vec             last_step; /* previous iterate */
     66:   PetscReal       tolerance; /* tolerance for changes between successive iterates */
     67:   ApplicationCtx *user;
     68: } StepCheckCtx;

     70: typedef struct {
     71:   PetscInt its0; /* num of previous outer KSP iterations */
     72: } SetSubKSPCtx;

     74: int main(int argc, char **argv)
     75: {
     76:   SNES           snes;       /* SNES context */
     77:   SNESLineSearch linesearch; /* SNESLineSearch context */
     78:   Mat            J;          /* Jacobian matrix */
     79:   ApplicationCtx ctx;        /* user-defined context */
     80:   Vec            x, r, U, F; /* vectors */
     81:   MonitorCtx     monP;       /* monitoring context */
     82:   StepCheckCtx   checkP;     /* step-checking context */
     83:   SetSubKSPCtx   checkP1;
     84:   PetscBool      pre_check, post_check, post_setsubksp; /* flag indicating whether we're checking candidate iterates */
     85:   PetscScalar    xp, *FF, *UU, none = -1.0;
     86:   PetscInt       its, N = 5, i, maxit, maxf, xs, xm;
     87:   PetscReal      abstol, rtol, stol, norm;
     88:   PetscBool      flg, viewinitial = PETSC_FALSE;

     90:   PetscFunctionBeginUser;
     91:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
     92:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &ctx.rank));
     93:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &ctx.size));
     94:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &N, NULL));
     95:   ctx.h     = 1.0 / (N - 1);
     96:   ctx.sjerr = PETSC_FALSE;
     97:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_jacobian_domain_error", &ctx.sjerr, NULL));
     98:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_initial", &viewinitial, NULL));

    100:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    101:      Create nonlinear solver context
    102:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    104:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));

    106:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    107:      Create vector data structures; set function evaluation routine
    108:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    110:   /*
    111:      Create distributed array (DMDA) to manage parallel grid and vectors
    112:   */
    113:   PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, N, 1, 1, NULL, &ctx.da));
    114:   PetscCall(DMSetFromOptions(ctx.da));
    115:   PetscCall(DMSetUp(ctx.da));

    117:   /*
    118:      Extract global and local vectors from DMDA; then duplicate for remaining
    119:      vectors that are the same types
    120:   */
    121:   PetscCall(DMCreateGlobalVector(ctx.da, &x));
    122:   PetscCall(VecDuplicate(x, &r));
    123:   PetscCall(VecDuplicate(x, &F));
    124:   ctx.F = F;
    125:   PetscCall(VecDuplicate(x, &U));

    127:   /*
    128:      Set function evaluation routine and vector.  Whenever the nonlinear
    129:      solver needs to compute the nonlinear function, it will call this
    130:      routine.
    131:       - Note that the final routine argument is the user-defined
    132:         context that provides application-specific data for the
    133:         function evaluation routine.
    134:   */
    135:   PetscCall(SNESSetFunction(snes, r, FormFunction, &ctx));

    137:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    138:      Create matrix data structure; set Jacobian evaluation routine
    139:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    141:   PetscCall(MatCreate(PETSC_COMM_WORLD, &J));
    142:   PetscCall(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, N, N));
    143:   PetscCall(MatSetFromOptions(J));
    144:   PetscCall(MatSeqAIJSetPreallocation(J, 3, NULL));
    145:   PetscCall(MatMPIAIJSetPreallocation(J, 3, NULL, 3, NULL));

    147:   /*
    148:      Set Jacobian matrix data structure and default Jacobian evaluation
    149:      routine.  Whenever the nonlinear solver needs to compute the
    150:      Jacobian matrix, it will call this routine.
    151:       - Note that the final routine argument is the user-defined
    152:         context that provides application-specific data for the
    153:         Jacobian evaluation routine.
    154:   */
    155:   PetscCall(SNESSetJacobian(snes, J, J, FormJacobian, &ctx));

    157:   /*
    158:      Optionally allow user-provided preconditioner
    159:    */
    160:   PetscCall(PetscOptionsHasName(NULL, NULL, "-user_precond", &flg));
    161:   if (flg) {
    162:     KSP ksp;
    163:     PC  pc;
    164:     PetscCall(SNESGetKSP(snes, &ksp));
    165:     PetscCall(KSPGetPC(ksp, &pc));
    166:     PetscCall(PCSetType(pc, PCSHELL));
    167:     PetscCall(PCShellSetApply(pc, MatrixFreePreconditioner));
    168:   }

    170:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    171:      Customize nonlinear solver; set runtime options
    172:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    174:   /*
    175:      Set an optional user-defined monitoring routine
    176:   */
    177:   PetscCall(PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, 0, 0, 0, 400, 400, &monP.viewer));
    178:   PetscCall(SNESMonitorSet(snes, Monitor, &monP, 0));

    180:   /*
    181:      Set names for some vectors to facilitate monitoring (optional)
    182:   */
    183:   PetscCall(PetscObjectSetName((PetscObject)x, "Approximate Solution"));
    184:   PetscCall(PetscObjectSetName((PetscObject)U, "Exact Solution"));

    186:   /*
    187:      Set SNES/KSP/KSP/PC runtime options, e.g.,
    188:          -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
    189:   */
    190:   PetscCall(SNESSetFromOptions(snes));

    192:   /*
    193:      Set an optional user-defined routine to check the validity of candidate
    194:      iterates that are determined by line search methods
    195:   */
    196:   PetscCall(SNESGetLineSearch(snes, &linesearch));
    197:   PetscCall(PetscOptionsHasName(NULL, NULL, "-post_check_iterates", &post_check));

    199:   if (post_check) {
    200:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Activating post step checking routine\n"));
    201:     PetscCall(SNESLineSearchSetPostCheck(linesearch, PostCheck, &checkP));
    202:     PetscCall(VecDuplicate(x, &checkP.last_step));

    204:     checkP.tolerance = 1.0;
    205:     checkP.user      = &ctx;

    207:     PetscCall(PetscOptionsGetReal(NULL, NULL, "-check_tol", &checkP.tolerance, NULL));
    208:   }

    210:   PetscCall(PetscOptionsHasName(NULL, NULL, "-post_setsubksp", &post_setsubksp));
    211:   if (post_setsubksp) {
    212:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Activating post setsubksp\n"));
    213:     PetscCall(SNESLineSearchSetPostCheck(linesearch, PostSetSubKSP, &checkP1));
    214:   }

    216:   PetscCall(PetscOptionsHasName(NULL, NULL, "-pre_check_iterates", &pre_check));
    217:   if (pre_check) {
    218:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Activating pre step checking routine\n"));
    219:     PetscCall(SNESLineSearchSetPreCheck(linesearch, PreCheck, &checkP));
    220:   }

    222:   /*
    223:      Print parameters used for convergence testing (optional) ... just
    224:      to demonstrate this routine; this information is also printed with
    225:      the option -snes_view
    226:   */
    227:   PetscCall(SNESGetTolerances(snes, &abstol, &rtol, &stol, &maxit, &maxf));
    228:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "atol=%g, rtol=%g, stol=%g, maxit=%" PetscInt_FMT ", maxf=%" PetscInt_FMT "\n", (double)abstol, (double)rtol, (double)stol, maxit, maxf));

    230:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    231:      Initialize application:
    232:      Store right-hand side of PDE and exact solution
    233:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    235:   /*
    236:      Get local grid boundaries (for 1-dimensional DMDA):
    237:        xs, xm - starting grid index, width of local grid (no ghost points)
    238:   */
    239:   PetscCall(DMDAGetCorners(ctx.da, &xs, NULL, NULL, &xm, NULL, NULL));

    241:   /*
    242:      Get pointers to vector data
    243:   */
    244:   PetscCall(DMDAVecGetArray(ctx.da, F, &FF));
    245:   PetscCall(DMDAVecGetArray(ctx.da, U, &UU));

    247:   /*
    248:      Compute local vector entries
    249:   */
    250:   xp = ctx.h * xs;
    251:   for (i = xs; i < xs + xm; i++) {
    252:     FF[i] = 6.0 * xp + PetscPowScalar(xp + 1.e-12, 6.0); /* +1.e-12 is to prevent 0^6 */
    253:     UU[i] = xp * xp * xp;
    254:     xp += ctx.h;
    255:   }

    257:   /*
    258:      Restore vectors
    259:   */
    260:   PetscCall(DMDAVecRestoreArray(ctx.da, F, &FF));
    261:   PetscCall(DMDAVecRestoreArray(ctx.da, U, &UU));
    262:   if (viewinitial) {
    263:     PetscCall(VecView(U, 0));
    264:     PetscCall(VecView(F, 0));
    265:   }

    267:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    268:      Evaluate initial guess; then solve nonlinear system
    269:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    271:   /*
    272:      Note: The user should initialize the vector, x, with the initial guess
    273:      for the nonlinear solver prior to calling SNESSolve().  In particular,
    274:      to employ an initial guess of zero, the user should explicitly set
    275:      this vector to zero by calling VecSet().
    276:   */
    277:   PetscCall(FormInitialGuess(x));
    278:   PetscCall(SNESSolve(snes, NULL, x));
    279:   PetscCall(SNESGetIterationNumber(snes, &its));
    280:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %" PetscInt_FMT "\n", its));

    282:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    283:      Check solution and clean up
    284:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    286:   /*
    287:      Check the error
    288:   */
    289:   PetscCall(VecAXPY(x, none, U));
    290:   PetscCall(VecNorm(x, NORM_2, &norm));
    291:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g Iterations %" PetscInt_FMT "\n", (double)norm, its));
    292:   if (ctx.sjerr) {
    293:     SNESType snestype;
    294:     PetscCall(SNESGetType(snes, &snestype));
    295:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "SNES Type %s\n", snestype));
    296:   }

    298:   /*
    299:      Free work space.  All PETSc objects should be destroyed when they
    300:      are no longer needed.
    301:   */
    302:   PetscCall(PetscViewerDestroy(&monP.viewer));
    303:   if (post_check) PetscCall(VecDestroy(&checkP.last_step));
    304:   PetscCall(VecDestroy(&x));
    305:   PetscCall(VecDestroy(&r));
    306:   PetscCall(VecDestroy(&U));
    307:   PetscCall(VecDestroy(&F));
    308:   PetscCall(MatDestroy(&J));
    309:   PetscCall(SNESDestroy(&snes));
    310:   PetscCall(DMDestroy(&ctx.da));
    311:   PetscCall(PetscFinalize());
    312:   return 0;
    313: }

    315: /* ------------------------------------------------------------------- */
    316: /*
    317:    FormInitialGuess - Computes initial guess.

    319:    Input/Output Parameter:
    320: .  x - the solution vector
    321: */
    322: PetscErrorCode FormInitialGuess(Vec x)
    323: {
    324:   PetscScalar pfive = .50;

    326:   PetscFunctionBeginUser;
    327:   PetscCall(VecSet(x, pfive));
    328:   PetscFunctionReturn(PETSC_SUCCESS);
    329: }

    331: /* ------------------------------------------------------------------- */
    332: /*
    333:    FormFunction - Evaluates nonlinear function, F(x).

    335:    Input Parameters:
    336: .  snes - the SNES context
    337: .  x - input vector
    338: .  ctx - optional user-defined context, as set by SNESSetFunction()

    340:    Output Parameter:
    341: .  f - function vector

    343:    Note:
    344:    The user-defined context can contain any application-specific
    345:    data needed for the function evaluation.
    346: */
    347: PetscErrorCode FormFunction(SNES snes, Vec x, Vec f, void *ctx)
    348: {
    349:   ApplicationCtx    *user = (ApplicationCtx *)ctx;
    350:   DM                 da   = user->da;
    351:   PetscScalar       *ff, d;
    352:   const PetscScalar *xx, *FF;
    353:   PetscInt           i, M, xs, xm;
    354:   Vec                xlocal;

    356:   PetscFunctionBeginUser;
    357:   PetscCall(DMGetLocalVector(da, &xlocal));
    358:   /*
    359:      Scatter ghost points to local vector, using the 2-step process
    360:         DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
    361:      By placing code between these two statements, computations can
    362:      be done while messages are in transition.
    363:   */
    364:   PetscCall(DMGlobalToLocalBegin(da, x, INSERT_VALUES, xlocal));
    365:   PetscCall(DMGlobalToLocalEnd(da, x, INSERT_VALUES, xlocal));

    367:   /*
    368:      Get pointers to vector data.
    369:        - The vector xlocal includes ghost point; the vectors x and f do
    370:          NOT include ghost points.
    371:        - Using DMDAVecGetArray() allows accessing the values using global ordering
    372:   */
    373:   PetscCall(DMDAVecGetArrayRead(da, xlocal, (void *)&xx));
    374:   PetscCall(DMDAVecGetArray(da, f, &ff));
    375:   PetscCall(DMDAVecGetArrayRead(da, user->F, (void *)&FF));

    377:   /*
    378:      Get local grid boundaries (for 1-dimensional DMDA):
    379:        xs, xm  - starting grid index, width of local grid (no ghost points)
    380:   */
    381:   PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
    382:   PetscCall(DMDAGetInfo(da, NULL, &M, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));

    384:   /*
    385:      Set function values for boundary points; define local interior grid point range:
    386:         xsi - starting interior grid index
    387:         xei - ending interior grid index
    388:   */
    389:   if (xs == 0) { /* left boundary */
    390:     ff[0] = xx[0];
    391:     xs++;
    392:     xm--;
    393:   }
    394:   if (xs + xm == M) { /* right boundary */
    395:     ff[xs + xm - 1] = xx[xs + xm - 1] - 1.0;
    396:     xm--;
    397:   }

    399:   /*
    400:      Compute function over locally owned part of the grid (interior points only)
    401:   */
    402:   d = 1.0 / (user->h * user->h);
    403:   for (i = xs; i < xs + xm; i++) ff[i] = d * (xx[i - 1] - 2.0 * xx[i] + xx[i + 1]) + xx[i] * xx[i] - FF[i];

    405:   /*
    406:      Restore vectors
    407:   */
    408:   PetscCall(DMDAVecRestoreArrayRead(da, xlocal, (void *)&xx));
    409:   PetscCall(DMDAVecRestoreArray(da, f, &ff));
    410:   PetscCall(DMDAVecRestoreArrayRead(da, user->F, (void *)&FF));
    411:   PetscCall(DMRestoreLocalVector(da, &xlocal));
    412:   PetscFunctionReturn(PETSC_SUCCESS);
    413: }

    415: /* ------------------------------------------------------------------- */
    416: /*
    417:    FormJacobian - Evaluates Jacobian matrix.

    419:    Input Parameters:
    420: .  snes - the SNES context
    421: .  x - input vector
    422: .  dummy - optional user-defined context (not used here)

    424:    Output Parameters:
    425: .  jac - Jacobian matrix
    426: .  B - optionally different matrix used to construct the preconditioner

    428: */
    429: PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat B, void *ctx)
    430: {
    431:   ApplicationCtx *user = (ApplicationCtx *)ctx;
    432:   PetscScalar    *xx, d, A[3];
    433:   PetscInt        i, j[3], M, xs, xm;
    434:   DM              da = user->da;

    436:   PetscFunctionBeginUser;
    437:   if (user->sjerr) {
    438:     PetscCall(SNESSetJacobianDomainError(snes));
    439:     PetscFunctionReturn(PETSC_SUCCESS);
    440:   }
    441:   /*
    442:      Get pointer to vector data
    443:   */
    444:   PetscCall(DMDAVecGetArrayRead(da, x, &xx));
    445:   PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));

    447:   /*
    448:     Get range of locally owned matrix
    449:   */
    450:   PetscCall(DMDAGetInfo(da, NULL, &M, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));

    452:   /*
    453:      Determine starting and ending local indices for interior grid points.
    454:      Set Jacobian entries for boundary points.
    455:   */

    457:   if (xs == 0) { /* left boundary */
    458:     i    = 0;
    459:     A[0] = 1.0;

    461:     PetscCall(MatSetValues(jac, 1, &i, 1, &i, A, INSERT_VALUES));
    462:     xs++;
    463:     xm--;
    464:   }
    465:   if (xs + xm == M) { /* right boundary */
    466:     i    = M - 1;
    467:     A[0] = 1.0;
    468:     PetscCall(MatSetValues(jac, 1, &i, 1, &i, A, INSERT_VALUES));
    469:     xm--;
    470:   }

    472:   /*
    473:      Interior grid points
    474:       - Note that in this case we set all elements for a particular
    475:         row at once.
    476:   */
    477:   d = 1.0 / (user->h * user->h);
    478:   for (i = xs; i < xs + xm; i++) {
    479:     j[0] = i - 1;
    480:     j[1] = i;
    481:     j[2] = i + 1;
    482:     A[0] = A[2] = d;
    483:     A[1]        = -2.0 * d + 2.0 * xx[i];
    484:     PetscCall(MatSetValues(jac, 1, &i, 3, j, A, INSERT_VALUES));
    485:   }

    487:   /*
    488:      Assemble matrix, using the 2-step process:
    489:        MatAssemblyBegin(), MatAssemblyEnd().
    490:      By placing code between these two statements, computations can be
    491:      done while messages are in transition.

    493:      Also, restore vector.
    494:   */

    496:   PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
    497:   PetscCall(DMDAVecRestoreArrayRead(da, x, &xx));
    498:   PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
    499:   PetscFunctionReturn(PETSC_SUCCESS);
    500: }

    502: /* ------------------------------------------------------------------- */
    503: /*
    504:    Monitor - Optional user-defined monitoring routine that views the
    505:    current iterate with an x-window plot. Set by SNESMonitorSet().

    507:    Input Parameters:
    508:    snes - the SNES context
    509:    its - iteration number
    510:    norm - 2-norm function value (may be estimated)
    511:    ctx - optional user-defined context for private data for the
    512:          monitor routine, as set by SNESMonitorSet()

    514:    Note:
    515:    See the manpage for PetscViewerDrawOpen() for useful runtime options,
    516:    such as -nox to deactivate all x-window output.
    517:  */
    518: PetscErrorCode Monitor(SNES snes, PetscInt its, PetscReal fnorm, void *ctx)
    519: {
    520:   MonitorCtx *monP = (MonitorCtx *)ctx;
    521:   Vec         x;

    523:   PetscFunctionBeginUser;
    524:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "iter = %" PetscInt_FMT ",SNES Function norm %g\n", its, (double)fnorm));
    525:   PetscCall(SNESGetSolution(snes, &x));
    526:   PetscCall(VecView(x, monP->viewer));
    527:   PetscFunctionReturn(PETSC_SUCCESS);
    528: }

    530: /* ------------------------------------------------------------------- */
    531: /*
    532:    PreCheck - Optional user-defined routine that checks the validity of
    533:    candidate steps of a line search method.  Set by SNESLineSearchSetPreCheck().

    535:    Input Parameters:
    536:    snes - the SNES context
    537:    xcurrent - current solution
    538:    y - search direction and length

    540:    Output Parameters:
    541:    y         - proposed step (search direction and length) (possibly changed)
    542:    changed_y - tells if the step has changed or not
    543:  */
    544: PetscErrorCode PreCheck(SNESLineSearch linesearch, Vec xcurrent, Vec y, PetscBool *changed_y, void *ctx)
    545: {
    546:   PetscFunctionBeginUser;
    547:   *changed_y = PETSC_FALSE;
    548:   PetscFunctionReturn(PETSC_SUCCESS);
    549: }

    551: /* ------------------------------------------------------------------- */
    552: /*
    553:    PostCheck - Optional user-defined routine that checks the validity of
    554:    candidate steps of a line search method.  Set by SNESLineSearchSetPostCheck().

    556:    Input Parameters:
    557:    snes - the SNES context
    558:    ctx  - optional user-defined context for private data for the
    559:           monitor routine, as set by SNESLineSearchSetPostCheck()
    560:    xcurrent - current solution
    561:    y - search direction and length
    562:    x    - the new candidate iterate

    564:    Output Parameters:
    565:    y    - proposed step (search direction and length) (possibly changed)
    566:    x    - current iterate (possibly modified)

    568:  */
    569: PetscErrorCode PostCheck(SNESLineSearch linesearch, Vec xcurrent, Vec y, Vec x, PetscBool *changed_y, PetscBool *changed_x, void *ctx)
    570: {
    571:   PetscInt        i, iter, xs, xm;
    572:   StepCheckCtx   *check;
    573:   ApplicationCtx *user;
    574:   PetscScalar    *xa, *xa_last, tmp;
    575:   PetscReal       rdiff;
    576:   DM              da;
    577:   SNES            snes;

    579:   PetscFunctionBeginUser;
    580:   *changed_x = PETSC_FALSE;
    581:   *changed_y = PETSC_FALSE;

    583:   PetscCall(SNESLineSearchGetSNES(linesearch, &snes));
    584:   check = (StepCheckCtx *)ctx;
    585:   user  = check->user;
    586:   PetscCall(SNESGetIterationNumber(snes, &iter));

    588:   /* iteration 1 indicates we are working on the second iteration */
    589:   if (iter > 0) {
    590:     da = user->da;
    591:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Checking candidate step at iteration %" PetscInt_FMT " with tolerance %g\n", iter, (double)check->tolerance));

    593:     /* Access local array data */
    594:     PetscCall(DMDAVecGetArray(da, check->last_step, &xa_last));
    595:     PetscCall(DMDAVecGetArray(da, x, &xa));
    596:     PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));

    598:     /*
    599:        If we fail the user-defined check for validity of the candidate iterate,
    600:        then modify the iterate as we like.  (Note that the particular modification
    601:        below is intended simply to demonstrate how to manipulate this data, not
    602:        as a meaningful or appropriate choice.)
    603:     */
    604:     for (i = xs; i < xs + xm; i++) {
    605:       if (!PetscAbsScalar(xa[i])) rdiff = 2 * check->tolerance;
    606:       else rdiff = PetscAbsScalar((xa[i] - xa_last[i]) / xa[i]);
    607:       if (rdiff > check->tolerance) {
    608:         tmp        = xa[i];
    609:         xa[i]      = .5 * (xa[i] + xa_last[i]);
    610:         *changed_x = PETSC_TRUE;
    611:         PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  Altering entry %" PetscInt_FMT ": x=%g, x_last=%g, diff=%g, x_new=%g\n", i, (double)PetscAbsScalar(tmp), (double)PetscAbsScalar(xa_last[i]), (double)rdiff, (double)PetscAbsScalar(xa[i])));
    612:       }
    613:     }
    614:     PetscCall(DMDAVecRestoreArray(da, check->last_step, &xa_last));
    615:     PetscCall(DMDAVecRestoreArray(da, x, &xa));
    616:   }
    617:   PetscCall(VecCopy(x, check->last_step));
    618:   PetscFunctionReturn(PETSC_SUCCESS);
    619: }

    621: /* ------------------------------------------------------------------- */
    622: /*
    623:    PostSetSubKSP - Optional user-defined routine that reset SubKSP options when hierarchical bjacobi PC is used
    624:    e.g,
    625:      mpiexec -n 8 ./ex3 -nox -n 10000 -ksp_type fgmres -pc_type bjacobi -pc_bjacobi_blocks 4 -sub_ksp_type gmres -sub_ksp_max_it 3 -post_setsubksp -sub_ksp_rtol 1.e-16
    626:    Set by SNESLineSearchSetPostCheck().

    628:    Input Parameters:
    629:    linesearch - the LineSearch context
    630:    xcurrent - current solution
    631:    y - search direction and length
    632:    x    - the new candidate iterate

    634:    Output Parameters:
    635:    y    - proposed step (search direction and length) (possibly changed)
    636:    x    - current iterate (possibly modified)

    638:  */
    639: PetscErrorCode PostSetSubKSP(SNESLineSearch linesearch, Vec xcurrent, Vec y, Vec x, PetscBool *changed_y, PetscBool *changed_x, void *ctx)
    640: {
    641:   SetSubKSPCtx *check;
    642:   PetscInt      iter, its, sub_its, maxit;
    643:   KSP           ksp, sub_ksp, *sub_ksps;
    644:   PC            pc;
    645:   PetscReal     ksp_ratio;
    646:   SNES          snes;

    648:   PetscFunctionBeginUser;
    649:   PetscCall(SNESLineSearchGetSNES(linesearch, &snes));
    650:   check = (SetSubKSPCtx *)ctx;
    651:   PetscCall(SNESGetIterationNumber(snes, &iter));
    652:   PetscCall(SNESGetKSP(snes, &ksp));
    653:   PetscCall(KSPGetPC(ksp, &pc));
    654:   PetscCall(PCBJacobiGetSubKSP(pc, NULL, NULL, &sub_ksps));
    655:   sub_ksp = sub_ksps[0];
    656:   PetscCall(KSPGetIterationNumber(ksp, &its));         /* outer KSP iteration number */
    657:   PetscCall(KSPGetIterationNumber(sub_ksp, &sub_its)); /* inner KSP iteration number */

    659:   if (iter) {
    660:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "    ...PostCheck snes iteration %" PetscInt_FMT ", ksp_it %" PetscInt_FMT " %" PetscInt_FMT ", subksp_it %" PetscInt_FMT "\n", iter, check->its0, its, sub_its));
    661:     ksp_ratio = (PetscReal)its / check->its0;
    662:     maxit     = (PetscInt)(ksp_ratio * sub_its + 0.5);
    663:     if (maxit < 2) maxit = 2;
    664:     PetscCall(KSPSetTolerances(sub_ksp, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT, maxit));
    665:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "    ...ksp_ratio %g, new maxit %" PetscInt_FMT "\n\n", (double)ksp_ratio, maxit));
    666:   }
    667:   check->its0 = its; /* save current outer KSP iteration number */
    668:   PetscFunctionReturn(PETSC_SUCCESS);
    669: }

    671: /* ------------------------------------------------------------------- */
    672: /*
    673:    MatrixFreePreconditioner - This routine demonstrates the use of a
    674:    user-provided preconditioner.  This code implements just the null
    675:    preconditioner, which of course is not recommended for general use.

    677:    Input Parameters:
    678: +  pc - preconditioner
    679: -  x - input vector

    681:    Output Parameter:
    682: .  y - preconditioned vector
    683: */
    684: PetscErrorCode MatrixFreePreconditioner(PC pc, Vec x, Vec y)
    685: {
    686:   PetscFunctionBeginUser;
    687:   PetscCall(VecCopy(x, y));
    688:   PetscFunctionReturn(PETSC_SUCCESS);
    689: }

    691: /*TEST

    693:    test:
    694:       args: -nox -snes_monitor_cancel -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always

    696:    test:
    697:       suffix: 2
    698:       nsize: 3
    699:       args: -nox -pc_type asm -mat_type mpiaij -snes_monitor_cancel -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always

    701:    test:
    702:       suffix: 3
    703:       nsize: 2
    704:       args: -nox -snes_monitor_cancel -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always

    706:    test:
    707:       suffix: 4
    708:       args: -nox -pre_check_iterates -post_check_iterates

    710:    test:
    711:       suffix: 5
    712:       requires: double !complex !single
    713:       nsize: 2
    714:       args: -nox -snes_test_jacobian -snes_test_jacobian_view

    716:    test:
    717:       suffix: 6
    718:       requires: double !complex !single
    719:       nsize: 4
    720:       args: -test_jacobian_domain_error -snes_converged_reason -snes_check_jacobian_domain_error 1

    722:    test:
    723:       suffix: 7
    724:       requires: double !complex !single
    725:       nsize: 4
    726:       args: -test_jacobian_domain_error -snes_converged_reason -snes_type newtontr -snes_check_jacobian_domain_error 1

    728:    test:
    729:       suffix: 8
    730:       requires: double !complex !single
    731:       nsize: 4
    732:       args: -test_jacobian_domain_error -snes_converged_reason -snes_type vinewtonrsls -snes_check_jacobian_domain_error 1

    734:    test:
    735:       suffix: 9
    736:       requires: double !complex !single
    737:       nsize: 4
    738:       args: -test_jacobian_domain_error -snes_converged_reason -snes_type vinewtonssls -snes_check_jacobian_domain_error 1

    740:    test:
    741:       suffix: 10
    742:       requires: double !complex !single
    743:       nsize: 4
    744:       args: -test_jacobian_domain_error -snes_converged_reason -snes_type qn -snes_qn_scale_type jacobian -snes_check_jacobian_domain_error 1

    746:    test:
    747:       suffix: 11
    748:       requires: double !complex !single
    749:       nsize: 4
    750:       args: -test_jacobian_domain_error -snes_converged_reason -snes_type ms -snes_ms_type m62 -snes_ms_damping 0.9 -snes_check_jacobian_domain_error 1

    752:    test:
    753:       suffix: 12
    754:       args: -view_initial
    755:       filter: grep -v "type:"

    757:    test:
    758:       suffix: 13
    759:       requires: double !complex !single
    760:       nsize: 4
    761:       args: -test_jacobian_domain_error -snes_converged_reason -snes_type newtontrdc -snes_check_jacobian_domain_error 1

    763: TEST*/