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
84 /* Initialize TAO,PETSc */
85 PetscFunctionBeginUser;
86 PetscCall(PetscInitialize(&argc, &argv, NULL, 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");
89
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));
96
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);
104
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));
109
110 /* Create TAO solver and set desired solution method */
111 PetscCall(TaoCreate(PETSC_COMM_SELF, &tao));
112 PetscCall(TaoSetType(tao, TAOLMVM));
113
114 /* Set solution vector with an initial guess */
115 PetscCall(FormInitialGuess(&user, x));
116 PetscCall(TaoSetSolution(tao, x));
117
118 /* Set routine for function and gradient evaluation */
119 PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user));
120
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, (PetscErrorCodeFn *)HessianProductMat));
126 PetscCall(MatSetOption(H, MAT_SYMMETRIC, PETSC_TRUE));
127
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 }
134
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 }
140
141 /* Check for any TAO command line options */
142 PetscCall(TaoSetFromOptions(tao));
143
144 /* SOLVE THE APPLICATION */
145 PetscCall(TaoSolve(tao));
146
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 }
165
166 PetscCall(TaoDestroy(&tao));
167 PetscCall(VecDestroy(&user.s));
168 PetscCall(VecDestroy(&user.y));
169 PetscCall(VecDestroy(&x));
170 PetscCall(MatDestroy(&H));
171
172 PetscCall(PetscFinalize());
173 return 0;
174 }
175
176 /* ------------------------------------------------------------------- */
177 /*
178 FormInitialGuess - Computes an initial approximation to the solution.
179
180 Input Parameters:
181 . user - user-defined application context
182 . X - vector
183
184 Output Parameters:
185 . X - vector
186 */
FormInitialGuess(AppCtx * user,Vec X)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;
192
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 }
207
208 /* ------------------------------------------------------------------- */
209 /*
210 FormFunctionGradient - Evaluates the function and corresponding gradient.
211
212 Input Parameters:
213 tao - the Tao context
214 X - the input vector
215 ptr - optional user-defined context, as set by TaoSetFunction()
216
217 Output Parameters:
218 f - the newly evaluated function
219 G - the newly evaluated gradient
220 */
FormFunctionGradient(Tao tao,Vec X,PetscReal * f,Vec G,void * ptr)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 }
228
229 /* ------------------------------------------------------------------- */
230 /*
231 FormFunction - Evaluates the function, f(X).
232
233 Input Parameters:
234 . tao - the Tao context
235 . X - the input vector
236 . ptr - optional user-defined context, as set by TaoSetFunction()
237
238 Output Parameters:
239 . f - the newly evaluated function
240 */
FormFunction(Tao tao,Vec X,PetscReal * f,void * ptr)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;
249
250 PetscFunctionBeginUser;
251 /* Get pointer to vector data */
252 PetscCall(VecGetArrayRead(X, &x));
253
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 }
270
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 }
287
288 /* Restore vector */
289 PetscCall(VecRestoreArrayRead(X, &x));
290
291 /* Scale the function */
292 area = p5 * hx * hy;
293 *f = area * (p5 * fquad + flin);
294
295 PetscCall(PetscLogFlops(24.0 * nx * ny));
296 PetscFunctionReturn(PETSC_SUCCESS);
297 }
298
299 /* ------------------------------------------------------------------- */
300 /*
301 FormGradient - Evaluates the gradient, G(X).
302
303 Input Parameters:
304 . tao - the Tao context
305 . X - input vector
306 . ptr - optional user-defined context
307
308 Output Parameters:
309 . G - vector containing the newly evaluated gradient
310 */
FormGradient(Tao tao,Vec X,Vec G,void * ptr)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;
320
321 PetscFunctionBeginUser;
322 /* Initialize gradient to zero */
323 PetscCall(VecSet(G, zero));
324
325 /* Get pointer to vector data */
326 PetscCall(VecGetArrayRead(X, &x));
327
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 }
357
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));
388
389 /* Assemble gradient vector */
390 PetscCall(VecAssemblyBegin(G));
391 PetscCall(VecAssemblyEnd(G));
392
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 }
399
400 /* ------------------------------------------------------------------- */
401 /*
402 FormHessian - Forms the Hessian matrix.
403
404 Input Parameters:
405 . tao - the Tao context
406 . X - the input vector
407 . ptr - optional user-defined context, as set by TaoSetHessian()
408
409 Output Parameters:
410 . H - Hessian matrix
411 . PrecH - optionally different preconditioning Hessian
412
413 Notes:
414 This routine is intended simply as an example of the interface
415 to a Hessian evaluation routine. Since this example compute the
416 Hessian a column at a time, it is not particularly efficient and
417 is not recommended.
418 */
FormHessian(Tao tao,Vec X,Mat H,Mat Hpre,void * ptr)419 PetscErrorCode FormHessian(Tao tao, Vec X, Mat H, Mat Hpre, void *ptr)
420 {
421 AppCtx *user = (AppCtx *)ptr;
422 PetscInt i, j, ndim = user->ndim;
423 PetscReal *y, zero = 0.0, one = 1.0;
424 PetscBool assembled;
425
426 PetscFunctionBeginUser;
427 user->xvec = X;
428
429 /* Initialize Hessian entries and work vector to zero */
430 PetscCall(MatAssembled(H, &assembled));
431 if (assembled) PetscCall(MatZeroEntries(H));
432
433 PetscCall(VecSet(user->s, zero));
434
435 /* Loop over matrix columns to compute entries of the Hessian */
436 for (j = 0; j < ndim; j++) {
437 PetscCall(VecSetValues(user->s, 1, &j, &one, INSERT_VALUES));
438 PetscCall(VecAssemblyBegin(user->s));
439 PetscCall(VecAssemblyEnd(user->s));
440
441 PetscCall(HessianProduct(ptr, user->s, user->y));
442
443 PetscCall(VecSetValues(user->s, 1, &j, &zero, INSERT_VALUES));
444 PetscCall(VecAssemblyBegin(user->s));
445 PetscCall(VecAssemblyEnd(user->s));
446
447 PetscCall(VecGetArray(user->y, &y));
448 for (i = 0; i < ndim; i++) {
449 if (y[i] != zero) PetscCall(MatSetValues(H, 1, &i, 1, &j, &y[i], ADD_VALUES));
450 }
451 PetscCall(VecRestoreArray(user->y, &y));
452 }
453 PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
454 PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
455 PetscFunctionReturn(PETSC_SUCCESS);
456 }
457
458 /* ------------------------------------------------------------------- */
459 /*
460 MatrixFreeHessian - Sets a pointer for use in computing Hessian-vector
461 products.
462
463 Input Parameters:
464 . tao - the Tao context
465 . X - the input vector
466 . ptr - optional user-defined context, as set by TaoSetHessian()
467
468 Output Parameters:
469 . H - Hessian matrix
470 . PrecH - optionally different preconditioning Hessian
471 */
MatrixFreeHessian(Tao tao,Vec X,Mat H,Mat PrecH,void * ptr)472 PetscErrorCode MatrixFreeHessian(Tao tao, Vec X, Mat H, Mat PrecH, void *ptr)
473 {
474 AppCtx *user = (AppCtx *)ptr;
475
476 /* Sets location of vector for use in computing matrix-vector products of the form H(X)*y */
477 PetscFunctionBeginUser;
478 user->xvec = X;
479 PetscFunctionReturn(PETSC_SUCCESS);
480 }
481
482 /* ------------------------------------------------------------------- */
483 /*
484 HessianProductMat - Computes the matrix-vector product
485 y = mat*svec.
486
487 Input Parameters:
488 . mat - input matrix
489 . svec - input vector
490
491 Output Parameters:
492 . y - solution vector
493 */
HessianProductMat(Mat mat,Vec svec,Vec y)494 PetscErrorCode HessianProductMat(Mat mat, Vec svec, Vec y)
495 {
496 void *ptr;
497
498 PetscFunctionBeginUser;
499 PetscCall(MatShellGetContext(mat, &ptr));
500 PetscCall(HessianProduct(ptr, svec, y));
501 PetscFunctionReturn(PETSC_SUCCESS);
502 }
503
504 /* ------------------------------------------------------------------- */
505 /*
506 Hessian Product - Computes the matrix-vector product:
507 y = f''(x)*svec.
508
509 Input Parameters:
510 . ptr - pointer to the user-defined context
511 . svec - input vector
512
513 Output Parameters:
514 . y - product vector
515 */
HessianProduct(void * ptr,Vec svec,Vec y)516 PetscErrorCode HessianProduct(void *ptr, Vec svec, Vec y)
517 {
518 AppCtx *user = (AppCtx *)ptr;
519 PetscReal p5 = 0.5, zero = 0.0, one = 1.0, hx, hy, val, area;
520 const PetscScalar *x, *s;
521 PetscReal v, vb, vl, vr, vt, hxhx, hyhy;
522 PetscInt nx, ny, i, j, k, ind;
523
524 PetscFunctionBeginUser;
525 nx = user->mx;
526 ny = user->my;
527 hx = user->hx;
528 hy = user->hy;
529 hxhx = one / (hx * hx);
530 hyhy = one / (hy * hy);
531
532 /* Get pointers to vector data */
533 PetscCall(VecGetArrayRead(user->xvec, &x));
534 PetscCall(VecGetArrayRead(svec, &s));
535
536 /* Initialize product vector to zero */
537 PetscCall(VecSet(y, zero));
538
539 /* Compute f''(x)*s over the lower triangular elements */
540 for (j = -1; j < ny; j++) {
541 for (i = -1; i < nx; i++) {
542 k = nx * j + i;
543 v = zero;
544 vr = zero;
545 vt = zero;
546 if (i != -1 && j != -1) v = s[k];
547 if (i != nx - 1 && j != -1) {
548 vr = s[k + 1];
549 ind = k + 1;
550 val = hxhx * (vr - v);
551 PetscCall(VecSetValues(y, 1, &ind, &val, ADD_VALUES));
552 }
553 if (i != -1 && j != ny - 1) {
554 vt = s[k + nx];
555 ind = k + nx;
556 val = hyhy * (vt - v);
557 PetscCall(VecSetValues(y, 1, &ind, &val, ADD_VALUES));
558 }
559 if (i != -1 && j != -1) {
560 ind = k;
561 val = hxhx * (v - vr) + hyhy * (v - vt);
562 PetscCall(VecSetValues(y, 1, &ind, &val, ADD_VALUES));
563 }
564 }
565 }
566
567 /* Compute f''(x)*s over the upper triangular elements */
568 for (j = 0; j <= ny; j++) {
569 for (i = 0; i <= nx; i++) {
570 k = nx * j + i;
571 v = zero;
572 vl = zero;
573 vb = zero;
574 if (i != nx && j != ny) v = s[k];
575 if (i != nx && j != 0) {
576 vb = s[k - nx];
577 ind = k - nx;
578 val = hyhy * (vb - v);
579 PetscCall(VecSetValues(y, 1, &ind, &val, ADD_VALUES));
580 }
581 if (i != 0 && j != ny) {
582 vl = s[k - 1];
583 ind = k - 1;
584 val = hxhx * (vl - v);
585 PetscCall(VecSetValues(y, 1, &ind, &val, ADD_VALUES));
586 }
587 if (i != nx && j != ny) {
588 ind = k;
589 val = hxhx * (v - vl) + hyhy * (v - vb);
590 PetscCall(VecSetValues(y, 1, &ind, &val, ADD_VALUES));
591 }
592 }
593 }
594 /* Restore vector data */
595 PetscCall(VecRestoreArrayRead(svec, &s));
596 PetscCall(VecRestoreArrayRead(user->xvec, &x));
597
598 /* Assemble vector */
599 PetscCall(VecAssemblyBegin(y));
600 PetscCall(VecAssemblyEnd(y));
601
602 /* Scale resulting vector by area */
603 area = p5 * hx * hy;
604 PetscCall(VecScale(y, area));
605 PetscCall(PetscLogFlops(18.0 * nx * ny));
606 PetscFunctionReturn(PETSC_SUCCESS);
607 }
608
609 /*TEST
610
611 build:
612 requires: !complex
613
614 test:
615 suffix: 1
616 args: -tao_monitor_short -tao_type ntl -tao_gatol 1.e-4
617
618 test:
619 suffix: 2
620 args: -tao_monitor_short -tao_type ntr -tao_gatol 1.e-4
621
622 test:
623 suffix: 3
624 args: -tao_monitor_short -tao_type bntr -tao_gatol 1.e-4 -my_tao_mf -tao_test_hessian
625
626 test:
627 suffix: 4
628 args: -tao_monitor_short -tao_gatol 1e-3 -tao_type bqnls
629
630 test:
631 suffix: 5
632 args: -tao_monitor_short -tao_gatol 1e-3 -tao_type blmvm
633
634 test:
635 suffix: 6
636 args: -tao_monitor_short -tao_gatol 1e-3 -tao_type bqnktr -tao_bqnk_mat_type lmvmsr1
637
638 test:
639 suffix: snes
640 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
641
642 test:
643 suffix: snes_2
644 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
645
646 test:
647 suffix: snes_3
648 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
649
650 TEST*/
651