Actual source code: plate2.c

```  1: #include <petscdmda.h>
2: #include <petsctao.h>

4: static char help[] = "This example demonstrates use of the TAO package to \n\
5: solve an unconstrained minimization problem.  This example is based on a \n\
6: problem from the MINPACK-2 test suite.  Given a rectangular 2-D domain, \n\
7: boundary values along the edges of the domain, and a plate represented by \n\
8: lower boundary conditions, the objective is to find the\n\
9: surface with the minimal area that satisfies the boundary conditions.\n\
10: The command line options are:\n\
11:   -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
12:   -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
13:   -bmx <bxg>, where <bxg> = number of grid points under plate in 1st direction\n\
14:   -bmy <byg>, where <byg> = number of grid points under plate in 2nd direction\n\
15:   -bheight <ht>, where <ht> = height of the plate\n\
16:   -start <st>, where <st> =0 for zero vector, <st> >0 for random start, and <st> <0 \n\
17:                for an average of the boundary conditions\n\n";

19: /*
20:    User-defined application context - contains data needed by the
22:    FormHessian().
23: */
24: typedef struct {
25:   /* problem parameters */
26:   PetscReal bheight;                  /* Height of plate under the surface */
27:   PetscInt  mx, my;                   /* discretization in x, y directions */
28:   PetscInt  bmx, bmy;                 /* Size of plate under the surface */
29:   Vec       Bottom, Top, Left, Right; /* boundary values */

31:   /* Working space */
32:   Vec localX, localV; /* ghosted local vector */
33:   DM  dm;             /* distributed array data structure */
34:   Mat H;
35: } AppCtx;

37: /* -------- User-defined Routines --------- */

39: static PetscErrorCode MSA_BoundaryConditions(AppCtx *);
40: static PetscErrorCode MSA_InitialPoint(AppCtx *, Vec);
41: static PetscErrorCode MSA_Plate(Vec, Vec, void *);
42: PetscErrorCode        FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);
43: PetscErrorCode        FormHessian(Tao, Vec, Mat, Mat, void *);

45: /* For testing matrix-free submatrices */
46: PetscErrorCode MatrixFreeHessian(Tao, Vec, Mat, Mat, void *);
47: PetscErrorCode MyMatMult(Mat, Vec, Vec);

49: int main(int argc, char **argv)
50: {
51:   PetscInt               Nx, Ny;    /* number of processors in x- and y- directions */
52:   PetscInt               m, N;      /* number of local and global elements in vectors */
53:   Vec                    x, xl, xu; /* solution vector  and bounds*/
54:   PetscBool              flg;       /* A return variable when checking for user options */
55:   Tao                    tao;       /* Tao solver context */
56:   ISLocalToGlobalMapping isltog;    /* local-to-global mapping object */
57:   Mat                    H_shell;   /* to test matrix-free submatrices */
58:   AppCtx                 user;      /* user-defined work context */

60:   /* Initialize PETSc, TAO */
61:   PetscFunctionBeginUser;
62:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));

64:   /* Specify default dimension of the problem */
65:   user.mx      = 10;
66:   user.my      = 10;
67:   user.bheight = 0.1;

69:   /* Check for any command line arguments that override defaults */
70:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &user.mx, &flg));
71:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &user.my, &flg));
72:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-bheight", &user.bheight, &flg));

74:   user.bmx = user.mx / 2;
75:   user.bmy = user.my / 2;
76:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-bmx", &user.bmx, &flg));
77:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-bmy", &user.bmy, &flg));

79:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n---- Minimum Surface Area With Plate Problem -----\n"));
80:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "mx:%" PetscInt_FMT ", my:%" PetscInt_FMT ", bmx:%" PetscInt_FMT ", bmy:%" PetscInt_FMT ", height:%g\n", user.mx, user.my, user.bmx, user.bmy, (double)user.bheight));

82:   /* Calculate any derived values from parameters */
83:   N = user.mx * user.my;

85:   /* Let Petsc determine the dimensions of the local vectors */
86:   Nx = PETSC_DECIDE;
87:   Ny = PETSC_DECIDE;

89:   /*
90:      A two dimensional distributed array will help define this problem,
91:      which derives from an elliptic PDE on two dimensional domain.  From
92:      the distributed array, Create the vectors.
93:   */
94:   PetscCall(DMDACreate2d(MPI_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, user.mx, user.my, Nx, Ny, 1, 1, NULL, NULL, &user.dm));
95:   PetscCall(DMSetFromOptions(user.dm));
96:   PetscCall(DMSetUp(user.dm));
97:   /*
98:      Extract global and local vectors from DM; The local vectors are
99:      used solely as work space for the evaluation of the function,
100:      gradient, and Hessian.  Duplicate for remaining vectors that are
101:      the same types.
102:   */
103:   PetscCall(DMCreateGlobalVector(user.dm, &x)); /* Solution */
104:   PetscCall(DMCreateLocalVector(user.dm, &user.localX));
105:   PetscCall(VecDuplicate(user.localX, &user.localV));

107:   PetscCall(VecDuplicate(x, &xl));
108:   PetscCall(VecDuplicate(x, &xu));

110:   /* The TAO code begins here */

112:   /*
113:      Create TAO solver and set desired solution method
114:      The method must either be TAOTRON or TAOBLMVM
115:      If TAOBLMVM is used, then hessian function is not called.
116:   */
117:   PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao));
118:   PetscCall(TaoSetType(tao, TAOBLMVM));

120:   /* Set initial solution guess; */
121:   PetscCall(MSA_BoundaryConditions(&user));
122:   PetscCall(MSA_InitialPoint(&user, x));
123:   PetscCall(TaoSetSolution(tao, x));

125:   /* Set routines for function, gradient and hessian evaluation */

128:   PetscCall(VecGetLocalSize(x, &m));
129:   PetscCall(DMCreateMatrix(user.dm, &user.H));
130:   PetscCall(MatSetOption(user.H, MAT_SYMMETRIC, PETSC_TRUE));

132:   PetscCall(DMGetLocalToGlobalMapping(user.dm, &isltog));
133:   PetscCall(MatSetLocalToGlobalMapping(user.H, isltog, isltog));
134:   PetscCall(PetscOptionsHasName(NULL, NULL, "-matrixfree", &flg));
135:   if (flg) {
136:     PetscCall(MatCreateShell(PETSC_COMM_WORLD, m, m, N, N, (void *)&user, &H_shell));
137:     PetscCall(MatShellSetOperation(H_shell, MATOP_MULT, (void (*)(void))MyMatMult));
138:     PetscCall(MatSetOption(H_shell, MAT_SYMMETRIC, PETSC_TRUE));
139:     PetscCall(TaoSetHessian(tao, H_shell, H_shell, MatrixFreeHessian, (void *)&user));
140:   } else {
141:     PetscCall(TaoSetHessian(tao, user.H, user.H, FormHessian, (void *)&user));
142:   }

144:   /* Set Variable bounds */
145:   PetscCall(MSA_Plate(xl, xu, (void *)&user));
146:   PetscCall(TaoSetVariableBounds(tao, xl, xu));

148:   /* Check for any tao command line options */
149:   PetscCall(TaoSetFromOptions(tao));

151:   /* SOLVE THE APPLICATION */
152:   PetscCall(TaoSolve(tao));

154:   PetscCall(TaoView(tao, PETSC_VIEWER_STDOUT_WORLD));

156:   /* Free TAO data structures */
157:   PetscCall(TaoDestroy(&tao));

159:   /* Free PETSc data structures */
160:   PetscCall(VecDestroy(&x));
161:   PetscCall(VecDestroy(&xl));
162:   PetscCall(VecDestroy(&xu));
163:   PetscCall(MatDestroy(&user.H));
164:   PetscCall(VecDestroy(&user.localX));
165:   PetscCall(VecDestroy(&user.localV));
166:   PetscCall(VecDestroy(&user.Bottom));
167:   PetscCall(VecDestroy(&user.Top));
168:   PetscCall(VecDestroy(&user.Left));
169:   PetscCall(VecDestroy(&user.Right));
170:   PetscCall(DMDestroy(&user.dm));
171:   if (flg) PetscCall(MatDestroy(&H_shell));
172:   PetscCall(PetscFinalize());
173:   return 0;
174: }

178:     Input Parameters:
179: .   tao     - the Tao context
180: .   X      - input vector
181: .   userCtx - optional user-defined context, as set by TaoSetObjectiveAndGradient()

183:     Output Parameters:
184: .   fcn     - the function value
185: .   G      - vector containing the newly evaluated gradient

187:    Notes:
188:    In this case, we discretize the domain and Create triangles. The
189:    surface of each triangle is planar, whose surface area can be easily
190:    computed.  The total surface area is found by sweeping through the grid
191:    and computing the surface area of the two triangles that have their
192:    right angle at the grid point.  The diagonal line segments on the
193:    grid that define the triangles run from top left to lower right.
194:    The numbering of points starts at the lower left and runs left to
195:    right, then bottom to top.
196: */
197: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn, Vec G, void *userCtx)
198: {
199:   AppCtx    *user = (AppCtx *)userCtx;
200:   PetscInt   i, j, row;
201:   PetscInt   mx = user->mx, my = user->my;
202:   PetscInt   xs, xm, gxs, gxm, ys, ym, gys, gym;
203:   PetscReal  ft   = 0;
204:   PetscReal  zero = 0.0;
205:   PetscReal  hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), hydhx = hy / hx, hxdhy = hx / hy, area = 0.5 * hx * hy;
206:   PetscReal  rhx = mx + 1, rhy = my + 1;
207:   PetscReal  f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb;
208:   PetscReal  df1dxc, df2dxc, df3dxc, df4dxc, df5dxc, df6dxc;
209:   PetscReal *g, *x, *left, *right, *bottom, *top;
210:   Vec        localX = user->localX, localG = user->localV;

212:   PetscFunctionBeginUser;
213:   /* Get local mesh boundaries */
214:   PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
215:   PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL));

217:   /* Scatter ghost points to local vector */
218:   PetscCall(DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, localX));
219:   PetscCall(DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, localX));

221:   /* Initialize vector to zero */
222:   PetscCall(VecSet(localG, zero));

224:   /* Get pointers to vector data */
225:   PetscCall(VecGetArray(localX, &x));
226:   PetscCall(VecGetArray(localG, &g));
227:   PetscCall(VecGetArray(user->Top, &top));
228:   PetscCall(VecGetArray(user->Bottom, &bottom));
229:   PetscCall(VecGetArray(user->Left, &left));
230:   PetscCall(VecGetArray(user->Right, &right));

232:   /* Compute function over the locally owned part of the mesh */
233:   for (j = ys; j < ys + ym; j++) {
234:     for (i = xs; i < xs + xm; i++) {
235:       row = (j - gys) * gxm + (i - gxs);

237:       xc  = x[row];
238:       xlt = xrb = xl = xr = xb = xt = xc;

240:       if (i == 0) { /* left side */
241:         xl  = left[j - ys + 1];
242:         xlt = left[j - ys + 2];
243:       } else {
244:         xl = x[row - 1];
245:       }

247:       if (j == 0) { /* bottom side */
248:         xb  = bottom[i - xs + 1];
249:         xrb = bottom[i - xs + 2];
250:       } else {
251:         xb = x[row - gxm];
252:       }

254:       if (i + 1 == gxs + gxm) { /* right side */
255:         xr  = right[j - ys + 1];
256:         xrb = right[j - ys];
257:       } else {
258:         xr = x[row + 1];
259:       }

261:       if (j + 1 == gys + gym) { /* top side */
262:         xt  = top[i - xs + 1];
263:         xlt = top[i - xs];
264:       } else {
265:         xt = x[row + gxm];
266:       }

268:       if (i > gxs && j + 1 < gys + gym) xlt = x[row - 1 + gxm];
269:       if (j > gys && i + 1 < gxs + gxm) xrb = x[row + 1 - gxm];

271:       d1 = (xc - xl);
272:       d2 = (xc - xr);
273:       d3 = (xc - xt);
274:       d4 = (xc - xb);
275:       d5 = (xr - xrb);
276:       d6 = (xrb - xb);
277:       d7 = (xlt - xl);
278:       d8 = (xt - xlt);

280:       df1dxc = d1 * hydhx;
281:       df2dxc = (d1 * hydhx + d4 * hxdhy);
282:       df3dxc = d3 * hxdhy;
283:       df4dxc = (d2 * hydhx + d3 * hxdhy);
284:       df5dxc = d2 * hydhx;
285:       df6dxc = d4 * hxdhy;

287:       d1 *= rhx;
288:       d2 *= rhx;
289:       d3 *= rhy;
290:       d4 *= rhy;
291:       d5 *= rhy;
292:       d6 *= rhx;
293:       d7 *= rhy;
294:       d8 *= rhx;

296:       f1 = PetscSqrtScalar(1.0 + d1 * d1 + d7 * d7);
297:       f2 = PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
298:       f3 = PetscSqrtScalar(1.0 + d3 * d3 + d8 * d8);
299:       f4 = PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
300:       f5 = PetscSqrtScalar(1.0 + d2 * d2 + d5 * d5);
301:       f6 = PetscSqrtScalar(1.0 + d4 * d4 + d6 * d6);

303:       ft = ft + (f2 + f4);

305:       df1dxc /= f1;
306:       df2dxc /= f2;
307:       df3dxc /= f3;
308:       df4dxc /= f4;
309:       df5dxc /= f5;
310:       df6dxc /= f6;

312:       g[row] = (df1dxc + df2dxc + df3dxc + df4dxc + df5dxc + df6dxc) * 0.5;
313:     }
314:   }

316:   /* Compute triangular areas along the border of the domain. */
317:   if (xs == 0) { /* left side */
318:     for (j = ys; j < ys + ym; j++) {
319:       d3 = (left[j - ys + 1] - left[j - ys + 2]) * rhy;
320:       d2 = (left[j - ys + 1] - x[(j - gys) * gxm]) * rhx;
321:       ft = ft + PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
322:     }
323:   }
324:   if (ys == 0) { /* bottom side */
325:     for (i = xs; i < xs + xm; i++) {
326:       d2 = (bottom[i + 1 - xs] - bottom[i - xs + 2]) * rhx;
327:       d3 = (bottom[i - xs + 1] - x[i - gxs]) * rhy;
328:       ft = ft + PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
329:     }
330:   }

332:   if (xs + xm == mx) { /* right side */
333:     for (j = ys; j < ys + ym; j++) {
334:       d1 = (x[(j + 1 - gys) * gxm - 1] - right[j - ys + 1]) * rhx;
335:       d4 = (right[j - ys] - right[j - ys + 1]) * rhy;
336:       ft = ft + PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
337:     }
338:   }
339:   if (ys + ym == my) { /* top side */
340:     for (i = xs; i < xs + xm; i++) {
341:       d1 = (x[(gym - 1) * gxm + i - gxs] - top[i - xs + 1]) * rhy;
342:       d4 = (top[i - xs + 1] - top[i - xs]) * rhx;
343:       ft = ft + PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
344:     }
345:   }

347:   if (ys == 0 && xs == 0) {
348:     d1 = (left[0] - left[1]) * rhy;
349:     d2 = (bottom[0] - bottom[1]) * rhx;
350:     ft += PetscSqrtScalar(1.0 + d1 * d1 + d2 * d2);
351:   }
352:   if (ys + ym == my && xs + xm == mx) {
353:     d1 = (right[ym + 1] - right[ym]) * rhy;
354:     d2 = (top[xm + 1] - top[xm]) * rhx;
355:     ft += PetscSqrtScalar(1.0 + d1 * d1 + d2 * d2);
356:   }

358:   ft = ft * area;
359:   PetscCall(MPIU_Allreduce(&ft, fcn, 1, MPIU_REAL, MPIU_SUM, MPI_COMM_WORLD));

361:   /* Restore vectors */
362:   PetscCall(VecRestoreArray(localX, &x));
363:   PetscCall(VecRestoreArray(localG, &g));
364:   PetscCall(VecRestoreArray(user->Left, &left));
365:   PetscCall(VecRestoreArray(user->Top, &top));
366:   PetscCall(VecRestoreArray(user->Bottom, &bottom));
367:   PetscCall(VecRestoreArray(user->Right, &right));

369:   /* Scatter values to global vector */
370:   PetscCall(DMLocalToGlobalBegin(user->dm, localG, INSERT_VALUES, G));
371:   PetscCall(DMLocalToGlobalEnd(user->dm, localG, INSERT_VALUES, G));

373:   PetscCall(PetscLogFlops(70.0 * xm * ym));
374:   PetscFunctionReturn(PETSC_SUCCESS);
375: }

377: /* ------------------------------------------------------------------- */
378: /*
379:    FormHessian - Evaluates Hessian matrix.

381:    Input Parameters:
382: .  tao  - the Tao context
383: .  x    - input vector
384: .  ptr  - optional user-defined context, as set by TaoSetHessian()

386:    Output Parameters:
387: .  A    - Hessian matrix
388: .  B    - optionally different preconditioning matrix

390:    Notes:
391:    Due to mesh point reordering with DMs, we must always work
392:    with the local mesh points, and then transform them to the new
393:    global numbering with the local-to-global mapping.  We cannot work
394:    directly with the global numbers for the original uniprocessor mesh!

396:    Two methods are available for imposing this transformation
397:    when setting matrix entries:
398:      (A) MatSetValuesLocal(), using the local ordering (including
399:          ghost points!)
400:          - Do the following two steps once, before calling TaoSolve()
401:            - Use DMGetISLocalToGlobalMapping() to extract the
402:              local-to-global map from the DM
403:            - Associate this map with the matrix by calling
404:              MatSetLocalToGlobalMapping()
405:          - Then set matrix entries using the local ordering
406:            by calling MatSetValuesLocal()
407:      (B) MatSetValues(), using the global ordering
408:          - Use DMGetGlobalIndices() to extract the local-to-global map
409:          - Then apply this map explicitly yourself
410:          - Set matrix entries using the global ordering by calling
411:            MatSetValues()
412:    Option (A) seems cleaner/easier in many cases, and is the procedure
413:    used in this example.
414: */
415: PetscErrorCode FormHessian(Tao tao, Vec X, Mat Hptr, Mat Hessian, void *ptr)
416: {
417:   AppCtx    *user = (AppCtx *)ptr;
418:   PetscInt   i, j, k, row;
419:   PetscInt   mx = user->mx, my = user->my;
420:   PetscInt   xs, xm, gxs, gxm, ys, ym, gys, gym, col[7];
421:   PetscReal  hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), hydhx = hy / hx, hxdhy = hx / hy;
422:   PetscReal  rhx = mx + 1, rhy = my + 1;
423:   PetscReal  f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb;
424:   PetscReal  hl, hr, ht, hb, hc, htl, hbr;
425:   PetscReal *x, *left, *right, *bottom, *top;
426:   PetscReal  v[7];
427:   Vec        localX = user->localX;
428:   PetscBool  assembled;

430:   PetscFunctionBeginUser;
431:   /* Set various matrix options */
432:   PetscCall(MatSetOption(Hessian, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE));

434:   /* Initialize matrix entries to zero */
435:   PetscCall(MatAssembled(Hessian, &assembled));
436:   if (assembled) PetscCall(MatZeroEntries(Hessian));

438:   /* Get local mesh boundaries */
439:   PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
440:   PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL));

442:   /* Scatter ghost points to local vector */
443:   PetscCall(DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, localX));
444:   PetscCall(DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, localX));

446:   /* Get pointers to vector data */
447:   PetscCall(VecGetArray(localX, &x));
448:   PetscCall(VecGetArray(user->Top, &top));
449:   PetscCall(VecGetArray(user->Bottom, &bottom));
450:   PetscCall(VecGetArray(user->Left, &left));
451:   PetscCall(VecGetArray(user->Right, &right));

453:   /* Compute Hessian over the locally owned part of the mesh */

455:   for (i = xs; i < xs + xm; i++) {
456:     for (j = ys; j < ys + ym; j++) {
457:       row = (j - gys) * gxm + (i - gxs);

459:       xc  = x[row];
460:       xlt = xrb = xl = xr = xb = xt = xc;

462:       /* Left side */
463:       if (i == gxs) {
464:         xl  = left[j - ys + 1];
465:         xlt = left[j - ys + 2];
466:       } else {
467:         xl = x[row - 1];
468:       }

470:       if (j == gys) {
471:         xb  = bottom[i - xs + 1];
472:         xrb = bottom[i - xs + 2];
473:       } else {
474:         xb = x[row - gxm];
475:       }

477:       if (i + 1 == gxs + gxm) {
478:         xr  = right[j - ys + 1];
479:         xrb = right[j - ys];
480:       } else {
481:         xr = x[row + 1];
482:       }

484:       if (j + 1 == gys + gym) {
485:         xt  = top[i - xs + 1];
486:         xlt = top[i - xs];
487:       } else {
488:         xt = x[row + gxm];
489:       }

491:       if (i > gxs && j + 1 < gys + gym) xlt = x[row - 1 + gxm];
492:       if (j > gys && i + 1 < gxs + gxm) xrb = x[row + 1 - gxm];

494:       d1 = (xc - xl) * rhx;
495:       d2 = (xc - xr) * rhx;
496:       d3 = (xc - xt) * rhy;
497:       d4 = (xc - xb) * rhy;
498:       d5 = (xrb - xr) * rhy;
499:       d6 = (xrb - xb) * rhx;
500:       d7 = (xlt - xl) * rhy;
501:       d8 = (xlt - xt) * rhx;

503:       f1 = PetscSqrtScalar(1.0 + d1 * d1 + d7 * d7);
504:       f2 = PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
505:       f3 = PetscSqrtScalar(1.0 + d3 * d3 + d8 * d8);
506:       f4 = PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
507:       f5 = PetscSqrtScalar(1.0 + d2 * d2 + d5 * d5);
508:       f6 = PetscSqrtScalar(1.0 + d4 * d4 + d6 * d6);

510:       hl = (-hydhx * (1.0 + d7 * d7) + d1 * d7) / (f1 * f1 * f1) + (-hydhx * (1.0 + d4 * d4) + d1 * d4) / (f2 * f2 * f2);
511:       hr = (-hydhx * (1.0 + d5 * d5) + d2 * d5) / (f5 * f5 * f5) + (-hydhx * (1.0 + d3 * d3) + d2 * d3) / (f4 * f4 * f4);
512:       ht = (-hxdhy * (1.0 + d8 * d8) + d3 * d8) / (f3 * f3 * f3) + (-hxdhy * (1.0 + d2 * d2) + d2 * d3) / (f4 * f4 * f4);
513:       hb = (-hxdhy * (1.0 + d6 * d6) + d4 * d6) / (f6 * f6 * f6) + (-hxdhy * (1.0 + d1 * d1) + d1 * d4) / (f2 * f2 * f2);

515:       hbr = -d2 * d5 / (f5 * f5 * f5) - d4 * d6 / (f6 * f6 * f6);
516:       htl = -d1 * d7 / (f1 * f1 * f1) - d3 * d8 / (f3 * f3 * f3);

518:       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);

520:       hl *= 0.5;
521:       hr *= 0.5;
522:       ht *= 0.5;
523:       hb *= 0.5;
524:       hbr *= 0.5;
525:       htl *= 0.5;
526:       hc *= 0.5;

528:       k = 0;
529:       if (j > 0) {
530:         v[k]   = hb;
531:         col[k] = row - gxm;
532:         k++;
533:       }

535:       if (j > 0 && i < mx - 1) {
536:         v[k]   = hbr;
537:         col[k] = row - gxm + 1;
538:         k++;
539:       }

541:       if (i > 0) {
542:         v[k]   = hl;
543:         col[k] = row - 1;
544:         k++;
545:       }

547:       v[k]   = hc;
548:       col[k] = row;
549:       k++;

551:       if (i < mx - 1) {
552:         v[k]   = hr;
553:         col[k] = row + 1;
554:         k++;
555:       }

557:       if (i > 0 && j < my - 1) {
558:         v[k]   = htl;
559:         col[k] = row + gxm - 1;
560:         k++;
561:       }

563:       if (j < my - 1) {
564:         v[k]   = ht;
565:         col[k] = row + gxm;
566:         k++;
567:       }

569:       /*
570:          Set matrix values using local numbering, which was defined
571:          earlier, in the main routine.
572:       */
573:       PetscCall(MatSetValuesLocal(Hessian, 1, &row, k, col, v, INSERT_VALUES));
574:     }
575:   }

577:   /* Restore vectors */
578:   PetscCall(VecRestoreArray(localX, &x));
579:   PetscCall(VecRestoreArray(user->Left, &left));
580:   PetscCall(VecRestoreArray(user->Top, &top));
581:   PetscCall(VecRestoreArray(user->Bottom, &bottom));
582:   PetscCall(VecRestoreArray(user->Right, &right));

584:   /* Assemble the matrix */
585:   PetscCall(MatAssemblyBegin(Hessian, MAT_FINAL_ASSEMBLY));
586:   PetscCall(MatAssemblyEnd(Hessian, MAT_FINAL_ASSEMBLY));

588:   PetscCall(PetscLogFlops(199.0 * xm * ym));
589:   PetscFunctionReturn(PETSC_SUCCESS);
590: }

592: /* ------------------------------------------------------------------- */
593: /*
594:    MSA_BoundaryConditions -  Calculates the boundary conditions for
595:    the region.

597:    Input Parameter:
598: .  user - user-defined application context

600:    Output Parameter:
601: .  user - user-defined application context
602: */
603: static PetscErrorCode MSA_BoundaryConditions(AppCtx *user)
604: {
605:   PetscInt   i, j, k, maxits = 5, limit = 0;
606:   PetscInt   xs, ys, xm, ym, gxs, gys, gxm, gym;
607:   PetscInt   mx = user->mx, my = user->my;
608:   PetscInt   bsize = 0, lsize = 0, tsize = 0, rsize = 0;
609:   PetscReal  one = 1.0, two = 2.0, three = 3.0, scl = 1.0, tol = 1e-10;
610:   PetscReal  fnorm, det, hx, hy, xt = 0, yt = 0;
611:   PetscReal  u1, u2, nf1, nf2, njac11, njac12, njac21, njac22;
612:   PetscReal  b = -0.5, t = 0.5, l = -0.5, r = 0.5;
613:   PetscReal *boundary;
614:   PetscBool  flg;
615:   Vec        Bottom, Top, Right, Left;

617:   PetscFunctionBeginUser;
618:   /* Get local mesh boundaries */
619:   PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
620:   PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL));

622:   bsize = xm + 2;
623:   lsize = ym + 2;
624:   rsize = ym + 2;
625:   tsize = xm + 2;

627:   PetscCall(VecCreateMPI(MPI_COMM_WORLD, bsize, PETSC_DECIDE, &Bottom));
628:   PetscCall(VecCreateMPI(MPI_COMM_WORLD, tsize, PETSC_DECIDE, &Top));
629:   PetscCall(VecCreateMPI(MPI_COMM_WORLD, lsize, PETSC_DECIDE, &Left));
630:   PetscCall(VecCreateMPI(MPI_COMM_WORLD, rsize, PETSC_DECIDE, &Right));

632:   user->Top    = Top;
633:   user->Left   = Left;
634:   user->Bottom = Bottom;
635:   user->Right  = Right;

637:   hx = (r - l) / (mx + 1);
638:   hy = (t - b) / (my + 1);

640:   for (j = 0; j < 4; j++) {
641:     if (j == 0) {
642:       yt    = b;
643:       xt    = l + hx * xs;
644:       limit = bsize;
645:       PetscCall(VecGetArray(Bottom, &boundary));
646:     } else if (j == 1) {
647:       yt    = t;
648:       xt    = l + hx * xs;
649:       limit = tsize;
650:       PetscCall(VecGetArray(Top, &boundary));
651:     } else if (j == 2) {
652:       yt    = b + hy * ys;
653:       xt    = l;
654:       limit = lsize;
655:       PetscCall(VecGetArray(Left, &boundary));
656:     } else if (j == 3) {
657:       yt    = b + hy * ys;
658:       xt    = r;
659:       limit = rsize;
660:       PetscCall(VecGetArray(Right, &boundary));
661:     }

663:     for (i = 0; i < limit; i++) {
664:       u1 = xt;
665:       u2 = -yt;
666:       for (k = 0; k < maxits; k++) {
667:         nf1   = u1 + u1 * u2 * u2 - u1 * u1 * u1 / three - xt;
668:         nf2   = -u2 - u1 * u1 * u2 + u2 * u2 * u2 / three - yt;
669:         fnorm = PetscSqrtScalar(nf1 * nf1 + nf2 * nf2);
670:         if (fnorm <= tol) break;
671:         njac11 = one + u2 * u2 - u1 * u1;
672:         njac12 = two * u1 * u2;
673:         njac21 = -two * u1 * u2;
674:         njac22 = -one - u1 * u1 + u2 * u2;
675:         det    = njac11 * njac22 - njac21 * njac12;
676:         u1     = u1 - (njac22 * nf1 - njac12 * nf2) / det;
677:         u2     = u2 - (njac11 * nf2 - njac21 * nf1) / det;
678:       }

680:       boundary[i] = u1 * u1 - u2 * u2;
681:       if (j == 0 || j == 1) {
682:         xt = xt + hx;
683:       } else if (j == 2 || j == 3) {
684:         yt = yt + hy;
685:       }
686:     }
687:     if (j == 0) {
688:       PetscCall(VecRestoreArray(Bottom, &boundary));
689:     } else if (j == 1) {
690:       PetscCall(VecRestoreArray(Top, &boundary));
691:     } else if (j == 2) {
692:       PetscCall(VecRestoreArray(Left, &boundary));
693:     } else if (j == 3) {
694:       PetscCall(VecRestoreArray(Right, &boundary));
695:     }
696:   }

698:   /* Scale the boundary if desired */

700:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-bottom", &scl, &flg));
701:   if (flg) PetscCall(VecScale(Bottom, scl));
702:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-top", &scl, &flg));
703:   if (flg) PetscCall(VecScale(Top, scl));
704:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-right", &scl, &flg));
705:   if (flg) PetscCall(VecScale(Right, scl));

707:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-left", &scl, &flg));
708:   if (flg) PetscCall(VecScale(Left, scl));
709:   PetscFunctionReturn(PETSC_SUCCESS);
710: }

712: /* ------------------------------------------------------------------- */
713: /*
714:    MSA_Plate -  Calculates an obstacle for surface to stretch over.

716:    Input Parameter:
717: .  user - user-defined application context

719:    Output Parameter:
720: .  user - user-defined application context
721: */
722: static PetscErrorCode MSA_Plate(Vec XL, Vec XU, void *ctx)
723: {
724:   AppCtx    *user = (AppCtx *)ctx;
725:   PetscInt   i, j, row;
726:   PetscInt   xs, ys, xm, ym;
727:   PetscInt   mx = user->mx, my = user->my, bmy, bmx;
728:   PetscReal  t1, t2, t3;
729:   PetscReal *xl, lb = PETSC_NINFINITY, ub = PETSC_INFINITY;
730:   PetscBool  cylinder;

732:   PetscFunctionBeginUser;
733:   user->bmy = PetscMax(0, user->bmy);
734:   user->bmy = PetscMin(my, user->bmy);
735:   user->bmx = PetscMax(0, user->bmx);
736:   user->bmx = PetscMin(mx, user->bmx);
737:   bmy       = user->bmy;
738:   bmx       = user->bmx;

740:   PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));

742:   PetscCall(VecSet(XL, lb));
743:   PetscCall(VecSet(XU, ub));

745:   PetscCall(VecGetArray(XL, &xl));

747:   PetscCall(PetscOptionsHasName(NULL, NULL, "-cylinder", &cylinder));
748:   /* Compute the optional lower box */
749:   if (cylinder) {
750:     for (i = xs; i < xs + xm; i++) {
751:       for (j = ys; j < ys + ym; j++) {
752:         row = (j - ys) * xm + (i - xs);
753:         t1  = (2.0 * i - mx) * bmy;
754:         t2  = (2.0 * j - my) * bmx;
755:         t3  = bmx * bmx * bmy * bmy;
756:         if (t1 * t1 + t2 * t2 <= t3) xl[row] = user->bheight;
757:       }
758:     }
759:   } else {
760:     /* Compute the optional lower box */
761:     for (i = xs; i < xs + xm; i++) {
762:       for (j = ys; j < ys + ym; j++) {
763:         row = (j - ys) * xm + (i - xs);
764:         if (i >= (mx - bmx) / 2 && i < mx - (mx - bmx) / 2 && j >= (my - bmy) / 2 && j < my - (my - bmy) / 2) xl[row] = user->bheight;
765:       }
766:     }
767:   }
768:   PetscCall(VecRestoreArray(XL, &xl));
769:   PetscFunctionReturn(PETSC_SUCCESS);
770: }

772: /* ------------------------------------------------------------------- */
773: /*
774:    MSA_InitialPoint - Calculates the initial guess in one of three ways.

776:    Input Parameters:
777: .  user - user-defined application context
778: .  X - vector for initial guess

780:    Output Parameters:
781: .  X - newly computed initial guess
782: */
783: static PetscErrorCode MSA_InitialPoint(AppCtx *user, Vec X)
784: {
785:   PetscInt  start = -1, i, j;
786:   PetscReal zero  = 0.0;
787:   PetscBool flg;

789:   PetscFunctionBeginUser;
790:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-start", &start, &flg));
791:   if (flg && start == 0) { /* The zero vector is reasonable */
792:     PetscCall(VecSet(X, zero));
793:   } else if (flg && start > 0) { /* Try a random start between -0.5 and 0.5 */
794:     PetscRandom rctx;
795:     PetscReal   np5 = -0.5;

797:     PetscCall(PetscRandomCreate(MPI_COMM_WORLD, &rctx));
798:     for (i = 0; i < start; i++) PetscCall(VecSetRandom(X, rctx));
799:     PetscCall(PetscRandomDestroy(&rctx));
800:     PetscCall(VecShift(X, np5));

802:   } else { /* Take an average of the boundary conditions */

804:     PetscInt   row, xs, xm, gxs, gxm, ys, ym, gys, gym;
805:     PetscInt   mx = user->mx, my = user->my;
806:     PetscReal *x, *left, *right, *bottom, *top;
807:     Vec        localX = user->localX;

809:     /* Get local mesh boundaries */
810:     PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
811:     PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL));

813:     /* Get pointers to vector data */
814:     PetscCall(VecGetArray(user->Top, &top));
815:     PetscCall(VecGetArray(user->Bottom, &bottom));
816:     PetscCall(VecGetArray(user->Left, &left));
817:     PetscCall(VecGetArray(user->Right, &right));

819:     PetscCall(VecGetArray(localX, &x));
820:     /* Perform local computations */
821:     for (j = ys; j < ys + ym; j++) {
822:       for (i = xs; i < xs + xm; i++) {
823:         row    = (j - gys) * gxm + (i - gxs);
824:         x[row] = ((j + 1) * bottom[i - xs + 1] / my + (my - j + 1) * top[i - xs + 1] / (my + 2) + (i + 1) * left[j - ys + 1] / mx + (mx - i + 1) * right[j - ys + 1] / (mx + 2)) / 2.0;
825:       }
826:     }

828:     /* Restore vectors */
829:     PetscCall(VecRestoreArray(localX, &x));

831:     PetscCall(VecRestoreArray(user->Left, &left));
832:     PetscCall(VecRestoreArray(user->Top, &top));
833:     PetscCall(VecRestoreArray(user->Bottom, &bottom));
834:     PetscCall(VecRestoreArray(user->Right, &right));

836:     /* Scatter values into global vector */
837:     PetscCall(DMLocalToGlobalBegin(user->dm, localX, INSERT_VALUES, X));
838:     PetscCall(DMLocalToGlobalEnd(user->dm, localX, INSERT_VALUES, X));
839:   }
840:   PetscFunctionReturn(PETSC_SUCCESS);
841: }

843: /* For testing matrix-free submatrices */
844: PetscErrorCode MatrixFreeHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ptr)
845: {
846:   AppCtx *user = (AppCtx *)ptr;

848:   PetscFunctionBegin;
849:   PetscCall(FormHessian(tao, x, user->H, user->H, ptr));
850:   PetscFunctionReturn(PETSC_SUCCESS);
851: }

853: PetscErrorCode MyMatMult(Mat H_shell, Vec X, Vec Y)
854: {
855:   void   *ptr;
856:   AppCtx *user;

858:   PetscFunctionBegin;
859:   PetscCall(MatShellGetContext(H_shell, &ptr));
860:   user = (AppCtx *)ptr;
861:   PetscCall(MatMult(user->H, X, Y));
862:   PetscFunctionReturn(PETSC_SUCCESS);
863: }

865: /*TEST

867:    build:
868:       requires: !complex

870:    test:
871:       args: -tao_monitor_short -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type tron -tao_gatol 1.e-5
872:       requires: !single

874:    test:
875:       suffix: 2
876:       nsize: 2
877:       args: -tao_monitor_short -mx 8 -my 8 -bmx 2 -bmy 5 -bheight 0.3 -tao_type blmvm -tao_gatol 1.e-4
878:       requires: !single

880:    test:
881:       suffix: 3
882:       nsize: 3
883:       args: -tao_monitor_short -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_type tron -tao_gatol 1.e-5
884:       requires: !single

886:    test:
887:       suffix: 4
888:       nsize: 3
889:       args: -tao_monitor_short -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_subset_type mask -tao_type tron -tao_gatol 1.e-5
890:       requires: !single

892:    test:
893:       suffix: 5
894:       nsize: 3
895:       args: -tao_monitor_short -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_subset_type matrixfree -matrixfree -pc_type none -tao_type tron -tao_gatol 1.e-5
896:       requires: !single

898:    test:
899:       suffix: 6
900:       nsize: 3
901:       args: -tao_monitor_short -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_subset_type matrixfree -matrixfree -tao_type blmvm -tao_gatol 1.e-4
902:       requires: !single

904:    test:
905:       suffix: 7
906:       nsize: 3
907:       args: -tao_monitor_short -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_subset_type matrixfree -pc_type none -tao_type gpcg -tao_gatol 1.e-5
908:       requires: !single

910:    test:
911:       suffix: 8
912:       args: -tao_monitor_short -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bncg -tao_bncg_type gd -tao_gatol 1e-4
913:       requires: !single

915:    test:
916:       suffix: 9
917:       args: -tao_monitor_short -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bncg -tao_gatol 1e-4
918:       requires: !single

920:    test:
921:       suffix: 10
922:       args: -tao_monitor_short -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bnls -tao_gatol 1e-5
923:       requires: !single

925:    test:
926:       suffix: 11
927:       args: -tao_monitor_short -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntr -tao_gatol 1e-5
928:       requires: !single

930:    test:
931:       suffix: 12
932:       args: -tao_monitor_short -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntl -tao_gatol 1e-5
933:       requires: !single

935:    test:
936:       suffix: 13
937:       args: -tao_monitor_short -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bnls -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
938:       requires: !single

940:    test:
941:       suffix: 14
942:       args: -tao_monitor_short -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntr -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
943:       requires: !single

945:    test:
946:       suffix: 15
947:       args: -tao_monitor_short -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntl -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
948:       requires: !single

950:    test:
951:      suffix: 16
952:      args: -tao_monitor_short -mx 8 -my 8 -bmx 2 -bmy 5 -bheight 0.3 -tao_gatol 1e-4 -tao_type bqnls
953:      requires: !single

955:    test:
956:      suffix: 17
957:      args: -tao_monitor_short -mx 8 -my 8 -bmx 2 -bmy 5 -bheight 0.3 -tao_gatol 1e-4 -tao_type bqnkls -tao_bqnk_mat_type lmvmbfgs
958:      requires: !single

960:    test:
961:      suffix: 18
962:      args: -tao_monitor_short -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bnls -tao_gatol 1e-5 -tao_mf_hessian
963:      requires: !single

965:    test:
966:      suffix: 19
967:      args: -tao_monitor_short -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntr -tao_gatol 1e-5 -tao_mf_hessian
968:      requires: !single

970:    test:
971:      suffix: 20
972:      args: -tao_monitor_short -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntl -tao_gatol 1e-5 -tao_mf_hessian
973:      requires: !single

975: TEST*/
```