Actual source code: eptorsion3.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;
     83:   Vec       ub, lb;
     84:   PetscInt  i = 3;

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

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

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

    107:   /* Allocate vectors */
    108:   PetscCall(VecCreateSeq(PETSC_COMM_SELF, user.ndim, &user.y));
    109:   PetscCall(VecDuplicate(user.y, &user.s));
    110:   PetscCall(VecDuplicate(user.y, &x));
    111:   PetscCall(VecDuplicate(user.y, &lb));
    112:   PetscCall(VecDuplicate(user.y, &ub));

    114:   PetscCall(VecSet(ub, 0.1));
    115:   PetscCall(VecSet(lb, -0.1));
    116:   PetscCall(VecSetValue(ub, i, 0., INSERT_VALUES));
    117:   PetscCall(VecSetValue(lb, i, 0., INSERT_VALUES));
    118:   PetscCall(VecAssemblyBegin(ub));
    119:   PetscCall(VecAssemblyBegin(lb));
    120:   PetscCall(VecAssemblyEnd(ub));
    121:   PetscCall(VecAssemblyEnd(lb));

    123:   /* Create TAO solver and set desired solution method */
    124:   PetscCall(TaoCreate(PETSC_COMM_SELF, &tao));
    125:   PetscCall(TaoSetType(tao, TAOLMVM));

    127:   /* Set solution vector with an initial guess */
    128:   PetscCall(FormInitialGuess(&user, x));
    129:   PetscCall(TaoSetSolution(tao, x));
    130:   PetscCall(TaoSetVariableBounds(tao, lb, ub));

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

    135:   /* From command line options, determine if using matrix-free hessian */
    136:   PetscCall(PetscOptionsHasName(NULL, NULL, "-my_tao_mf", &flg));
    137:   if (flg) {
    138:     PetscCall(MatCreateShell(PETSC_COMM_SELF, user.ndim, user.ndim, user.ndim, user.ndim, (void *)&user, &H));
    139:     PetscCall(MatShellSetOperation(H, MATOP_MULT, (void (*)(void))HessianProductMat));
    140:     PetscCall(MatSetOption(H, MAT_SYMMETRIC, PETSC_TRUE));

    142:     PetscCall(TaoSetHessian(tao, H, H, MatrixFreeHessian, (void *)&user));
    143:   } else {
    144:     PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, user.ndim, user.ndim, 5, NULL, &H));
    145:     PetscCall(MatSetOption(H, MAT_SYMMETRIC, PETSC_TRUE));
    146:     PetscCall(TaoSetHessian(tao, H, H, FormHessian, (void *)&user));
    147:   }

    149:   /* Test the LMVM matrix */
    150:   if (test_lmvm) {
    151:     PetscCall(PetscOptionsSetValue(NULL, "-tao_type", "bntr"));
    152:     PetscCall(PetscOptionsSetValue(NULL, "-tao_bnk_pc_type", "lmvm"));
    153:   }

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

    158:   /* SOLVE THE APPLICATION */
    159:   PetscCall(TaoSolve(tao));

    161:   /* Test the LMVM matrix */
    162:   if (test_lmvm) {
    163:     PetscCall(TaoGetKSP(tao, &ksp));
    164:     PetscCall(KSPGetPC(ksp, &pc));
    165:     PetscCall(PCLMVMGetMatLMVM(pc, &M));
    166:     PetscCall(VecDuplicate(x, &in));
    167:     PetscCall(VecDuplicate(x, &out));
    168:     PetscCall(VecDuplicate(x, &out2));
    169:     PetscCall(VecSet(in, 5.0));
    170:     PetscCall(MatMult(M, in, out));
    171:     PetscCall(MatSolve(M, out, out2));
    172:     PetscCall(VecAXPY(out2, -1.0, in));
    173:     PetscCall(VecNorm(out2, NORM_2, &mult_solve_dist));
    174:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)tao), "error between MatMult and MatSolve: %e\n", (double)mult_solve_dist));
    175:     PetscCall(VecDestroy(&in));
    176:     PetscCall(VecDestroy(&out));
    177:     PetscCall(VecDestroy(&out2));
    178:   }

    180:   PetscCall(TaoDestroy(&tao));
    181:   PetscCall(VecDestroy(&user.s));
    182:   PetscCall(VecDestroy(&user.y));
    183:   PetscCall(VecDestroy(&x));
    184:   PetscCall(MatDestroy(&H));
    185:   PetscCall(VecDestroy(&lb));
    186:   PetscCall(VecDestroy(&ub));

    188:   PetscCall(PetscFinalize());
    189:   return 0;
    190: }

    192: /* ------------------------------------------------------------------- */
    193: /*
    194:     FormInitialGuess - Computes an initial approximation to the solution.

    196:     Input Parameters:
    197: .   user - user-defined application context
    198: .   X    - vector

    200:     Output Parameters:
    201: .   X    - vector
    202: */
    203: PetscErrorCode FormInitialGuess(AppCtx *user, Vec X)
    204: {
    205:   PetscReal hx = user->hx, hy = user->hy, temp;
    206:   PetscReal val;
    207:   PetscInt  i, j, k, nx = user->mx, ny = user->my;

    209:   /* Compute initial guess */
    210:   PetscFunctionBeginUser;
    211:   for (j = 0; j < ny; j++) {
    212:     temp = PetscMin(j + 1, ny - j) * hy;
    213:     for (i = 0; i < nx; i++) {
    214:       k   = nx * j + i;
    215:       val = PetscMin((PetscMin(i + 1, nx - i)) * hx, temp);
    216:       PetscCall(VecSetValues(X, 1, &k, &val, ADD_VALUES));
    217:     }
    218:   }
    219:   PetscCall(VecAssemblyBegin(X));
    220:   PetscCall(VecAssemblyEnd(X));
    221:   PetscFunctionReturn(PETSC_SUCCESS);
    222: }

    224: /* ------------------------------------------------------------------- */
    225: /*
    226:    FormFunctionGradient - Evaluates the function and corresponding gradient.

    228:    Input Parameters:
    229:    tao - the Tao context
    230:    X   - the input vector
    231:    ptr - optional user-defined context, as set by TaoSetFunction()

    233:    Output Parameters:
    234:    f   - the newly evaluated function
    235:    G   - the newly evaluated gradient
    236: */
    237: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ptr)
    238: {
    239:   PetscFunctionBeginUser;
    240:   PetscCall(FormFunction(tao, X, f, ptr));
    241:   PetscCall(FormGradient(tao, X, G, ptr));
    242:   PetscFunctionReturn(PETSC_SUCCESS);
    243: }

    245: /* ------------------------------------------------------------------- */
    246: /*
    247:    FormFunction - Evaluates the function, f(X).

    249:    Input Parameters:
    250: .  tao - the Tao context
    251: .  X   - the input vector
    252: .  ptr - optional user-defined context, as set by TaoSetFunction()

    254:    Output Parameters:
    255: .  f    - the newly evaluated function
    256: */
    257: PetscErrorCode FormFunction(Tao tao, Vec X, PetscReal *f, void *ptr)
    258: {
    259:   AppCtx            *user = (AppCtx *)ptr;
    260:   PetscReal          hx = user->hx, hy = user->hy, area, three = 3.0, p5 = 0.5;
    261:   PetscReal          zero = 0.0, vb, vl, vr, vt, dvdx, dvdy, flin = 0.0, fquad = 0.0;
    262:   PetscReal          v, cdiv3 = user->param / three;
    263:   const PetscScalar *x;
    264:   PetscInt           nx = user->mx, ny = user->my, i, j, k;

    266:   PetscFunctionBeginUser;
    267:   /* Get pointer to vector data */
    268:   PetscCall(VecGetArrayRead(X, &x));

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

    287:   /* Compute function contributions over the upper triangular elements */
    288:   for (j = 0; j <= ny; j++) {
    289:     for (i = 0; i <= nx; i++) {
    290:       k  = nx * j + i;
    291:       vb = zero;
    292:       vl = zero;
    293:       v  = zero;
    294:       if (i < nx && j > 0) vb = x[k - nx];
    295:       if (i > 0 && j < ny) vl = x[k - 1];
    296:       if (i < nx && j < ny) v = x[k];
    297:       dvdx  = (v - vl) / hx;
    298:       dvdy  = (v - vb) / hy;
    299:       fquad = fquad + dvdx * dvdx + dvdy * dvdy;
    300:       flin  = flin - cdiv3 * (vb + vl + v);
    301:     }
    302:   }

    304:   /* Restore vector */
    305:   PetscCall(VecRestoreArrayRead(X, &x));

    307:   /* Scale the function */
    308:   area = p5 * hx * hy;
    309:   *f   = area * (p5 * fquad + flin);

    311:   PetscCall(PetscLogFlops(24.0 * nx * ny));
    312:   PetscFunctionReturn(PETSC_SUCCESS);
    313: }

    315: /* ------------------------------------------------------------------- */
    316: /*
    317:     FormGradient - Evaluates the gradient, G(X).

    319:     Input Parameters:
    320: .   tao  - the Tao context
    321: .   X    - input vector
    322: .   ptr  - optional user-defined context

    324:     Output Parameters:
    325: .   G - vector containing the newly evaluated gradient
    326: */
    327: PetscErrorCode FormGradient(Tao tao, Vec X, Vec G, void *ptr)
    328: {
    329:   AppCtx            *user = (AppCtx *)ptr;
    330:   PetscReal          zero = 0.0, p5 = 0.5, three = 3.0, area, val;
    331:   PetscInt           nx = user->mx, ny = user->my, ind, i, j, k;
    332:   PetscReal          hx = user->hx, hy = user->hy;
    333:   PetscReal          vb, vl, vr, vt, dvdx, dvdy;
    334:   PetscReal          v, cdiv3 = user->param / three;
    335:   const PetscScalar *x;

    337:   PetscFunctionBeginUser;
    338:   /* Initialize gradient to zero */
    339:   PetscCall(VecSet(G, zero));

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

    344:   /* Compute gradient contributions over the lower triangular elements */
    345:   for (j = -1; j < ny; j++) {
    346:     for (i = -1; i < nx; i++) {
    347:       k  = nx * j + i;
    348:       v  = zero;
    349:       vr = zero;
    350:       vt = zero;
    351:       if (i >= 0 && j >= 0) v = x[k];
    352:       if (i < nx - 1 && j > -1) vr = x[k + 1];
    353:       if (i > -1 && j < ny - 1) vt = x[k + nx];
    354:       dvdx = (vr - v) / hx;
    355:       dvdy = (vt - v) / hy;
    356:       if (i != -1 && j != -1) {
    357:         ind = k;
    358:         val = -dvdx / hx - dvdy / hy - cdiv3;
    359:         PetscCall(VecSetValues(G, 1, &ind, &val, ADD_VALUES));
    360:       }
    361:       if (i != nx - 1 && j != -1) {
    362:         ind = k + 1;
    363:         val = dvdx / hx - cdiv3;
    364:         PetscCall(VecSetValues(G, 1, &ind, &val, ADD_VALUES));
    365:       }
    366:       if (i != -1 && j != ny - 1) {
    367:         ind = k + nx;
    368:         val = dvdy / hy - cdiv3;
    369:         PetscCall(VecSetValues(G, 1, &ind, &val, ADD_VALUES));
    370:       }
    371:     }
    372:   }

    374:   /* Compute gradient contributions over the upper triangular elements */
    375:   for (j = 0; j <= ny; j++) {
    376:     for (i = 0; i <= nx; i++) {
    377:       k  = nx * j + i;
    378:       vb = zero;
    379:       vl = zero;
    380:       v  = zero;
    381:       if (i < nx && j > 0) vb = x[k - nx];
    382:       if (i > 0 && j < ny) vl = x[k - 1];
    383:       if (i < nx && j < ny) v = x[k];
    384:       dvdx = (v - vl) / hx;
    385:       dvdy = (v - vb) / hy;
    386:       if (i != nx && j != 0) {
    387:         ind = k - nx;
    388:         val = -dvdy / hy - cdiv3;
    389:         PetscCall(VecSetValues(G, 1, &ind, &val, ADD_VALUES));
    390:       }
    391:       if (i != 0 && j != ny) {
    392:         ind = k - 1;
    393:         val = -dvdx / hx - cdiv3;
    394:         PetscCall(VecSetValues(G, 1, &ind, &val, ADD_VALUES));
    395:       }
    396:       if (i != nx && j != ny) {
    397:         ind = k;
    398:         val = dvdx / hx + dvdy / hy - cdiv3;
    399:         PetscCall(VecSetValues(G, 1, &ind, &val, ADD_VALUES));
    400:       }
    401:     }
    402:   }
    403:   PetscCall(VecRestoreArrayRead(X, &x));

    405:   /* Assemble gradient vector */
    406:   PetscCall(VecAssemblyBegin(G));
    407:   PetscCall(VecAssemblyEnd(G));

    409:   /* Scale the gradient */
    410:   area = p5 * hx * hy;
    411:   PetscCall(VecScale(G, area));
    412:   PetscCall(PetscLogFlops(24.0 * nx * ny));
    413:   PetscFunctionReturn(PETSC_SUCCESS);
    414: }

    416: /* ------------------------------------------------------------------- */
    417: /*
    418:    FormHessian - Forms the Hessian matrix.

    420:    Input Parameters:
    421: .  tao - the Tao context
    422: .  X    - the input vector
    423: .  ptr  - optional user-defined context, as set by TaoSetHessian()

    425:    Output Parameters:
    426: .  H     - Hessian matrix
    427: .  PrecH - optionally different preconditioning Hessian

    429:    Notes:
    430:    This routine is intended simply as an example of the interface
    431:    to a Hessian evaluation routine.  Since this example compute the
    432:    Hessian a column at a time, it is not particularly efficient and
    433:    is not recommended.
    434: */
    435: PetscErrorCode FormHessian(Tao tao, Vec X, Mat H, Mat Hpre, void *ptr)
    436: {
    437:   AppCtx    *user = (AppCtx *)ptr;
    438:   PetscInt   i, j, ndim = user->ndim;
    439:   PetscReal *y, zero = 0.0, one = 1.0;
    440:   PetscBool  assembled;

    442:   PetscFunctionBeginUser;
    443:   user->xvec = X;

    445:   /* Initialize Hessian entries and work vector to zero */
    446:   PetscCall(MatAssembled(H, &assembled));
    447:   if (assembled) PetscCall(MatZeroEntries(H));

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

    451:   /* Loop over matrix columns to compute entries of the Hessian */
    452:   for (j = 0; j < ndim; j++) {
    453:     PetscCall(VecSetValues(user->s, 1, &j, &one, INSERT_VALUES));
    454:     PetscCall(VecAssemblyBegin(user->s));
    455:     PetscCall(VecAssemblyEnd(user->s));

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

    459:     PetscCall(VecSetValues(user->s, 1, &j, &zero, INSERT_VALUES));
    460:     PetscCall(VecAssemblyBegin(user->s));
    461:     PetscCall(VecAssemblyEnd(user->s));

    463:     PetscCall(VecGetArray(user->y, &y));
    464:     for (i = 0; i < ndim; i++) {
    465:       if (y[i] != zero) PetscCall(MatSetValues(H, 1, &i, 1, &j, &y[i], ADD_VALUES));
    466:     }
    467:     PetscCall(VecRestoreArray(user->y, &y));
    468:   }
    469:   PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
    470:   PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
    471:   PetscFunctionReturn(PETSC_SUCCESS);
    472: }

    474: /* ------------------------------------------------------------------- */
    475: /*
    476:    MatrixFreeHessian - Sets a pointer for use in computing Hessian-vector
    477:    products.

    479:    Input Parameters:
    480: .  tao - the Tao context
    481: .  X    - the input vector
    482: .  ptr  - optional user-defined context, as set by TaoSetHessian()

    484:    Output Parameters:
    485: .  H     - Hessian matrix
    486: .  PrecH - optionally different preconditioning Hessian
    487: */
    488: PetscErrorCode MatrixFreeHessian(Tao tao, Vec X, Mat H, Mat PrecH, void *ptr)
    489: {
    490:   AppCtx *user = (AppCtx *)ptr;

    492:   /* Sets location of vector for use in computing matrix-vector products  of the form H(X)*y  */
    493:   PetscFunctionBeginUser;
    494:   user->xvec = X;
    495:   PetscFunctionReturn(PETSC_SUCCESS);
    496: }

    498: /* ------------------------------------------------------------------- */
    499: /*
    500:    HessianProductMat - Computes the matrix-vector product
    501:    y = mat*svec.

    503:    Input Parameters:
    504: .  mat  - input matrix
    505: .  svec - input vector

    507:    Output Parameters:
    508: .  y    - solution vector
    509: */
    510: PetscErrorCode HessianProductMat(Mat mat, Vec svec, Vec y)
    511: {
    512:   void *ptr;

    514:   PetscFunctionBeginUser;
    515:   PetscCall(MatShellGetContext(mat, &ptr));
    516:   PetscCall(HessianProduct(ptr, svec, y));
    517:   PetscFunctionReturn(PETSC_SUCCESS);
    518: }

    520: /* ------------------------------------------------------------------- */
    521: /*
    522:    Hessian Product - Computes the matrix-vector product:
    523:    y = f''(x)*svec.

    525:    Input Parameters:
    526: .  ptr  - pointer to the user-defined context
    527: .  svec - input vector

    529:    Output Parameters:
    530: .  y    - product vector
    531: */
    532: PetscErrorCode HessianProduct(void *ptr, Vec svec, Vec y)
    533: {
    534:   AppCtx            *user = (AppCtx *)ptr;
    535:   PetscReal          p5 = 0.5, zero = 0.0, one = 1.0, hx, hy, val, area;
    536:   const PetscScalar *x, *s;
    537:   PetscReal          v, vb, vl, vr, vt, hxhx, hyhy;
    538:   PetscInt           nx, ny, i, j, k, ind;

    540:   PetscFunctionBeginUser;
    541:   nx   = user->mx;
    542:   ny   = user->my;
    543:   hx   = user->hx;
    544:   hy   = user->hy;
    545:   hxhx = one / (hx * hx);
    546:   hyhy = one / (hy * hy);

    548:   /* Get pointers to vector data */
    549:   PetscCall(VecGetArrayRead(user->xvec, &x));
    550:   PetscCall(VecGetArrayRead(svec, &s));

    552:   /* Initialize product vector to zero */
    553:   PetscCall(VecSet(y, zero));

    555:   /* Compute f''(x)*s over the lower triangular elements */
    556:   for (j = -1; j < ny; j++) {
    557:     for (i = -1; i < nx; i++) {
    558:       k  = nx * j + i;
    559:       v  = zero;
    560:       vr = zero;
    561:       vt = zero;
    562:       if (i != -1 && j != -1) v = s[k];
    563:       if (i != nx - 1 && j != -1) {
    564:         vr  = s[k + 1];
    565:         ind = k + 1;
    566:         val = hxhx * (vr - v);
    567:         PetscCall(VecSetValues(y, 1, &ind, &val, ADD_VALUES));
    568:       }
    569:       if (i != -1 && j != ny - 1) {
    570:         vt  = s[k + nx];
    571:         ind = k + nx;
    572:         val = hyhy * (vt - v);
    573:         PetscCall(VecSetValues(y, 1, &ind, &val, ADD_VALUES));
    574:       }
    575:       if (i != -1 && j != -1) {
    576:         ind = k;
    577:         val = hxhx * (v - vr) + hyhy * (v - vt);
    578:         PetscCall(VecSetValues(y, 1, &ind, &val, ADD_VALUES));
    579:       }
    580:     }
    581:   }

    583:   /* Compute f''(x)*s over the upper triangular elements */
    584:   for (j = 0; j <= ny; j++) {
    585:     for (i = 0; i <= nx; i++) {
    586:       k  = nx * j + i;
    587:       v  = zero;
    588:       vl = zero;
    589:       vb = zero;
    590:       if (i != nx && j != ny) v = s[k];
    591:       if (i != nx && j != 0) {
    592:         vb  = s[k - nx];
    593:         ind = k - nx;
    594:         val = hyhy * (vb - v);
    595:         PetscCall(VecSetValues(y, 1, &ind, &val, ADD_VALUES));
    596:       }
    597:       if (i != 0 && j != ny) {
    598:         vl  = s[k - 1];
    599:         ind = k - 1;
    600:         val = hxhx * (vl - v);
    601:         PetscCall(VecSetValues(y, 1, &ind, &val, ADD_VALUES));
    602:       }
    603:       if (i != nx && j != ny) {
    604:         ind = k;
    605:         val = hxhx * (v - vl) + hyhy * (v - vb);
    606:         PetscCall(VecSetValues(y, 1, &ind, &val, ADD_VALUES));
    607:       }
    608:     }
    609:   }
    610:   /* Restore vector data */
    611:   PetscCall(VecRestoreArrayRead(svec, &s));
    612:   PetscCall(VecRestoreArrayRead(user->xvec, &x));

    614:   /* Assemble vector */
    615:   PetscCall(VecAssemblyBegin(y));
    616:   PetscCall(VecAssemblyEnd(y));

    618:   /* Scale resulting vector by area */
    619:   area = p5 * hx * hy;
    620:   PetscCall(VecScale(y, area));
    621:   PetscCall(PetscLogFlops(18.0 * nx * ny));
    622:   PetscFunctionReturn(PETSC_SUCCESS);
    623: }

    625: /*TEST

    627:    build:
    628:       requires: !complex

    630:    test:
    631:       args: -tao_monitor -tao_type bntr -tao_view -tao_bnk_ksp_monitor_short

    633: TEST*/