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
 21:    application-provided call-back routines, FormFunctionGradient(),
 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 */
 62:   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:   PetscOptionsGetInt(NULL, NULL, "-mx", &user.mx, &flg);
 71:   PetscOptionsGetInt(NULL, NULL, "-my", &user.my, &flg);
 72:   PetscOptionsGetReal(NULL, NULL, "-bheight", &user.bheight, &flg);

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

 79:   PetscPrintf(PETSC_COMM_WORLD, "\n---- Minimum Surface Area With Plate Problem -----\n");
 80:   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:   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:   DMSetFromOptions(user.dm);
 96:   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:   DMCreateGlobalVector(user.dm, &x); /* Solution */
104:   DMCreateLocalVector(user.dm, &user.localX);
105:   VecDuplicate(user.localX, &user.localV);

107:   VecDuplicate(x, &xl);
108:   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:   TaoCreate(PETSC_COMM_WORLD, &tao);
118:   TaoSetType(tao, TAOBLMVM);

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

125:   /* Set routines for function, gradient and hessian evaluation */
126:   TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user);

128:   VecGetLocalSize(x, &m);
129:   MatCreateAIJ(MPI_COMM_WORLD, m, m, N, N, 7, NULL, 3, NULL, &(user.H));
130:   MatSetOption(user.H, MAT_SYMMETRIC, PETSC_TRUE);

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

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

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

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

154:   TaoView(tao, PETSC_VIEWER_STDOUT_WORLD);

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

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

176: /*  FormFunctionGradient - Evaluates f(x) and gradient g(x).

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:   /* Get local mesh boundaries */
213:   DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL);
214:   DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

357:   ft = ft * area;
358:   MPI_Allreduce(&ft, fcn, 1, MPIU_REAL, MPIU_SUM, MPI_COMM_WORLD);

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

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

372:   PetscLogFlops(70.0 * xm * ym);

374:   return 0;
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:   /* Set various matrix options */
431:   MatSetOption(Hessian, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

587:   PetscLogFlops(199.0 * xm * ym);
588:   return 0;
589: }

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

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

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

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

620:   bsize = xm + 2;
621:   lsize = ym + 2;
622:   rsize = ym + 2;
623:   tsize = xm + 2;

625:   VecCreateMPI(MPI_COMM_WORLD, bsize, PETSC_DECIDE, &Bottom);
626:   VecCreateMPI(MPI_COMM_WORLD, tsize, PETSC_DECIDE, &Top);
627:   VecCreateMPI(MPI_COMM_WORLD, lsize, PETSC_DECIDE, &Left);
628:   VecCreateMPI(MPI_COMM_WORLD, rsize, PETSC_DECIDE, &Right);

630:   user->Top    = Top;
631:   user->Left   = Left;
632:   user->Bottom = Bottom;
633:   user->Right  = Right;

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

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

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

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

696:   /* Scale the boundary if desired */

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

705:   PetscOptionsGetReal(NULL, NULL, "-left", &scl, &flg);
706:   if (flg) VecScale(Left, scl);
707:   return 0;
708: }

710: /* ------------------------------------------------------------------- */
711: /*
712:    MSA_Plate -  Calculates an obstacle for surface to stretch over.

714:    Input Parameter:
715: .  user - user-defined application context

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

730:   user->bmy = PetscMax(0, user->bmy);
731:   user->bmy = PetscMin(my, user->bmy);
732:   user->bmx = PetscMax(0, user->bmx);
733:   user->bmx = PetscMin(mx, user->bmx);
734:   bmy       = user->bmy;
735:   bmx       = user->bmx;

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

739:   VecSet(XL, lb);
740:   VecSet(XU, ub);

742:   VecGetArray(XL, &xl);

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

767:   return 0;
768: }

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

774:    Input Parameters:
775: .  user - user-defined application context
776: .  X - vector for initial guess

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

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

794:     PetscRandomCreate(MPI_COMM_WORLD, &rctx);
795:     for (i = 0; i < start; i++) VecSetRandom(X, rctx);
796:     PetscRandomDestroy(&rctx);
797:     VecShift(X, np5);

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

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

806:     /* Get local mesh boundaries */
807:     DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL);
808:     DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL);

810:     /* Get pointers to vector data */
811:     VecGetArray(user->Top, &top);
812:     VecGetArray(user->Bottom, &bottom);
813:     VecGetArray(user->Left, &left);
814:     VecGetArray(user->Right, &right);

816:     VecGetArray(localX, &x);
817:     /* Perform local computations */
818:     for (j = ys; j < ys + ym; j++) {
819:       for (i = xs; i < xs + xm; i++) {
820:         row    = (j - gys) * gxm + (i - gxs);
821:         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;
822:       }
823:     }

825:     /* Restore vectors */
826:     VecRestoreArray(localX, &x);

828:     VecRestoreArray(user->Left, &left);
829:     VecRestoreArray(user->Top, &top);
830:     VecRestoreArray(user->Bottom, &bottom);
831:     VecRestoreArray(user->Right, &right);

833:     /* Scatter values into global vector */
834:     DMLocalToGlobalBegin(user->dm, localX, INSERT_VALUES, X);
835:     DMLocalToGlobalEnd(user->dm, localX, INSERT_VALUES, X);
836:   }
837:   return 0;
838: }

840: /* For testing matrix free submatrices */
841: PetscErrorCode MatrixFreeHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ptr)
842: {
843:   AppCtx *user = (AppCtx *)ptr;
844:   FormHessian(tao, x, user->H, user->H, ptr);
845:   return 0;
846: }
847: PetscErrorCode MyMatMult(Mat H_shell, Vec X, Vec Y)
848: {
849:   void   *ptr;
850:   AppCtx *user;
851:   MatShellGetContext(H_shell, &ptr);
852:   user = (AppCtx *)ptr;
853:   MatMult(user->H, X, Y);
854:   return 0;
855: }

857: /*TEST

859:    build:
860:       requires: !complex

862:    test:
863:       args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type tron -tao_gatol 1.e-5
864:       requires: !single

866:    test:
867:       suffix: 2
868:       nsize: 2
869:       args: -tao_smonitor -mx 8 -my 8 -bmx 2 -bmy 5 -bheight 0.3 -tao_type blmvm -tao_gatol 1.e-4
870:       requires: !single

872:    test:
873:       suffix: 3
874:       nsize: 3
875:       args: -tao_smonitor -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_type tron -tao_gatol 1.e-5
876:       requires: !single

878:    test:
879:       suffix: 4
880:       nsize: 3
881:       args: -tao_smonitor -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_subset_type mask -tao_type tron -tao_gatol 1.e-5
882:       requires: !single

884:    test:
885:       suffix: 5
886:       nsize: 3
887:       args: -tao_smonitor -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
888:       requires: !single

890:    test:
891:       suffix: 6
892:       nsize: 3
893:       args: -tao_smonitor -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_subset_type matrixfree -matrixfree -tao_type blmvm -tao_gatol 1.e-4
894:       requires: !single

896:    test:
897:       suffix: 7
898:       nsize: 3
899:       args: -tao_smonitor -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
900:       requires: !single

902:    test:
903:       suffix: 8
904:       args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bncg -tao_bncg_type gd -tao_gatol 1e-4
905:       requires: !single

907:    test:
908:       suffix: 9
909:       args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bncg -tao_gatol 1e-4
910:       requires: !single

912:    test:
913:       suffix: 10
914:       args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bnls -tao_gatol 1e-5
915:       requires: !single

917:    test:
918:       suffix: 11
919:       args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntr -tao_gatol 1e-5
920:       requires: !single

922:    test:
923:       suffix: 12
924:       args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntl -tao_gatol 1e-5
925:       requires: !single

927:    test:
928:       suffix: 13
929:       args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bnls -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
930:       requires: !single

932:    test:
933:       suffix: 14
934:       args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntr -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
935:       requires: !single

937:    test:
938:       suffix: 15
939:       args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntl -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
940:       requires: !single

942:    test:
943:      suffix: 16
944:      args: -tao_smonitor -mx 8 -my 8 -bmx 2 -bmy 5 -bheight 0.3 -tao_gatol 1e-4 -tao_type bqnls
945:      requires: !single

947:    test:
948:      suffix: 17
949:      args: -tao_smonitor -mx 8 -my 8 -bmx 2 -bmy 5 -bheight 0.3 -tao_gatol 1e-4 -tao_type bqnkls -tao_bqnk_mat_type lmvmbfgs
950:      requires: !single

952:    test:
953:      suffix: 18
954:      args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bnls -tao_gatol 1e-5 -tao_mf_hessian
955:      requires: !single

957:    test:
958:      suffix: 19
959:      args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntr -tao_gatol 1e-5 -tao_mf_hessian
960:      requires: !single

962:    test:
963:      suffix: 20
964:      args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntl -tao_gatol 1e-5 -tao_mf_hessian

966: TEST*/