1 /* Program usage: mpiexec -n 1 minsurf1 [-help] [all TAO options] */
2
3 /* Include "petsctao.h" so we can use TAO solvers. */
4 #include <petsctao.h>
5
6 static char help[] = "This example demonstrates use of the TAO package to \n\
7 solve an unconstrained minimization problem. This example is based on a \n\
8 problem from the MINPACK-2 test suite. Given a rectangular 2-D domain and \n\
9 boundary values along the edges of the domain, the objective is to find the\n\
10 surface with the minimal area that satisfies the boundary conditions.\n\
11 The command line options are:\n\
12 -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
13 -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
14 -start <st>, where <st> =0 for zero vector, and an average of the boundary conditions otherwise \n\n";
15
16 /*
17 User-defined application context - contains data needed by the
18 application-provided call-back routines, FormFunctionGradient()
19 and FormHessian().
20 */
21 typedef struct {
22 PetscInt mx, my; /* discretization in x, y directions */
23 PetscReal *bottom, *top, *left, *right; /* boundary values */
24 Mat H;
25 } AppCtx;
26
27 /* -------- User-defined Routines --------- */
28
29 static PetscErrorCode MSA_BoundaryConditions(AppCtx *);
30 static PetscErrorCode MSA_InitialPoint(AppCtx *, Vec);
31 static PetscErrorCode QuadraticH(AppCtx *, Vec, Mat);
32 PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);
33 PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void *);
34
main(int argc,char ** argv)35 int main(int argc, char **argv)
36 {
37 PetscInt N; /* Size of vector */
38 PetscMPIInt size; /* Number of processors */
39 Vec x; /* solution, gradient vectors */
40 KSP ksp; /* PETSc Krylov subspace method */
41 PetscBool flg; /* A return value when checking for user options */
42 Tao tao; /* Tao solver context */
43 AppCtx user; /* user-defined work context */
44
45 /* Initialize TAO,PETSc */
46 PetscFunctionBeginUser;
47 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
48
49 PetscCallMPI(MPI_Comm_size(MPI_COMM_WORLD, &size));
50 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Incorrect number of processors");
51
52 /* Specify default dimension of the problem */
53 user.mx = 4;
54 user.my = 4;
55
56 /* Check for any command line arguments that override defaults */
57 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &user.mx, &flg));
58 PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &user.my, &flg));
59
60 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n---- Minimum Surface Area Problem -----\n"));
61 PetscCall(PetscPrintf(PETSC_COMM_SELF, "mx: %" PetscInt_FMT " my: %" PetscInt_FMT " \n\n", user.mx, user.my));
62
63 /* Calculate any derived values from parameters */
64 N = user.mx * user.my;
65
66 /* Create TAO solver and set desired solution method */
67 PetscCall(TaoCreate(PETSC_COMM_SELF, &tao));
68 PetscCall(TaoSetType(tao, TAOLMVM));
69
70 /* Initialize minsurf application data structure for use in the function evaluations */
71 PetscCall(MSA_BoundaryConditions(&user)); /* Application specific routine */
72
73 /*
74 Create a vector to hold the variables. Compute an initial solution.
75 Set this vector, which will be used by TAO.
76 */
77 PetscCall(VecCreateSeq(PETSC_COMM_SELF, N, &x));
78 PetscCall(MSA_InitialPoint(&user, x)); /* Application specific routine */
79 PetscCall(TaoSetSolution(tao, x)); /* A TAO routine */
80
81 /* Provide TAO routines for function, gradient, and Hessian evaluation */
82 PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user));
83
84 /* Create a matrix data structure to store the Hessian. This structure will be used by TAO */
85 PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, N, N, 7, NULL, &user.H));
86 PetscCall(MatSetOption(user.H, MAT_SYMMETRIC, PETSC_TRUE));
87 PetscCall(TaoSetHessian(tao, user.H, user.H, FormHessian, (void *)&user));
88
89 /* Check for any TAO command line options */
90 PetscCall(TaoSetFromOptions(tao));
91
92 /* Limit the number of iterations in the KSP linear solver */
93 PetscCall(TaoGetKSP(tao, &ksp));
94 if (ksp) PetscCall(KSPSetTolerances(ksp, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT, user.mx * user.my));
95
96 /* SOLVE THE APPLICATION */
97 PetscCall(TaoSolve(tao));
98
99 PetscCall(TaoDestroy(&tao));
100 PetscCall(VecDestroy(&x));
101 PetscCall(MatDestroy(&user.H));
102 PetscCall(PetscFree(user.bottom));
103 PetscCall(PetscFree(user.top));
104 PetscCall(PetscFree(user.left));
105 PetscCall(PetscFree(user.right));
106
107 PetscCall(PetscFinalize());
108 return 0;
109 }
110
111 /* -------------------------------------------------------------------- */
112
113 /* FormFunctionGradient - Evaluates function and corresponding gradient.
114
115 Input Parameters:
116 . tao - the Tao context
117 . X - input vector
118 . userCtx - optional user-defined context, as set by TaoSetFunctionGradient()
119
120 Output Parameters:
121 . fcn - the newly evaluated function
122 . G - vector containing the newly evaluated gradient
123 */
FormFunctionGradient(Tao tao,Vec X,PetscReal * fcn,Vec G,void * userCtx)124 PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn, Vec G, void *userCtx)
125 {
126 AppCtx *user = (AppCtx *)userCtx;
127 PetscInt i, j, row;
128 PetscInt mx = user->mx, my = user->my;
129 PetscReal rhx = mx + 1, rhy = my + 1;
130 PetscReal hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), hydhx = hy / hx, hxdhy = hx / hy, area = 0.5 * hx * hy, ft = 0;
131 PetscReal f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb;
132 PetscReal df1dxc, df2dxc, df3dxc, df4dxc, df5dxc, df6dxc;
133 PetscReal zero = 0.0;
134 PetscScalar *g;
135 const PetscScalar *x;
136
137 PetscFunctionBeginUser;
138 PetscCall(VecSet(G, zero));
139
140 PetscCall(VecGetArrayRead(X, &x));
141 PetscCall(VecGetArray(G, &g));
142
143 /* Compute function over the locally owned part of the mesh */
144 for (j = 0; j < my; j++) {
145 for (i = 0; i < mx; i++) {
146 row = (j)*mx + (i);
147 xc = x[row];
148 xlt = xrb = xl = xr = xb = xt = xc;
149 if (i == 0) { /* left side */
150 xl = user->left[j + 1];
151 xlt = user->left[j + 2];
152 } else {
153 xl = x[row - 1];
154 }
155
156 if (j == 0) { /* bottom side */
157 xb = user->bottom[i + 1];
158 xrb = user->bottom[i + 2];
159 } else {
160 xb = x[row - mx];
161 }
162
163 if (i + 1 == mx) { /* right side */
164 xr = user->right[j + 1];
165 xrb = user->right[j];
166 } else {
167 xr = x[row + 1];
168 }
169
170 if (j + 1 == 0 + my) { /* top side */
171 xt = user->top[i + 1];
172 xlt = user->top[i];
173 } else {
174 xt = x[row + mx];
175 }
176
177 if (i > 0 && j + 1 < my) xlt = x[row - 1 + mx];
178 if (j > 0 && i + 1 < mx) xrb = x[row + 1 - mx];
179
180 d1 = (xc - xl);
181 d2 = (xc - xr);
182 d3 = (xc - xt);
183 d4 = (xc - xb);
184 d5 = (xr - xrb);
185 d6 = (xrb - xb);
186 d7 = (xlt - xl);
187 d8 = (xt - xlt);
188
189 df1dxc = d1 * hydhx;
190 df2dxc = (d1 * hydhx + d4 * hxdhy);
191 df3dxc = d3 * hxdhy;
192 df4dxc = (d2 * hydhx + d3 * hxdhy);
193 df5dxc = d2 * hydhx;
194 df6dxc = d4 * hxdhy;
195
196 d1 *= rhx;
197 d2 *= rhx;
198 d3 *= rhy;
199 d4 *= rhy;
200 d5 *= rhy;
201 d6 *= rhx;
202 d7 *= rhy;
203 d8 *= rhx;
204
205 f1 = PetscSqrtScalar(1.0 + d1 * d1 + d7 * d7);
206 f2 = PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
207 f3 = PetscSqrtScalar(1.0 + d3 * d3 + d8 * d8);
208 f4 = PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
209 f5 = PetscSqrtScalar(1.0 + d2 * d2 + d5 * d5);
210 f6 = PetscSqrtScalar(1.0 + d4 * d4 + d6 * d6);
211
212 ft = ft + (f2 + f4);
213
214 df1dxc /= f1;
215 df2dxc /= f2;
216 df3dxc /= f3;
217 df4dxc /= f4;
218 df5dxc /= f5;
219 df6dxc /= f6;
220
221 g[row] = (df1dxc + df2dxc + df3dxc + df4dxc + df5dxc + df6dxc) / 2.0;
222 }
223 }
224
225 for (j = 0; j < my; j++) { /* left side */
226 d3 = (user->left[j + 1] - user->left[j + 2]) * rhy;
227 d2 = (user->left[j + 1] - x[j * mx]) * rhx;
228 ft = ft + PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
229 }
230
231 for (i = 0; i < mx; i++) { /* bottom */
232 d2 = (user->bottom[i + 1] - user->bottom[i + 2]) * rhx;
233 d3 = (user->bottom[i + 1] - x[i]) * rhy;
234 ft = ft + PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
235 }
236
237 for (j = 0; j < my; j++) { /* right side */
238 d1 = (x[(j + 1) * mx - 1] - user->right[j + 1]) * rhx;
239 d4 = (user->right[j] - user->right[j + 1]) * rhy;
240 ft = ft + PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
241 }
242
243 for (i = 0; i < mx; i++) { /* top side */
244 d1 = (x[(my - 1) * mx + i] - user->top[i + 1]) * rhy;
245 d4 = (user->top[i + 1] - user->top[i]) * rhx;
246 ft = ft + PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
247 }
248
249 /* Bottom left corner */
250 d1 = (user->left[0] - user->left[1]) * rhy;
251 d2 = (user->bottom[0] - user->bottom[1]) * rhx;
252 ft += PetscSqrtScalar(1.0 + d1 * d1 + d2 * d2);
253
254 /* Top right corner */
255 d1 = (user->right[my + 1] - user->right[my]) * rhy;
256 d2 = (user->top[mx + 1] - user->top[mx]) * rhx;
257 ft += PetscSqrtScalar(1.0 + d1 * d1 + d2 * d2);
258
259 (*fcn) = ft * area;
260
261 /* Restore vectors */
262 PetscCall(VecRestoreArrayRead(X, &x));
263 PetscCall(VecRestoreArray(G, &g));
264 PetscCall(PetscLogFlops(67.0 * mx * my));
265 PetscFunctionReturn(PETSC_SUCCESS);
266 }
267
268 /* ------------------------------------------------------------------- */
269 /*
270 FormHessian - Evaluates the Hessian matrix.
271
272 Input Parameters:
273 . tao - the Tao context
274 . x - input vector
275 . ptr - optional user-defined context, as set by TaoSetHessian()
276
277 Output Parameters:
278 . H - Hessian matrix
279 . Hpre - optionally different matrix used to compute the preconditioner
280
281 */
FormHessian(Tao tao,Vec X,Mat H,Mat Hpre,void * ptr)282 PetscErrorCode FormHessian(Tao tao, Vec X, Mat H, Mat Hpre, void *ptr)
283 {
284 AppCtx *user = (AppCtx *)ptr;
285
286 PetscFunctionBeginUser;
287 /* Evaluate the Hessian entries*/
288 PetscCall(QuadraticH(user, X, H));
289 PetscFunctionReturn(PETSC_SUCCESS);
290 }
291
292 /* ------------------------------------------------------------------- */
293 /*
294 QuadraticH - Evaluates the Hessian matrix.
295
296 Input Parameters:
297 . user - user-defined context, as set by TaoSetHessian()
298 . X - input vector
299
300 Output Parameter:
301 . H - Hessian matrix
302 */
QuadraticH(AppCtx * user,Vec X,Mat Hessian)303 PetscErrorCode QuadraticH(AppCtx *user, Vec X, Mat Hessian)
304 {
305 PetscInt i, j, k, row;
306 PetscInt mx = user->mx, my = user->my;
307 PetscInt col[7];
308 PetscReal hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), hydhx = hy / hx, hxdhy = hx / hy;
309 PetscReal rhx = mx + 1, rhy = my + 1;
310 PetscReal f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb;
311 PetscReal hl, hr, ht, hb, hc, htl, hbr;
312 const PetscScalar *x;
313 PetscReal v[7];
314
315 PetscFunctionBeginUser;
316 /* Get pointers to vector data */
317 PetscCall(VecGetArrayRead(X, &x));
318
319 /* Initialize matrix entries to zero */
320 PetscCall(MatZeroEntries(Hessian));
321
322 /* Set various matrix options */
323 PetscCall(MatSetOption(Hessian, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE));
324
325 /* Compute Hessian over the locally owned part of the mesh */
326 for (i = 0; i < mx; i++) {
327 for (j = 0; j < my; j++) {
328 row = (j)*mx + (i);
329
330 xc = x[row];
331 xlt = xrb = xl = xr = xb = xt = xc;
332
333 /* Left side */
334 if (i == 0) {
335 xl = user->left[j + 1];
336 xlt = user->left[j + 2];
337 } else {
338 xl = x[row - 1];
339 }
340
341 if (j == 0) {
342 xb = user->bottom[i + 1];
343 xrb = user->bottom[i + 2];
344 } else {
345 xb = x[row - mx];
346 }
347
348 if (i + 1 == mx) {
349 xr = user->right[j + 1];
350 xrb = user->right[j];
351 } else {
352 xr = x[row + 1];
353 }
354
355 if (j + 1 == my) {
356 xt = user->top[i + 1];
357 xlt = user->top[i];
358 } else {
359 xt = x[row + mx];
360 }
361
362 if (i > 0 && j + 1 < my) xlt = x[row - 1 + mx];
363 if (j > 0 && i + 1 < mx) xrb = x[row + 1 - mx];
364
365 d1 = (xc - xl) * rhx;
366 d2 = (xc - xr) * rhx;
367 d3 = (xc - xt) * rhy;
368 d4 = (xc - xb) * rhy;
369 d5 = (xrb - xr) * rhy;
370 d6 = (xrb - xb) * rhx;
371 d7 = (xlt - xl) * rhy;
372 d8 = (xlt - xt) * rhx;
373
374 f1 = PetscSqrtScalar(1.0 + d1 * d1 + d7 * d7);
375 f2 = PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
376 f3 = PetscSqrtScalar(1.0 + d3 * d3 + d8 * d8);
377 f4 = PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
378 f5 = PetscSqrtScalar(1.0 + d2 * d2 + d5 * d5);
379 f6 = PetscSqrtScalar(1.0 + d4 * d4 + d6 * d6);
380
381 hl = (-hydhx * (1.0 + d7 * d7) + d1 * d7) / (f1 * f1 * f1) + (-hydhx * (1.0 + d4 * d4) + d1 * d4) / (f2 * f2 * f2);
382 hr = (-hydhx * (1.0 + d5 * d5) + d2 * d5) / (f5 * f5 * f5) + (-hydhx * (1.0 + d3 * d3) + d2 * d3) / (f4 * f4 * f4);
383 ht = (-hxdhy * (1.0 + d8 * d8) + d3 * d8) / (f3 * f3 * f3) + (-hxdhy * (1.0 + d2 * d2) + d2 * d3) / (f4 * f4 * f4);
384 hb = (-hxdhy * (1.0 + d6 * d6) + d4 * d6) / (f6 * f6 * f6) + (-hxdhy * (1.0 + d1 * d1) + d1 * d4) / (f2 * f2 * f2);
385
386 hbr = -d2 * d5 / (f5 * f5 * f5) - d4 * d6 / (f6 * f6 * f6);
387 htl = -d1 * d7 / (f1 * f1 * f1) - d3 * d8 / (f3 * f3 * f3);
388
389 hc = hydhx * (1.0 + d7 * d7) / (f1 * f1 * f1) + hxdhy * (1.0 + d8 * d8) / (f3 * f3 * f3) + hydhx * (1.0 + d5 * d5) / (f5 * f5 * f5) + hxdhy * (1.0 + d6 * d6) / (f6 * f6 * f6) + (hxdhy * (1.0 + d1 * d1) + hydhx * (1.0 + d4 * d4) - 2 * d1 * d4) / (f2 * f2 * f2) + (hxdhy * (1.0 + d2 * d2) + hydhx * (1.0 + d3 * d3) - 2 * d2 * d3) / (f4 * f4 * f4);
390
391 hl *= 0.5;
392 hr *= 0.5;
393 ht *= 0.5;
394 hb *= 0.5;
395 hbr *= 0.5;
396 htl *= 0.5;
397 hc *= 0.5;
398
399 k = 0;
400 if (j > 0) {
401 v[k] = hb;
402 col[k] = row - mx;
403 k++;
404 }
405
406 if (j > 0 && i < mx - 1) {
407 v[k] = hbr;
408 col[k] = row - mx + 1;
409 k++;
410 }
411
412 if (i > 0) {
413 v[k] = hl;
414 col[k] = row - 1;
415 k++;
416 }
417
418 v[k] = hc;
419 col[k] = row;
420 k++;
421
422 if (i < mx - 1) {
423 v[k] = hr;
424 col[k] = row + 1;
425 k++;
426 }
427
428 if (i > 0 && j < my - 1) {
429 v[k] = htl;
430 col[k] = row + mx - 1;
431 k++;
432 }
433
434 if (j < my - 1) {
435 v[k] = ht;
436 col[k] = row + mx;
437 k++;
438 }
439
440 /*
441 Set matrix values using local numbering, which was defined
442 earlier, in the main routine.
443 */
444 PetscCall(MatSetValues(Hessian, 1, &row, k, col, v, INSERT_VALUES));
445 }
446 }
447
448 /* Restore vectors */
449 PetscCall(VecRestoreArrayRead(X, &x));
450
451 /* Assemble the matrix */
452 PetscCall(MatAssemblyBegin(Hessian, MAT_FINAL_ASSEMBLY));
453 PetscCall(MatAssemblyEnd(Hessian, MAT_FINAL_ASSEMBLY));
454
455 PetscCall(PetscLogFlops(199.0 * mx * my));
456 PetscFunctionReturn(PETSC_SUCCESS);
457 }
458
459 /* ------------------------------------------------------------------- */
460 /*
461 MSA_BoundaryConditions - Calculates the boundary conditions for
462 the region.
463
464 Input Parameter:
465 . user - user-defined application context
466
467 Output Parameter:
468 . user - user-defined application context
469 */
MSA_BoundaryConditions(AppCtx * user)470 static PetscErrorCode MSA_BoundaryConditions(AppCtx *user)
471 {
472 PetscInt i, j, k, limit = 0;
473 PetscInt maxits = 5;
474 PetscInt mx = user->mx, my = user->my;
475 PetscInt bsize = 0, lsize = 0, tsize = 0, rsize = 0;
476 PetscReal one = 1.0, two = 2.0, three = 3.0, tol = 1e-10;
477 PetscReal fnorm, det, hx, hy, xt = 0, yt = 0;
478 PetscReal u1, u2, nf1, nf2, njac11, njac12, njac21, njac22;
479 PetscReal b = -0.5, t = 0.5, l = -0.5, r = 0.5;
480 PetscReal *boundary;
481
482 PetscFunctionBeginUser;
483 bsize = mx + 2;
484 lsize = my + 2;
485 rsize = my + 2;
486 tsize = mx + 2;
487
488 PetscCall(PetscMalloc1(bsize, &user->bottom));
489 PetscCall(PetscMalloc1(tsize, &user->top));
490 PetscCall(PetscMalloc1(lsize, &user->left));
491 PetscCall(PetscMalloc1(rsize, &user->right));
492
493 hx = (r - l) / (mx + 1);
494 hy = (t - b) / (my + 1);
495
496 for (j = 0; j < 4; j++) {
497 if (j == 0) {
498 yt = b;
499 xt = l;
500 limit = bsize;
501 boundary = user->bottom;
502 } else if (j == 1) {
503 yt = t;
504 xt = l;
505 limit = tsize;
506 boundary = user->top;
507 } else if (j == 2) {
508 yt = b;
509 xt = l;
510 limit = lsize;
511 boundary = user->left;
512 } else { /* if (j==3) */
513 yt = b;
514 xt = r;
515 limit = rsize;
516 boundary = user->right;
517 }
518
519 for (i = 0; i < limit; i++) {
520 u1 = xt;
521 u2 = -yt;
522 for (k = 0; k < maxits; k++) {
523 nf1 = u1 + u1 * u2 * u2 - u1 * u1 * u1 / three - xt;
524 nf2 = -u2 - u1 * u1 * u2 + u2 * u2 * u2 / three - yt;
525 fnorm = PetscSqrtScalar(nf1 * nf1 + nf2 * nf2);
526 if (fnorm <= tol) break;
527 njac11 = one + u2 * u2 - u1 * u1;
528 njac12 = two * u1 * u2;
529 njac21 = -two * u1 * u2;
530 njac22 = -one - u1 * u1 + u2 * u2;
531 det = njac11 * njac22 - njac21 * njac12;
532 u1 = u1 - (njac22 * nf1 - njac12 * nf2) / det;
533 u2 = u2 - (njac11 * nf2 - njac21 * nf1) / det;
534 }
535
536 boundary[i] = u1 * u1 - u2 * u2;
537 if (j == 0 || j == 1) {
538 xt = xt + hx;
539 } else { /* if (j==2 || j==3) */
540 yt = yt + hy;
541 }
542 }
543 }
544 PetscFunctionReturn(PETSC_SUCCESS);
545 }
546
547 /* ------------------------------------------------------------------- */
548 /*
549 MSA_InitialPoint - Calculates the initial guess in one of three ways.
550
551 Input Parameters:
552 . user - user-defined application context
553 . X - vector for initial guess
554
555 Output Parameters:
556 . X - newly computed initial guess
557 */
MSA_InitialPoint(AppCtx * user,Vec X)558 static PetscErrorCode MSA_InitialPoint(AppCtx *user, Vec X)
559 {
560 PetscInt start = -1, i, j;
561 PetscReal zero = 0.0;
562 PetscBool flg;
563
564 PetscFunctionBeginUser;
565 PetscCall(VecSet(X, zero));
566 PetscCall(PetscOptionsGetInt(NULL, NULL, "-start", &start, &flg));
567
568 if (flg && start == 0) { /* The zero vector is reasonable */
569 PetscCall(VecSet(X, zero));
570 } else { /* Take an average of the boundary conditions */
571 PetscInt row;
572 PetscInt mx = user->mx, my = user->my;
573 PetscReal *x;
574
575 /* Get pointers to vector data */
576 PetscCall(VecGetArray(X, &x));
577 /* Perform local computations */
578 for (j = 0; j < my; j++) {
579 for (i = 0; i < mx; i++) {
580 row = (j)*mx + (i);
581 x[row] = (((j + 1) * user->bottom[i + 1] + (my - j + 1) * user->top[i + 1]) / (my + 2) + ((i + 1) * user->left[j + 1] + (mx - i + 1) * user->right[j + 1]) / (mx + 2)) / 2.0;
582 }
583 }
584 /* Restore vectors */
585 PetscCall(VecRestoreArray(X, &x));
586 }
587 PetscFunctionReturn(PETSC_SUCCESS);
588 }
589
590 /*TEST
591
592 build:
593 requires: !complex
594
595 test:
596 args: -tao_monitor_short -tao_type nls -mx 10 -my 8 -tao_gatol 1.e-4
597 requires: !single
598
599 test:
600 suffix: 2
601 args: -tao_monitor_short -tao_type bmrm -mx 10 -my 8 -tao_gatol 1.e-3
602 requires: !single
603
604 test:
605 suffix: 3
606 args: -tao_monitor_short -tao_type lmvm -mx 10 -my 8 -tao_gatol 1.e-3
607 requires: !single
608
609 test:
610 suffix: 4
611 args: -tao_monitor_short -tao_type bntr -mx 10 -my 8 -tao_gatol 1.e-4
612 requires: !single
613
614 test:
615 suffix: 4_ew
616 args: -tao_monitor_short -tao_type bntr -tao_ksp_ew -mx 10 -my 8 -tao_gatol 1.e-4
617 requires: !single
618
619 test:
620 suffix: 5
621 args: -tao_monitor_short -tao_type bntl -mx 10 -my 8 -tao_gatol 1.e-4
622 requires: !single
623
624 test:
625 suffix: 5_ew
626 args: -tao_monitor_short -tao_type bntl -tao_ksp_ew -mx 10 -my 8 -tao_gatol 1.e-4
627 requires: !single
628
629 test:
630 suffix: 6
631 args: -tao_monitor_short -tao_type bnls -mx 10 -my 8 -tao_gatol 1.e-4
632 requires: !single
633
634 test:
635 suffix: 6_ew
636 args: -tao_monitor_short -tao_type bnls -tao_ksp_ew -tao_bnk_ksp_ew_version 3 -mx 10 -my 8 -tao_gatol 1.e-4
637 requires: !single
638
639 test:
640 suffix: 7
641 args: -tao_monitor_short -tao_type bntr -mx 10 -my 8 -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4 -tao_bnk_cg_tao_monitor_short
642 requires: !single
643
644 test:
645 suffix: 8
646 args: -tao_monitor_short -tao_type bntl -mx 10 -my 8 -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4 -tao_bnk_cg_tao_monitor_short
647 requires: !single
648
649 test:
650 suffix: 9
651 args: -tao_monitor_short -tao_type bnls -mx 10 -my 8 -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4 -tao_bnk_cg_tao_monitor_short
652 requires: !single
653
654 test:
655 suffix: 10
656 args: -tao_monitor_short -tao_type bnls -mx 10 -my 8 -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4 -tao_mf_hessian -tao_bnk_cg_tao_monitor_short
657
658 test:
659 suffix: 11
660 args: -tao_monitor_short -tao_type bntr -mx 10 -my 8 -tao_gatol 1.e-4 -tao_mf_hessian
661 requires: !single
662
663 test:
664 suffix: 12
665 args: -tao_monitor_short -tao_type bntl -mx 10 -my 8 -tao_gatol 1.e-4 -tao_mf_hessian
666 requires: !single
667
668 TEST*/
669