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: petscsys.h - system routines
24: petscis.h - index sets
25: petscvec.h - vectors
26: petscmat.h - matrices
27: petscpc.h - preconditioners
28: petscksp.h - Krylov subspace methods
29: */
31: #include <petsctao.h>
33: static char help[] = "Demonstrates use of the TAO package to solve \n\
34: unconstrained minimization problems on a single processor. This example \n\
35: is based on the Elastic-Plastic Torsion (dept) problem from the MINPACK-2 \n\
36: test suite.\n\
37: The command line options are:\n\
38: -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
39: -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
40: -par <param>, where <param> = angle of twist per unit length\n\n";
42: /*
43: User-defined application context - contains data needed by the
44: application-provided call-back routines, FormFunction(),
45: FormGradient(), and FormHessian().
46: */
48: typedef struct {
49: PetscReal param; /* nonlinearity parameter */
50: PetscInt mx, my; /* discretization in x- and y-directions */
51: PetscInt ndim; /* problem dimension */
52: Vec s, y, xvec; /* work space for computing Hessian */
53: PetscReal hx, hy; /* mesh spacing in x- and y-directions */
54: } AppCtx;
56: /* -------- User-defined Routines --------- */
58: PetscErrorCode FormInitialGuess(AppCtx *, Vec);
59: PetscErrorCode FormFunction(Tao, Vec, PetscReal *, void *);
60: PetscErrorCode FormGradient(Tao, Vec, Vec, void *);
61: PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void *);
62: PetscErrorCode HessianProductMat(Mat, Vec, Vec);
63: PetscErrorCode HessianProduct(void *, Vec, Vec);
64: PetscErrorCode MatrixFreeHessian(Tao, Vec, Mat, Mat, void *);
65: PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);
67: int main(int argc, char **argv)
68: {
69: PetscInt mx = 10; /* discretization in x-direction */
70: PetscInt my = 10; /* discretization in y-direction */
71: Vec x; /* solution, gradient vectors */
72: PetscBool flg; /* A return value when checking for use options */
73: Tao tao; /* Tao solver context */
74: Mat H; /* Hessian matrix */
75: AppCtx user; /* application context */
76: PetscMPIInt size; /* number of processes */
77: PetscReal one = 1.0;
79: PetscBool test_lmvm = PETSC_FALSE;
80: KSP ksp;
81: PC pc;
82: Mat M;
83: Vec in, out, out2;
84: PetscReal mult_solve_dist;
85: Vec ub, lb;
86: PetscInt i = 3;
88: /* Initialize TAO,PETSc */
89: PetscFunctionBeginUser;
90: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
91: PetscCallMPI(MPI_Comm_size(MPI_COMM_WORLD, &size));
92: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Incorrect number of processors");
94: /* Specify default parameters for the problem, check for command-line overrides */
95: user.param = 5.0;
96: PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &my, &flg));
97: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &mx, &flg));
98: PetscCall(PetscOptionsGetReal(NULL, NULL, "-par", &user.param, &flg));
99: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_lmvm", &test_lmvm, &flg));
101: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n---- Elastic-Plastic Torsion Problem -----\n"));
102: PetscCall(PetscPrintf(PETSC_COMM_SELF, "mx: %" PetscInt_FMT " my: %" PetscInt_FMT " \n\n", mx, my));
103: user.ndim = mx * my;
104: user.mx = mx;
105: user.my = my;
106: user.hx = one / (mx + 1);
107: user.hy = one / (my + 1);
109: /* Allocate vectors */
110: PetscCall(VecCreateSeq(PETSC_COMM_SELF, user.ndim, &user.y));
111: PetscCall(VecDuplicate(user.y, &user.s));
112: PetscCall(VecDuplicate(user.y, &x));
113: PetscCall(VecDuplicate(user.y, &lb));
114: PetscCall(VecDuplicate(user.y, &ub));
116: PetscCall(VecSet(ub, 0.1));
117: PetscCall(VecSet(lb, -0.1));
118: PetscCall(VecSetValue(ub, i, 0., INSERT_VALUES));
119: PetscCall(VecSetValue(lb, i, 0., INSERT_VALUES));
120: PetscCall(VecAssemblyBegin(ub));
121: PetscCall(VecAssemblyBegin(lb));
122: PetscCall(VecAssemblyEnd(ub));
123: PetscCall(VecAssemblyEnd(lb));
125: /* Create TAO solver and set desired solution method */
126: PetscCall(TaoCreate(PETSC_COMM_SELF, &tao));
127: PetscCall(TaoSetType(tao, TAOLMVM));
129: /* Set solution vector with an initial guess */
130: PetscCall(FormInitialGuess(&user, x));
131: PetscCall(TaoSetSolution(tao, x));
132: PetscCall(TaoSetVariableBounds(tao, lb, ub));
134: /* Set routine for function and gradient evaluation */
135: PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user));
137: /* From command line options, determine if using matrix-free hessian */
138: PetscCall(PetscOptionsHasName(NULL, NULL, "-my_tao_mf", &flg));
139: if (flg) {
140: PetscCall(MatCreateShell(PETSC_COMM_SELF, user.ndim, user.ndim, user.ndim, user.ndim, (void *)&user, &H));
141: PetscCall(MatShellSetOperation(H, MATOP_MULT, (PetscErrorCodeFn *)HessianProductMat));
142: PetscCall(MatSetOption(H, MAT_SYMMETRIC, PETSC_TRUE));
144: PetscCall(TaoSetHessian(tao, H, H, MatrixFreeHessian, (void *)&user));
145: } else {
146: PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, user.ndim, user.ndim, 5, NULL, &H));
147: PetscCall(MatSetOption(H, MAT_SYMMETRIC, PETSC_TRUE));
148: PetscCall(TaoSetHessian(tao, H, H, FormHessian, (void *)&user));
149: }
151: /* Test the LMVM matrix */
152: if (test_lmvm) {
153: PetscCall(PetscOptionsSetValue(NULL, "-tao_type", "bntr"));
154: PetscCall(PetscOptionsSetValue(NULL, "-tao_bnk_pc_type", "lmvm"));
155: }
157: /* Check for any TAO command line options */
158: PetscCall(TaoSetFromOptions(tao));
160: /* SOLVE THE APPLICATION */
161: PetscCall(TaoSolve(tao));
163: /* Test the LMVM matrix */
164: if (test_lmvm) {
165: PetscCall(TaoGetKSP(tao, &ksp));
166: PetscCall(KSPGetPC(ksp, &pc));
167: PetscCall(PCLMVMGetMatLMVM(pc, &M));
168: PetscCall(VecDuplicate(x, &in));
169: PetscCall(VecDuplicate(x, &out));
170: PetscCall(VecDuplicate(x, &out2));
171: PetscCall(VecSet(in, 5.0));
172: PetscCall(MatMult(M, in, out));
173: PetscCall(MatSolve(M, out, out2));
174: PetscCall(VecAXPY(out2, -1.0, in));
175: PetscCall(VecNorm(out2, NORM_2, &mult_solve_dist));
176: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)tao), "error between MatMult and MatSolve: %e\n", (double)mult_solve_dist));
177: PetscCall(VecDestroy(&in));
178: PetscCall(VecDestroy(&out));
179: PetscCall(VecDestroy(&out2));
180: }
182: PetscCall(TaoDestroy(&tao));
183: PetscCall(VecDestroy(&user.s));
184: PetscCall(VecDestroy(&user.y));
185: PetscCall(VecDestroy(&x));
186: PetscCall(MatDestroy(&H));
187: PetscCall(VecDestroy(&lb));
188: PetscCall(VecDestroy(&ub));
190: PetscCall(PetscFinalize());
191: return 0;
192: }
194: /* ------------------------------------------------------------------- */
195: /*
196: FormInitialGuess - Computes an initial approximation to the solution.
198: Input Parameters:
199: . user - user-defined application context
200: . X - vector
202: Output Parameters:
203: . X - vector
204: */
205: PetscErrorCode FormInitialGuess(AppCtx *user, Vec X)
206: {
207: PetscReal hx = user->hx, hy = user->hy, temp;
208: PetscReal val;
209: PetscInt i, j, k, nx = user->mx, ny = user->my;
211: /* Compute initial guess */
212: PetscFunctionBeginUser;
213: for (j = 0; j < ny; j++) {
214: temp = PetscMin(j + 1, ny - j) * hy;
215: for (i = 0; i < nx; i++) {
216: k = nx * j + i;
217: val = PetscMin((PetscMin(i + 1, nx - i)) * hx, temp);
218: PetscCall(VecSetValues(X, 1, &k, &val, ADD_VALUES));
219: }
220: }
221: PetscCall(VecAssemblyBegin(X));
222: PetscCall(VecAssemblyEnd(X));
223: PetscFunctionReturn(PETSC_SUCCESS);
224: }
226: /* ------------------------------------------------------------------- */
227: /*
228: FormFunctionGradient - Evaluates the function and corresponding gradient.
230: Input Parameters:
231: tao - the Tao context
232: X - the input vector
233: ptr - optional user-defined context, as set by TaoSetFunction()
235: Output Parameters:
236: f - the newly evaluated function
237: G - the newly evaluated gradient
238: */
239: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ptr)
240: {
241: PetscFunctionBeginUser;
242: PetscCall(FormFunction(tao, X, f, ptr));
243: PetscCall(FormGradient(tao, X, G, ptr));
244: PetscFunctionReturn(PETSC_SUCCESS);
245: }
247: /* ------------------------------------------------------------------- */
248: /*
249: FormFunction - Evaluates the function, f(X).
251: Input Parameters:
252: . tao - the Tao context
253: . X - the input vector
254: . ptr - optional user-defined context, as set by TaoSetFunction()
256: Output Parameters:
257: . f - the newly evaluated function
258: */
259: PetscErrorCode FormFunction(Tao tao, Vec X, PetscReal *f, void *ptr)
260: {
261: AppCtx *user = (AppCtx *)ptr;
262: PetscReal hx = user->hx, hy = user->hy, area, three = 3.0, p5 = 0.5;
263: PetscReal zero = 0.0, vb, vl, vr, vt, dvdx, dvdy, flin = 0.0, fquad = 0.0;
264: PetscReal v, cdiv3 = user->param / three;
265: const PetscScalar *x;
266: PetscInt nx = user->mx, ny = user->my, i, j, k;
268: PetscFunctionBeginUser;
269: /* Get pointer to vector data */
270: PetscCall(VecGetArrayRead(X, &x));
272: /* Compute function contributions over the lower triangular elements */
273: for (j = -1; j < ny; j++) {
274: for (i = -1; i < nx; i++) {
275: k = nx * j + i;
276: v = zero;
277: vr = zero;
278: vt = zero;
279: if (i >= 0 && j >= 0) v = x[k];
280: if (i < nx - 1 && j > -1) vr = x[k + 1];
281: if (i > -1 && j < ny - 1) vt = x[k + nx];
282: dvdx = (vr - v) / hx;
283: dvdy = (vt - v) / hy;
284: fquad += dvdx * dvdx + dvdy * dvdy;
285: flin -= cdiv3 * (v + vr + vt);
286: }
287: }
289: /* Compute function contributions over the upper triangular elements */
290: for (j = 0; j <= ny; j++) {
291: for (i = 0; i <= nx; i++) {
292: k = nx * j + i;
293: vb = zero;
294: vl = zero;
295: v = zero;
296: if (i < nx && j > 0) vb = x[k - nx];
297: if (i > 0 && j < ny) vl = x[k - 1];
298: if (i < nx && j < ny) v = x[k];
299: dvdx = (v - vl) / hx;
300: dvdy = (v - vb) / hy;
301: fquad = fquad + dvdx * dvdx + dvdy * dvdy;
302: flin = flin - cdiv3 * (vb + vl + v);
303: }
304: }
306: /* Restore vector */
307: PetscCall(VecRestoreArrayRead(X, &x));
309: /* Scale the function */
310: area = p5 * hx * hy;
311: *f = area * (p5 * fquad + flin);
313: PetscCall(PetscLogFlops(24.0 * nx * ny));
314: PetscFunctionReturn(PETSC_SUCCESS);
315: }
317: /* ------------------------------------------------------------------- */
318: /*
319: FormGradient - Evaluates the gradient, G(X).
321: Input Parameters:
322: . tao - the Tao context
323: . X - input vector
324: . ptr - optional user-defined context
326: Output Parameters:
327: . G - vector containing the newly evaluated gradient
328: */
329: PetscErrorCode FormGradient(Tao tao, Vec X, Vec G, void *ptr)
330: {
331: AppCtx *user = (AppCtx *)ptr;
332: PetscReal zero = 0.0, p5 = 0.5, three = 3.0, area, val;
333: PetscInt nx = user->mx, ny = user->my, ind, i, j, k;
334: PetscReal hx = user->hx, hy = user->hy;
335: PetscReal vb, vl, vr, vt, dvdx, dvdy;
336: PetscReal v, cdiv3 = user->param / three;
337: const PetscScalar *x;
339: PetscFunctionBeginUser;
340: /* Initialize gradient to zero */
341: PetscCall(VecSet(G, zero));
343: /* Get pointer to vector data */
344: PetscCall(VecGetArrayRead(X, &x));
346: /* Compute gradient contributions over the lower triangular elements */
347: for (j = -1; j < ny; j++) {
348: for (i = -1; i < nx; i++) {
349: k = nx * j + i;
350: v = zero;
351: vr = zero;
352: vt = zero;
353: if (i >= 0 && j >= 0) v = x[k];
354: if (i < nx - 1 && j > -1) vr = x[k + 1];
355: if (i > -1 && j < ny - 1) vt = x[k + nx];
356: dvdx = (vr - v) / hx;
357: dvdy = (vt - v) / hy;
358: if (i != -1 && j != -1) {
359: ind = k;
360: val = -dvdx / hx - dvdy / hy - cdiv3;
361: PetscCall(VecSetValues(G, 1, &ind, &val, ADD_VALUES));
362: }
363: if (i != nx - 1 && j != -1) {
364: ind = k + 1;
365: val = dvdx / hx - cdiv3;
366: PetscCall(VecSetValues(G, 1, &ind, &val, ADD_VALUES));
367: }
368: if (i != -1 && j != ny - 1) {
369: ind = k + nx;
370: val = dvdy / hy - cdiv3;
371: PetscCall(VecSetValues(G, 1, &ind, &val, ADD_VALUES));
372: }
373: }
374: }
376: /* Compute gradient contributions over the upper triangular elements */
377: for (j = 0; j <= ny; j++) {
378: for (i = 0; i <= nx; i++) {
379: k = nx * j + i;
380: vb = zero;
381: vl = zero;
382: v = zero;
383: if (i < nx && j > 0) vb = x[k - nx];
384: if (i > 0 && j < ny) vl = x[k - 1];
385: if (i < nx && j < ny) v = x[k];
386: dvdx = (v - vl) / hx;
387: dvdy = (v - vb) / hy;
388: if (i != nx && j != 0) {
389: ind = k - nx;
390: val = -dvdy / hy - cdiv3;
391: PetscCall(VecSetValues(G, 1, &ind, &val, ADD_VALUES));
392: }
393: if (i != 0 && j != ny) {
394: ind = k - 1;
395: val = -dvdx / hx - cdiv3;
396: PetscCall(VecSetValues(G, 1, &ind, &val, ADD_VALUES));
397: }
398: if (i != nx && j != ny) {
399: ind = k;
400: val = dvdx / hx + dvdy / hy - cdiv3;
401: PetscCall(VecSetValues(G, 1, &ind, &val, ADD_VALUES));
402: }
403: }
404: }
405: PetscCall(VecRestoreArrayRead(X, &x));
407: /* Assemble gradient vector */
408: PetscCall(VecAssemblyBegin(G));
409: PetscCall(VecAssemblyEnd(G));
411: /* Scale the gradient */
412: area = p5 * hx * hy;
413: PetscCall(VecScale(G, area));
414: PetscCall(PetscLogFlops(24.0 * nx * ny));
415: PetscFunctionReturn(PETSC_SUCCESS);
416: }
418: /* ------------------------------------------------------------------- */
419: /*
420: FormHessian - Forms the Hessian matrix.
422: Input Parameters:
423: . tao - the Tao context
424: . X - the input vector
425: . ptr - optional user-defined context, as set by TaoSetHessian()
427: Output Parameters:
428: . H - Hessian matrix
429: . PrecH - optionally different preconditioning Hessian
431: Notes:
432: This routine is intended simply as an example of the interface
433: to a Hessian evaluation routine. Since this example compute the
434: Hessian a column at a time, it is not particularly efficient and
435: is not recommended.
436: */
437: PetscErrorCode FormHessian(Tao tao, Vec X, Mat H, Mat Hpre, void *ptr)
438: {
439: AppCtx *user = (AppCtx *)ptr;
440: PetscInt i, j, ndim = user->ndim;
441: PetscReal *y, zero = 0.0, one = 1.0;
442: PetscBool assembled;
444: PetscFunctionBeginUser;
445: user->xvec = X;
447: /* Initialize Hessian entries and work vector to zero */
448: PetscCall(MatAssembled(H, &assembled));
449: if (assembled) PetscCall(MatZeroEntries(H));
451: PetscCall(VecSet(user->s, zero));
453: /* Loop over matrix columns to compute entries of the Hessian */
454: for (j = 0; j < ndim; j++) {
455: PetscCall(VecSetValues(user->s, 1, &j, &one, INSERT_VALUES));
456: PetscCall(VecAssemblyBegin(user->s));
457: PetscCall(VecAssemblyEnd(user->s));
459: PetscCall(HessianProduct(ptr, user->s, user->y));
461: PetscCall(VecSetValues(user->s, 1, &j, &zero, INSERT_VALUES));
462: PetscCall(VecAssemblyBegin(user->s));
463: PetscCall(VecAssemblyEnd(user->s));
465: PetscCall(VecGetArray(user->y, &y));
466: for (i = 0; i < ndim; i++) {
467: if (y[i] != zero) PetscCall(MatSetValues(H, 1, &i, 1, &j, &y[i], ADD_VALUES));
468: }
469: PetscCall(VecRestoreArray(user->y, &y));
470: }
471: PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
472: PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
473: PetscFunctionReturn(PETSC_SUCCESS);
474: }
476: /* ------------------------------------------------------------------- */
477: /*
478: MatrixFreeHessian - Sets a pointer for use in computing Hessian-vector
479: products.
481: Input Parameters:
482: . tao - the Tao context
483: . X - the input vector
484: . ptr - optional user-defined context, as set by TaoSetHessian()
486: Output Parameters:
487: . H - Hessian matrix
488: . PrecH - optionally different preconditioning Hessian
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
635: TEST*/