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, (char *)0, 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
428: .  flag  - flag indicating matrix structure

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

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

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

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

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

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

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

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

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

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

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

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

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

505:    Input Parameters:
506: .  mat  - input matrix
507: .  svec - input vector

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

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

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

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

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

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

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

554:   /* Initialize product vector to zero */
555:   PetscCall(VecSet(y, zero));

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

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

616:   /* Assemble vector */
617:   PetscCall(VecAssemblyBegin(y));
618:   PetscCall(VecAssemblyEnd(y));

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

627: /*TEST

629:    build:
630:       requires: !complex

632:    test:
633:       args: -tao_monitor -tao_type bntr  -tao_view -tao_bnk_ksp_monitor_short

635: TEST*/