Actual source code: eptorsion1.c

      1: /* Program usage: mpiexec -n 1 eptorsion1 [-help] [all TAO options] */

      3: /* ----------------------------------------------------------------------

      5:   Elastic-plastic torsion problem.

      7:   The elastic plastic torsion problem arises from the determination
      8:   of the stress field on an infinitely long cylindrical bar, which is
      9:   equivalent to the solution of the following problem:

     11:   min{ .5 * integral(||gradient(v(x))||^2 dx) - C * integral(v(x) dx)}

     13:   where C is the torsion angle per unit length.

     15:   The multiprocessor version of this code is eptorsion2.c.

     17: ---------------------------------------------------------------------- */

     19: /*
     20:   Include "petsctao.h" so that we can use TAO solvers.  Note that this
     21:   file automatically includes files for lower-level support, such as those
     22:   provided by the PETSc library:
     23:      petsc.h       - base PETSc routines   petscvec.h - vectors
     24:      petscsys.h    - system routines        petscmat.h - matrices
     25:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
     26:      petscviewer.h - viewers               petscpc.h  - preconditioners
     27: */

     29: #include <petsctao.h>

     31: static char help[] = "Demonstrates use of the TAO package to solve \n\
     32: unconstrained minimization problems on a single processor.  This example \n\
     33: is based on the Elastic-Plastic Torsion (dept) problem from the MINPACK-2 \n\
     34: test suite.\n\
     35: The command line options are:\n\
     36:   -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
     37:   -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
     38:   -par <param>, where <param> = angle of twist per unit length\n\n";

     40: /*
     41:    User-defined application context - contains data needed by the
     42:    application-provided call-back routines, FormFunction(),
     43:    FormGradient(), and FormHessian().
     44: */

     46: typedef struct {
     47:   PetscReal param;      /* nonlinearity parameter */
     48:   PetscInt  mx, my;     /* discretization in x- and y-directions */
     49:   PetscInt  ndim;       /* problem dimension */
     50:   Vec       s, y, xvec; /* work space for computing Hessian */
     51:   PetscReal hx, hy;     /* mesh spacing in x- and y-directions */
     52: } AppCtx;

     54: /* -------- User-defined Routines --------- */

     56: PetscErrorCode FormInitialGuess(AppCtx *, Vec);
     57: PetscErrorCode FormFunction(Tao, Vec, PetscReal *, void *);
     58: PetscErrorCode FormGradient(Tao, Vec, Vec, void *);
     59: PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void *);
     60: PetscErrorCode HessianProductMat(Mat, Vec, Vec);
     61: PetscErrorCode HessianProduct(void *, Vec, Vec);
     62: PetscErrorCode MatrixFreeHessian(Tao, Vec, Mat, Mat, void *);
     63: PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);

     65: int main(int argc, char **argv)
     66: {
     67:   PetscInt    mx = 10; /* discretization in x-direction */
     68:   PetscInt    my = 10; /* discretization in y-direction */
     69:   Vec         x;       /* solution, gradient vectors */
     70:   PetscBool   flg;     /* A return value when checking for use options */
     71:   Tao         tao;     /* Tao solver context */
     72:   Mat         H;       /* Hessian matrix */
     73:   AppCtx      user;    /* application context */
     74:   PetscMPIInt size;    /* number of processes */
     75:   PetscReal   one = 1.0;

     77:   PetscBool test_lmvm = PETSC_FALSE;
     78:   KSP       ksp;
     79:   PC        pc;
     80:   Mat       M;
     81:   Vec       in, out, out2;
     82:   PetscReal mult_solve_dist;

     84:   /* Initialize TAO,PETSc */
     85:   PetscFunctionBeginUser;
     86:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
     87:   PetscCallMPI(MPI_Comm_size(MPI_COMM_WORLD, &size));
     88:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Incorrect number of processors");

     90:   /* Specify default parameters for the problem, check for command-line overrides */
     91:   user.param = 5.0;
     92:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &my, &flg));
     93:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &mx, &flg));
     94:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-par", &user.param, &flg));
     95:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_lmvm", &test_lmvm, &flg));

     97:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n---- Elastic-Plastic Torsion Problem -----\n"));
     98:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "mx: %" PetscInt_FMT "     my: %" PetscInt_FMT "   \n\n", mx, my));
     99:   user.ndim = mx * my;
    100:   user.mx   = mx;
    101:   user.my   = my;
    102:   user.hx   = one / (mx + 1);
    103:   user.hy   = one / (my + 1);

    105:   /* Allocate vectors */
    106:   PetscCall(VecCreateSeq(PETSC_COMM_SELF, user.ndim, &user.y));
    107:   PetscCall(VecDuplicate(user.y, &user.s));
    108:   PetscCall(VecDuplicate(user.y, &x));

    110:   /* Create TAO solver and set desired solution method */
    111:   PetscCall(TaoCreate(PETSC_COMM_SELF, &tao));
    112:   PetscCall(TaoSetType(tao, TAOLMVM));

    114:   /* Set solution vector with an initial guess */
    115:   PetscCall(FormInitialGuess(&user, x));
    116:   PetscCall(TaoSetSolution(tao, x));

    118:   /* Set routine for function and gradient evaluation */
    119:   PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user));

    121:   /* From command line options, determine if using matrix-free hessian */
    122:   PetscCall(PetscOptionsHasName(NULL, NULL, "-my_tao_mf", &flg));
    123:   if (flg) {
    124:     PetscCall(MatCreateShell(PETSC_COMM_SELF, user.ndim, user.ndim, user.ndim, user.ndim, (void *)&user, &H));
    125:     PetscCall(MatShellSetOperation(H, MATOP_MULT, (void (*)(void))HessianProductMat));
    126:     PetscCall(MatSetOption(H, MAT_SYMMETRIC, PETSC_TRUE));

    128:     PetscCall(TaoSetHessian(tao, H, H, MatrixFreeHessian, (void *)&user));
    129:   } else {
    130:     PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, user.ndim, user.ndim, 5, NULL, &H));
    131:     PetscCall(MatSetOption(H, MAT_SYMMETRIC, PETSC_TRUE));
    132:     PetscCall(TaoSetHessian(tao, H, H, FormHessian, (void *)&user));
    133:   }

    135:   /* Test the LMVM matrix */
    136:   if (test_lmvm) {
    137:     PetscCall(PetscOptionsSetValue(NULL, "-tao_type", "bntr"));
    138:     PetscCall(PetscOptionsSetValue(NULL, "-tao_bnk_pc_type", "lmvm"));
    139:   }

    141:   /* Check for any TAO command line options */
    142:   PetscCall(TaoSetFromOptions(tao));

    144:   /* SOLVE THE APPLICATION */
    145:   PetscCall(TaoSolve(tao));

    147:   /* Test the LMVM matrix */
    148:   if (test_lmvm) {
    149:     PetscCall(TaoGetKSP(tao, &ksp));
    150:     PetscCall(KSPGetPC(ksp, &pc));
    151:     PetscCall(PCLMVMGetMatLMVM(pc, &M));
    152:     PetscCall(VecDuplicate(x, &in));
    153:     PetscCall(VecDuplicate(x, &out));
    154:     PetscCall(VecDuplicate(x, &out2));
    155:     PetscCall(VecSet(in, 5.0));
    156:     PetscCall(MatMult(M, in, out));
    157:     PetscCall(MatSolve(M, out, out2));
    158:     PetscCall(VecAXPY(out2, -1.0, in));
    159:     PetscCall(VecNorm(out2, NORM_2, &mult_solve_dist));
    160:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)tao), "error between MatMult and MatSolve: %e\n", (double)mult_solve_dist));
    161:     PetscCall(VecDestroy(&in));
    162:     PetscCall(VecDestroy(&out));
    163:     PetscCall(VecDestroy(&out2));
    164:   }

    166:   PetscCall(TaoDestroy(&tao));
    167:   PetscCall(VecDestroy(&user.s));
    168:   PetscCall(VecDestroy(&user.y));
    169:   PetscCall(VecDestroy(&x));
    170:   PetscCall(MatDestroy(&H));

    172:   PetscCall(PetscFinalize());
    173:   return 0;
    174: }

    176: /* ------------------------------------------------------------------- */
    177: /*
    178:     FormInitialGuess - Computes an initial approximation to the solution.

    180:     Input Parameters:
    181: .   user - user-defined application context
    182: .   X    - vector

    184:     Output Parameters:
    185: .   X    - vector
    186: */
    187: PetscErrorCode FormInitialGuess(AppCtx *user, Vec X)
    188: {
    189:   PetscReal hx = user->hx, hy = user->hy, temp;
    190:   PetscReal val;
    191:   PetscInt  i, j, k, nx = user->mx, ny = user->my;

    193:   /* Compute initial guess */
    194:   PetscFunctionBeginUser;
    195:   for (j = 0; j < ny; j++) {
    196:     temp = PetscMin(j + 1, ny - j) * hy;
    197:     for (i = 0; i < nx; i++) {
    198:       k   = nx * j + i;
    199:       val = PetscMin((PetscMin(i + 1, nx - i)) * hx, temp);
    200:       PetscCall(VecSetValues(X, 1, &k, &val, ADD_VALUES));
    201:     }
    202:   }
    203:   PetscCall(VecAssemblyBegin(X));
    204:   PetscCall(VecAssemblyEnd(X));
    205:   PetscFunctionReturn(PETSC_SUCCESS);
    206: }

    208: /* ------------------------------------------------------------------- */
    209: /*
    210:    FormFunctionGradient - Evaluates the function and corresponding gradient.

    212:    Input Parameters:
    213:    tao - the Tao context
    214:    X   - the input vector
    215:    ptr - optional user-defined context, as set by TaoSetFunction()

    217:    Output Parameters:
    218:    f   - the newly evaluated function
    219:    G   - the newly evaluated gradient
    220: */
    221: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ptr)
    222: {
    223:   PetscFunctionBeginUser;
    224:   PetscCall(FormFunction(tao, X, f, ptr));
    225:   PetscCall(FormGradient(tao, X, G, ptr));
    226:   PetscFunctionReturn(PETSC_SUCCESS);
    227: }

    229: /* ------------------------------------------------------------------- */
    230: /*
    231:    FormFunction - Evaluates the function, f(X).

    233:    Input Parameters:
    234: .  tao - the Tao context
    235: .  X   - the input vector
    236: .  ptr - optional user-defined context, as set by TaoSetFunction()

    238:    Output Parameters:
    239: .  f    - the newly evaluated function
    240: */
    241: PetscErrorCode FormFunction(Tao tao, Vec X, PetscReal *f, void *ptr)
    242: {
    243:   AppCtx            *user = (AppCtx *)ptr;
    244:   PetscReal          hx = user->hx, hy = user->hy, area, three = 3.0, p5 = 0.5;
    245:   PetscReal          zero = 0.0, vb, vl, vr, vt, dvdx, dvdy, flin = 0.0, fquad = 0.0;
    246:   PetscReal          v, cdiv3 = user->param / three;
    247:   const PetscScalar *x;
    248:   PetscInt           nx = user->mx, ny = user->my, i, j, k;

    250:   PetscFunctionBeginUser;
    251:   /* Get pointer to vector data */
    252:   PetscCall(VecGetArrayRead(X, &x));

    254:   /* Compute function contributions over the lower triangular elements */
    255:   for (j = -1; j < ny; j++) {
    256:     for (i = -1; i < nx; i++) {
    257:       k  = nx * j + i;
    258:       v  = zero;
    259:       vr = zero;
    260:       vt = zero;
    261:       if (i >= 0 && j >= 0) v = x[k];
    262:       if (i < nx - 1 && j > -1) vr = x[k + 1];
    263:       if (i > -1 && j < ny - 1) vt = x[k + nx];
    264:       dvdx = (vr - v) / hx;
    265:       dvdy = (vt - v) / hy;
    266:       fquad += dvdx * dvdx + dvdy * dvdy;
    267:       flin -= cdiv3 * (v + vr + vt);
    268:     }
    269:   }

    271:   /* Compute function contributions over the upper triangular elements */
    272:   for (j = 0; j <= ny; j++) {
    273:     for (i = 0; i <= nx; i++) {
    274:       k  = nx * j + i;
    275:       vb = zero;
    276:       vl = zero;
    277:       v  = zero;
    278:       if (i < nx && j > 0) vb = x[k - nx];
    279:       if (i > 0 && j < ny) vl = x[k - 1];
    280:       if (i < nx && j < ny) v = x[k];
    281:       dvdx  = (v - vl) / hx;
    282:       dvdy  = (v - vb) / hy;
    283:       fquad = fquad + dvdx * dvdx + dvdy * dvdy;
    284:       flin  = flin - cdiv3 * (vb + vl + v);
    285:     }
    286:   }

    288:   /* Restore vector */
    289:   PetscCall(VecRestoreArrayRead(X, &x));

    291:   /* Scale the function */
    292:   area = p5 * hx * hy;
    293:   *f   = area * (p5 * fquad + flin);

    295:   PetscCall(PetscLogFlops(24.0 * nx * ny));
    296:   PetscFunctionReturn(PETSC_SUCCESS);
    297: }

    299: /* ------------------------------------------------------------------- */
    300: /*
    301:     FormGradient - Evaluates the gradient, G(X).

    303:     Input Parameters:
    304: .   tao  - the Tao context
    305: .   X    - input vector
    306: .   ptr  - optional user-defined context

    308:     Output Parameters:
    309: .   G - vector containing the newly evaluated gradient
    310: */
    311: PetscErrorCode FormGradient(Tao tao, Vec X, Vec G, void *ptr)
    312: {
    313:   AppCtx            *user = (AppCtx *)ptr;
    314:   PetscReal          zero = 0.0, p5 = 0.5, three = 3.0, area, val;
    315:   PetscInt           nx = user->mx, ny = user->my, ind, i, j, k;
    316:   PetscReal          hx = user->hx, hy = user->hy;
    317:   PetscReal          vb, vl, vr, vt, dvdx, dvdy;
    318:   PetscReal          v, cdiv3 = user->param / three;
    319:   const PetscScalar *x;

    321:   PetscFunctionBeginUser;
    322:   /* Initialize gradient to zero */
    323:   PetscCall(VecSet(G, zero));

    325:   /* Get pointer to vector data */
    326:   PetscCall(VecGetArrayRead(X, &x));

    328:   /* Compute gradient contributions over the lower triangular elements */
    329:   for (j = -1; j < ny; j++) {
    330:     for (i = -1; i < nx; i++) {
    331:       k  = nx * j + i;
    332:       v  = zero;
    333:       vr = zero;
    334:       vt = zero;
    335:       if (i >= 0 && j >= 0) v = x[k];
    336:       if (i < nx - 1 && j > -1) vr = x[k + 1];
    337:       if (i > -1 && j < ny - 1) vt = x[k + nx];
    338:       dvdx = (vr - v) / hx;
    339:       dvdy = (vt - v) / hy;
    340:       if (i != -1 && j != -1) {
    341:         ind = k;
    342:         val = -dvdx / hx - dvdy / hy - cdiv3;
    343:         PetscCall(VecSetValues(G, 1, &ind, &val, ADD_VALUES));
    344:       }
    345:       if (i != nx - 1 && j != -1) {
    346:         ind = k + 1;
    347:         val = dvdx / hx - cdiv3;
    348:         PetscCall(VecSetValues(G, 1, &ind, &val, ADD_VALUES));
    349:       }
    350:       if (i != -1 && j != ny - 1) {
    351:         ind = k + nx;
    352:         val = dvdy / hy - cdiv3;
    353:         PetscCall(VecSetValues(G, 1, &ind, &val, ADD_VALUES));
    354:       }
    355:     }
    356:   }

    358:   /* Compute gradient contributions over the upper triangular elements */
    359:   for (j = 0; j <= ny; j++) {
    360:     for (i = 0; i <= nx; i++) {
    361:       k  = nx * j + i;
    362:       vb = zero;
    363:       vl = zero;
    364:       v  = zero;
    365:       if (i < nx && j > 0) vb = x[k - nx];
    366:       if (i > 0 && j < ny) vl = x[k - 1];
    367:       if (i < nx && j < ny) v = x[k];
    368:       dvdx = (v - vl) / hx;
    369:       dvdy = (v - vb) / hy;
    370:       if (i != nx && j != 0) {
    371:         ind = k - nx;
    372:         val = -dvdy / hy - cdiv3;
    373:         PetscCall(VecSetValues(G, 1, &ind, &val, ADD_VALUES));
    374:       }
    375:       if (i != 0 && j != ny) {
    376:         ind = k - 1;
    377:         val = -dvdx / hx - cdiv3;
    378:         PetscCall(VecSetValues(G, 1, &ind, &val, ADD_VALUES));
    379:       }
    380:       if (i != nx && j != ny) {
    381:         ind = k;
    382:         val = dvdx / hx + dvdy / hy - cdiv3;
    383:         PetscCall(VecSetValues(G, 1, &ind, &val, ADD_VALUES));
    384:       }
    385:     }
    386:   }
    387:   PetscCall(VecRestoreArrayRead(X, &x));

    389:   /* Assemble gradient vector */
    390:   PetscCall(VecAssemblyBegin(G));
    391:   PetscCall(VecAssemblyEnd(G));

    393:   /* Scale the gradient */
    394:   area = p5 * hx * hy;
    395:   PetscCall(VecScale(G, area));
    396:   PetscCall(PetscLogFlops(24.0 * nx * ny));
    397:   PetscFunctionReturn(PETSC_SUCCESS);
    398: }

    400: /* ------------------------------------------------------------------- */
    401: /*
    402:    FormHessian - Forms the Hessian matrix.

    404:    Input Parameters:
    405: .  tao - the Tao context
    406: .  X    - the input vector
    407: .  ptr  - optional user-defined context, as set by TaoSetHessian()

    409:    Output Parameters:
    410: .  H     - Hessian matrix
    411: .  PrecH - optionally different preconditioning Hessian

    413:    Notes:
    414:    This routine is intended simply as an example of the interface
    415:    to a Hessian evaluation routine.  Since this example compute the
    416:    Hessian a column at a time, it is not particularly efficient and
    417:    is not recommended.
    418: */
    419: PetscErrorCode FormHessian(Tao tao, Vec X, Mat H, Mat Hpre, void *ptr)
    420: {
    421:   AppCtx    *user = (AppCtx *)ptr;
    422:   PetscInt   i, j, ndim = user->ndim;
    423:   PetscReal *y, zero = 0.0, one = 1.0;
    424:   PetscBool  assembled;

    426:   PetscFunctionBeginUser;
    427:   user->xvec = X;

    429:   /* Initialize Hessian entries and work vector to zero */
    430:   PetscCall(MatAssembled(H, &assembled));
    431:   if (assembled) PetscCall(MatZeroEntries(H));

    433:   PetscCall(VecSet(user->s, zero));

    435:   /* Loop over matrix columns to compute entries of the Hessian */
    436:   for (j = 0; j < ndim; j++) {
    437:     PetscCall(VecSetValues(user->s, 1, &j, &one, INSERT_VALUES));
    438:     PetscCall(VecAssemblyBegin(user->s));
    439:     PetscCall(VecAssemblyEnd(user->s));

    441:     PetscCall(HessianProduct(ptr, user->s, user->y));

    443:     PetscCall(VecSetValues(user->s, 1, &j, &zero, INSERT_VALUES));
    444:     PetscCall(VecAssemblyBegin(user->s));
    445:     PetscCall(VecAssemblyEnd(user->s));

    447:     PetscCall(VecGetArray(user->y, &y));
    448:     for (i = 0; i < ndim; i++) {
    449:       if (y[i] != zero) PetscCall(MatSetValues(H, 1, &i, 1, &j, &y[i], ADD_VALUES));
    450:     }
    451:     PetscCall(VecRestoreArray(user->y, &y));
    452:   }
    453:   PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
    454:   PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
    455:   PetscFunctionReturn(PETSC_SUCCESS);
    456: }

    458: /* ------------------------------------------------------------------- */
    459: /*
    460:    MatrixFreeHessian - Sets a pointer for use in computing Hessian-vector
    461:    products.

    463:    Input Parameters:
    464: .  tao - the Tao context
    465: .  X    - the input vector
    466: .  ptr  - optional user-defined context, as set by TaoSetHessian()

    468:    Output Parameters:
    469: .  H     - Hessian matrix
    470: .  PrecH - optionally different preconditioning Hessian
    471: */
    472: PetscErrorCode MatrixFreeHessian(Tao tao, Vec X, Mat H, Mat PrecH, void *ptr)
    473: {
    474:   AppCtx *user = (AppCtx *)ptr;

    476:   /* Sets location of vector for use in computing matrix-vector products  of the form H(X)*y  */
    477:   PetscFunctionBeginUser;
    478:   user->xvec = X;
    479:   PetscFunctionReturn(PETSC_SUCCESS);
    480: }

    482: /* ------------------------------------------------------------------- */
    483: /*
    484:    HessianProductMat - Computes the matrix-vector product
    485:    y = mat*svec.

    487:    Input Parameters:
    488: .  mat  - input matrix
    489: .  svec - input vector

    491:    Output Parameters:
    492: .  y    - solution vector
    493: */
    494: PetscErrorCode HessianProductMat(Mat mat, Vec svec, Vec y)
    495: {
    496:   void *ptr;

    498:   PetscFunctionBeginUser;
    499:   PetscCall(MatShellGetContext(mat, &ptr));
    500:   PetscCall(HessianProduct(ptr, svec, y));
    501:   PetscFunctionReturn(PETSC_SUCCESS);
    502: }

    504: /* ------------------------------------------------------------------- */
    505: /*
    506:    Hessian Product - Computes the matrix-vector product:
    507:    y = f''(x)*svec.

    509:    Input Parameters:
    510: .  ptr  - pointer to the user-defined context
    511: .  svec - input vector

    513:    Output Parameters:
    514: .  y    - product vector
    515: */
    516: PetscErrorCode HessianProduct(void *ptr, Vec svec, Vec y)
    517: {
    518:   AppCtx            *user = (AppCtx *)ptr;
    519:   PetscReal          p5 = 0.5, zero = 0.0, one = 1.0, hx, hy, val, area;
    520:   const PetscScalar *x, *s;
    521:   PetscReal          v, vb, vl, vr, vt, hxhx, hyhy;
    522:   PetscInt           nx, ny, i, j, k, ind;

    524:   PetscFunctionBeginUser;
    525:   nx   = user->mx;
    526:   ny   = user->my;
    527:   hx   = user->hx;
    528:   hy   = user->hy;
    529:   hxhx = one / (hx * hx);
    530:   hyhy = one / (hy * hy);

    532:   /* Get pointers to vector data */
    533:   PetscCall(VecGetArrayRead(user->xvec, &x));
    534:   PetscCall(VecGetArrayRead(svec, &s));

    536:   /* Initialize product vector to zero */
    537:   PetscCall(VecSet(y, zero));

    539:   /* Compute f''(x)*s over the lower triangular elements */
    540:   for (j = -1; j < ny; j++) {
    541:     for (i = -1; i < nx; i++) {
    542:       k  = nx * j + i;
    543:       v  = zero;
    544:       vr = zero;
    545:       vt = zero;
    546:       if (i != -1 && j != -1) v = s[k];
    547:       if (i != nx - 1 && j != -1) {
    548:         vr  = s[k + 1];
    549:         ind = k + 1;
    550:         val = hxhx * (vr - v);
    551:         PetscCall(VecSetValues(y, 1, &ind, &val, ADD_VALUES));
    552:       }
    553:       if (i != -1 && j != ny - 1) {
    554:         vt  = s[k + nx];
    555:         ind = k + nx;
    556:         val = hyhy * (vt - v);
    557:         PetscCall(VecSetValues(y, 1, &ind, &val, ADD_VALUES));
    558:       }
    559:       if (i != -1 && j != -1) {
    560:         ind = k;
    561:         val = hxhx * (v - vr) + hyhy * (v - vt);
    562:         PetscCall(VecSetValues(y, 1, &ind, &val, ADD_VALUES));
    563:       }
    564:     }
    565:   }

    567:   /* Compute f''(x)*s over the upper triangular elements */
    568:   for (j = 0; j <= ny; j++) {
    569:     for (i = 0; i <= nx; i++) {
    570:       k  = nx * j + i;
    571:       v  = zero;
    572:       vl = zero;
    573:       vb = zero;
    574:       if (i != nx && j != ny) v = s[k];
    575:       if (i != nx && j != 0) {
    576:         vb  = s[k - nx];
    577:         ind = k - nx;
    578:         val = hyhy * (vb - v);
    579:         PetscCall(VecSetValues(y, 1, &ind, &val, ADD_VALUES));
    580:       }
    581:       if (i != 0 && j != ny) {
    582:         vl  = s[k - 1];
    583:         ind = k - 1;
    584:         val = hxhx * (vl - v);
    585:         PetscCall(VecSetValues(y, 1, &ind, &val, ADD_VALUES));
    586:       }
    587:       if (i != nx && j != ny) {
    588:         ind = k;
    589:         val = hxhx * (v - vl) + hyhy * (v - vb);
    590:         PetscCall(VecSetValues(y, 1, &ind, &val, ADD_VALUES));
    591:       }
    592:     }
    593:   }
    594:   /* Restore vector data */
    595:   PetscCall(VecRestoreArrayRead(svec, &s));
    596:   PetscCall(VecRestoreArrayRead(user->xvec, &x));

    598:   /* Assemble vector */
    599:   PetscCall(VecAssemblyBegin(y));
    600:   PetscCall(VecAssemblyEnd(y));

    602:   /* Scale resulting vector by area */
    603:   area = p5 * hx * hy;
    604:   PetscCall(VecScale(y, area));
    605:   PetscCall(PetscLogFlops(18.0 * nx * ny));
    606:   PetscFunctionReturn(PETSC_SUCCESS);
    607: }

    609: /*TEST

    611:    build:
    612:       requires: !complex

    614:    test:
    615:       suffix: 1
    616:       args: -tao_monitor_short -tao_type ntl -tao_gatol 1.e-4

    618:    test:
    619:       suffix: 2
    620:       args: -tao_monitor_short -tao_type ntr -tao_gatol 1.e-4

    622:    test:
    623:       suffix: 3
    624:       args: -tao_monitor_short -tao_type bntr -tao_gatol 1.e-4 -my_tao_mf -tao_test_hessian

    626:    test:
    627:      suffix: 4
    628:      args: -tao_monitor_short -tao_gatol 1e-3 -tao_type bqnls

    630:    test:
    631:      suffix: 5
    632:      args: -tao_monitor_short -tao_gatol 1e-3 -tao_type blmvm

    634:    test:
    635:      suffix: 6
    636:      args: -tao_monitor_short -tao_gatol 1e-3 -tao_type bqnktr -tao_bqnk_mat_type lmvmsr1

    638:    test:
    639:      suffix: snes
    640:      args: -snes_monitor ::ascii_info_detail -tao_type snes -snes_type newtontr -ksp_type cg -snes_atol 1.e-4 -tao_mf_hessian {{0 1}} -pc_type none

    642:    test:
    643:      suffix: snes_2
    644:      args: -snes_monitor ::ascii_info_detail -tao_type snes -snes_type newtontr -snes_atol 5.e-4 -tao_mf_hessian -pc_type none -snes_tr_fallback_type cauchy

    646:    test:
    647:      suffix: snes_3
    648:      args: -snes_monitor ::ascii_info_detail -tao_type snes -snes_type newtontr -snes_atol 5.e-4 -tao_mf_hessian -pc_type lmvm -snes_tr_fallback_type cauchy

    650: TEST*/