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*/