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

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

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

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

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

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

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

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

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

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

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

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

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

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

489:    Input Parameters:
490: .  mat  - input matrix
491: .  svec - input vector

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

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

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

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

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

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

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

538:   /* Initialize product vector to zero */
539:   PetscCall(VecSet(y, zero));

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

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

600:   /* Assemble vector */
601:   PetscCall(VecAssemblyBegin(y));
602:   PetscCall(VecAssemblyEnd(y));

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

611: /*TEST

613:    build:
614:       requires: !complex

616:    test:
617:       suffix: 1
618:       args: -tao_monitor_short -tao_type ntl -tao_gatol 1.e-4

620:    test:
621:       suffix: 2
622:       args: -tao_monitor_short -tao_type ntr -tao_gatol 1.e-4

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

628:    test:
629:      suffix: 4
630:      args: -tao_monitor_short -tao_gatol 1e-3 -tao_type bqnls

632:    test:
633:      suffix: 5
634:      args: -tao_monitor_short -tao_gatol 1e-3 -tao_type blmvm

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

640:    test:
641:      suffix: snes
642:      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

644:    test:
645:      suffix: snes_2
646:      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

648:    test:
649:      suffix: snes_3
650:      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

652: TEST*/