1 /* Program usage: mpiexec -n 1 eptorsion1 [-help] [all TAO options] */
2
3 /* ----------------------------------------------------------------------
4
5 Elastic-plastic torsion problem.
6
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:
10
11 min{ .5 * integral(||gradient(v(x))||^2 dx) - C * integral(v(x) dx)}
12
13 where C is the torsion angle per unit length.
14
15 The multiprocessor version of this code is eptorsion2.c.
16
17 ---------------------------------------------------------------------- */
18
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 */
28
29 #include <petsctao.h>
30
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";
39
40 /*
41 User-defined application context - contains data needed by the
42 application-provided call-back routines, FormFunction(),
43 FormGradient(), and FormHessian().
44 */
45
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;
53
54 /* -------- User-defined Routines --------- */
55
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 *);
64
main(int argc,char ** argv)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;
76
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;
85
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");
91
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));
98
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);
106
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));
113
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));
122
123 /* Create TAO solver and set desired solution method */
124 PetscCall(TaoCreate(PETSC_COMM_SELF, &tao));
125 PetscCall(TaoSetType(tao, TAOLMVM));
126
127 /* Set solution vector with an initial guess */
128 PetscCall(FormInitialGuess(&user, x));
129 PetscCall(TaoSetSolution(tao, x));
130 PetscCall(TaoSetVariableBounds(tao, lb, ub));
131
132 /* Set routine for function and gradient evaluation */
133 PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user));
134
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, (PetscErrorCodeFn *)HessianProductMat));
140 PetscCall(MatSetOption(H, MAT_SYMMETRIC, PETSC_TRUE));
141
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 }
148
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 }
154
155 /* Check for any TAO command line options */
156 PetscCall(TaoSetFromOptions(tao));
157
158 /* SOLVE THE APPLICATION */
159 PetscCall(TaoSolve(tao));
160
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 }
179
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));
187
188 PetscCall(PetscFinalize());
189 return 0;
190 }
191
192 /* ------------------------------------------------------------------- */
193 /*
194 FormInitialGuess - Computes an initial approximation to the solution.
195
196 Input Parameters:
197 . user - user-defined application context
198 . X - vector
199
200 Output Parameters:
201 . X - vector
202 */
FormInitialGuess(AppCtx * user,Vec X)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;
208
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 }
223
224 /* ------------------------------------------------------------------- */
225 /*
226 FormFunctionGradient - Evaluates the function and corresponding gradient.
227
228 Input Parameters:
229 tao - the Tao context
230 X - the input vector
231 ptr - optional user-defined context, as set by TaoSetFunction()
232
233 Output Parameters:
234 f - the newly evaluated function
235 G - the newly evaluated gradient
236 */
FormFunctionGradient(Tao tao,Vec X,PetscReal * f,Vec G,void * ptr)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 }
244
245 /* ------------------------------------------------------------------- */
246 /*
247 FormFunction - Evaluates the function, f(X).
248
249 Input Parameters:
250 . tao - the Tao context
251 . X - the input vector
252 . ptr - optional user-defined context, as set by TaoSetFunction()
253
254 Output Parameters:
255 . f - the newly evaluated function
256 */
FormFunction(Tao tao,Vec X,PetscReal * f,void * ptr)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;
265
266 PetscFunctionBeginUser;
267 /* Get pointer to vector data */
268 PetscCall(VecGetArrayRead(X, &x));
269
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 }
286
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 }
303
304 /* Restore vector */
305 PetscCall(VecRestoreArrayRead(X, &x));
306
307 /* Scale the function */
308 area = p5 * hx * hy;
309 *f = area * (p5 * fquad + flin);
310
311 PetscCall(PetscLogFlops(24.0 * nx * ny));
312 PetscFunctionReturn(PETSC_SUCCESS);
313 }
314
315 /* ------------------------------------------------------------------- */
316 /*
317 FormGradient - Evaluates the gradient, G(X).
318
319 Input Parameters:
320 . tao - the Tao context
321 . X - input vector
322 . ptr - optional user-defined context
323
324 Output Parameters:
325 . G - vector containing the newly evaluated gradient
326 */
FormGradient(Tao tao,Vec X,Vec G,void * ptr)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;
336
337 PetscFunctionBeginUser;
338 /* Initialize gradient to zero */
339 PetscCall(VecSet(G, zero));
340
341 /* Get pointer to vector data */
342 PetscCall(VecGetArrayRead(X, &x));
343
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 }
373
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));
404
405 /* Assemble gradient vector */
406 PetscCall(VecAssemblyBegin(G));
407 PetscCall(VecAssemblyEnd(G));
408
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 }
415
416 /* ------------------------------------------------------------------- */
417 /*
418 FormHessian - Forms the Hessian matrix.
419
420 Input Parameters:
421 . tao - the Tao context
422 . X - the input vector
423 . ptr - optional user-defined context, as set by TaoSetHessian()
424
425 Output Parameters:
426 . H - Hessian matrix
427 . PrecH - optionally different preconditioning Hessian
428
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 */
FormHessian(Tao tao,Vec X,Mat H,Mat Hpre,void * ptr)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;
441
442 PetscFunctionBeginUser;
443 user->xvec = X;
444
445 /* Initialize Hessian entries and work vector to zero */
446 PetscCall(MatAssembled(H, &assembled));
447 if (assembled) PetscCall(MatZeroEntries(H));
448
449 PetscCall(VecSet(user->s, zero));
450
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));
456
457 PetscCall(HessianProduct(ptr, user->s, user->y));
458
459 PetscCall(VecSetValues(user->s, 1, &j, &zero, INSERT_VALUES));
460 PetscCall(VecAssemblyBegin(user->s));
461 PetscCall(VecAssemblyEnd(user->s));
462
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 }
473
474 /* ------------------------------------------------------------------- */
475 /*
476 MatrixFreeHessian - Sets a pointer for use in computing Hessian-vector
477 products.
478
479 Input Parameters:
480 . tao - the Tao context
481 . X - the input vector
482 . ptr - optional user-defined context, as set by TaoSetHessian()
483
484 Output Parameters:
485 . H - Hessian matrix
486 . PrecH - optionally different preconditioning Hessian
487 */
MatrixFreeHessian(Tao tao,Vec X,Mat H,Mat PrecH,void * ptr)488 PetscErrorCode MatrixFreeHessian(Tao tao, Vec X, Mat H, Mat PrecH, void *ptr)
489 {
490 AppCtx *user = (AppCtx *)ptr;
491
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 }
497
498 /* ------------------------------------------------------------------- */
499 /*
500 HessianProductMat - Computes the matrix-vector product
501 y = mat*svec.
502
503 Input Parameters:
504 . mat - input matrix
505 . svec - input vector
506
507 Output Parameters:
508 . y - solution vector
509 */
HessianProductMat(Mat mat,Vec svec,Vec y)510 PetscErrorCode HessianProductMat(Mat mat, Vec svec, Vec y)
511 {
512 void *ptr;
513
514 PetscFunctionBeginUser;
515 PetscCall(MatShellGetContext(mat, &ptr));
516 PetscCall(HessianProduct(ptr, svec, y));
517 PetscFunctionReturn(PETSC_SUCCESS);
518 }
519
520 /* ------------------------------------------------------------------- */
521 /*
522 Hessian Product - Computes the matrix-vector product:
523 y = f''(x)*svec.
524
525 Input Parameters:
526 . ptr - pointer to the user-defined context
527 . svec - input vector
528
529 Output Parameters:
530 . y - product vector
531 */
HessianProduct(void * ptr,Vec svec,Vec y)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;
539
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);
547
548 /* Get pointers to vector data */
549 PetscCall(VecGetArrayRead(user->xvec, &x));
550 PetscCall(VecGetArrayRead(svec, &s));
551
552 /* Initialize product vector to zero */
553 PetscCall(VecSet(y, zero));
554
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 }
582
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));
613
614 /* Assemble vector */
615 PetscCall(VecAssemblyBegin(y));
616 PetscCall(VecAssemblyEnd(y));
617
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 }
624
625 /*TEST
626
627 build:
628 requires: !complex
629
630 test:
631 args: -tao_monitor -tao_type bntr -tao_view -tao_bnk_ksp_monitor_short
632
633 TEST*/
634