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 */
     61:   PetscFunctionBeginUser;
     62:   PetscCall(PetscInitialize(&argc, &argv, NULL, 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 */
    126:   PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user));

    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, (PetscErrorCodeFn *)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: }

    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:   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:   PetscCallMPI(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 matrix used to construct the preconditioner

    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, PetscCtx 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 -tao_bnk_cg_tao_monitor_short
    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 -tao_bnk_cg_tao_monitor_short
    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 -tao_bnk_cg_tao_monitor_short
    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*/