Actual source code: minsurf1.c

      1: /* Program usage: mpiexec -n 1 minsurf1 [-help] [all TAO options] */

      3: /*  Include "petsctao.h" so we can use TAO solvers.  */
      4: #include <petsctao.h>

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

     16: /*
     17:    User-defined application context - contains data needed by the
     18:    application-provided call-back routines, FormFunctionGradient()
     19:    and FormHessian().
     20: */
     21: typedef struct {
     22:   PetscInt   mx, my;                      /* discretization in x, y directions */
     23:   PetscReal *bottom, *top, *left, *right; /* boundary values */
     24:   Mat        H;
     25: } AppCtx;

     27: /* -------- User-defined Routines --------- */

     29: static PetscErrorCode MSA_BoundaryConditions(AppCtx *);
     30: static PetscErrorCode MSA_InitialPoint(AppCtx *, Vec);
     31: static PetscErrorCode QuadraticH(AppCtx *, Vec, Mat);
     32: PetscErrorCode        FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);
     33: PetscErrorCode        FormHessian(Tao, Vec, Mat, Mat, void *);

     35: int main(int argc, char **argv)
     36: {
     37:   PetscInt    N;    /* Size of vector */
     38:   PetscMPIInt size; /* Number of processors */
     39:   Vec         x;    /* solution, gradient vectors */
     40:   KSP         ksp;  /*  PETSc Krylov subspace method */
     41:   PetscBool   flg;  /* A return value when checking for user options */
     42:   Tao         tao;  /* Tao solver context */
     43:   AppCtx      user; /* user-defined work context */

     45:   /* Initialize TAO,PETSc */
     46:   PetscFunctionBeginUser;
     47:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));

     49:   PetscCallMPI(MPI_Comm_size(MPI_COMM_WORLD, &size));
     50:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Incorrect number of processors");

     52:   /* Specify default dimension of the problem */
     53:   user.mx = 4;
     54:   user.my = 4;

     56:   /* Check for any command line arguments that override defaults */
     57:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &user.mx, &flg));
     58:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &user.my, &flg));

     60:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n---- Minimum Surface Area Problem -----\n"));
     61:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "mx: %" PetscInt_FMT "     my: %" PetscInt_FMT "   \n\n", user.mx, user.my));

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

     66:   /* Create TAO solver and set desired solution method  */
     67:   PetscCall(TaoCreate(PETSC_COMM_SELF, &tao));
     68:   PetscCall(TaoSetType(tao, TAOLMVM));

     70:   /* Initialize minsurf application data structure for use in the function evaluations  */
     71:   PetscCall(MSA_BoundaryConditions(&user)); /* Application specific routine */

     73:   /*
     74:      Create a vector to hold the variables.  Compute an initial solution.
     75:      Set this vector, which will be used by TAO.
     76:   */
     77:   PetscCall(VecCreateSeq(PETSC_COMM_SELF, N, &x));
     78:   PetscCall(MSA_InitialPoint(&user, x)); /* Application specific routine */
     79:   PetscCall(TaoSetSolution(tao, x));     /* A TAO routine                */

     81:   /* Provide TAO routines for function, gradient, and Hessian evaluation */
     82:   PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user));

     84:   /* Create a matrix data structure to store the Hessian.  This structure will be used by TAO */
     85:   PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, N, N, 7, NULL, &user.H));
     86:   PetscCall(MatSetOption(user.H, MAT_SYMMETRIC, PETSC_TRUE));
     87:   PetscCall(TaoSetHessian(tao, user.H, user.H, FormHessian, (void *)&user));

     89:   /* Check for any TAO command line options */
     90:   PetscCall(TaoSetFromOptions(tao));

     92:   /* Limit the number of iterations in the KSP linear solver */
     93:   PetscCall(TaoGetKSP(tao, &ksp));
     94:   if (ksp) PetscCall(KSPSetTolerances(ksp, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT, user.mx * user.my));

     96:   /* SOLVE THE APPLICATION */
     97:   PetscCall(TaoSolve(tao));

     99:   PetscCall(TaoDestroy(&tao));
    100:   PetscCall(VecDestroy(&x));
    101:   PetscCall(MatDestroy(&user.H));
    102:   PetscCall(PetscFree(user.bottom));
    103:   PetscCall(PetscFree(user.top));
    104:   PetscCall(PetscFree(user.left));
    105:   PetscCall(PetscFree(user.right));

    107:   PetscCall(PetscFinalize());
    108:   return 0;
    109: }

    111: /* -------------------------------------------------------------------- */

    113: /*  FormFunctionGradient - Evaluates function and corresponding gradient.

    115:     Input Parameters:
    116: .   tao     - the Tao context
    117: .   X       - input vector
    118: .   userCtx - optional user-defined context, as set by TaoSetFunctionGradient()

    120:     Output Parameters:
    121: .   fcn     - the newly evaluated function
    122: .   G       - vector containing the newly evaluated gradient
    123: */
    124: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn, Vec G, void *userCtx)
    125: {
    126:   AppCtx            *user = (AppCtx *)userCtx;
    127:   PetscInt           i, j, row;
    128:   PetscInt           mx = user->mx, my = user->my;
    129:   PetscReal          rhx = mx + 1, rhy = my + 1;
    130:   PetscReal          hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), hydhx = hy / hx, hxdhy = hx / hy, area = 0.5 * hx * hy, ft = 0;
    131:   PetscReal          f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb;
    132:   PetscReal          df1dxc, df2dxc, df3dxc, df4dxc, df5dxc, df6dxc;
    133:   PetscReal          zero = 0.0;
    134:   PetscScalar       *g;
    135:   const PetscScalar *x;

    137:   PetscFunctionBeginUser;
    138:   PetscCall(VecSet(G, zero));

    140:   PetscCall(VecGetArrayRead(X, &x));
    141:   PetscCall(VecGetArray(G, &g));

    143:   /* Compute function over the locally owned part of the mesh */
    144:   for (j = 0; j < my; j++) {
    145:     for (i = 0; i < mx; i++) {
    146:       row = (j)*mx + (i);
    147:       xc  = x[row];
    148:       xlt = xrb = xl = xr = xb = xt = xc;
    149:       if (i == 0) { /* left side */
    150:         xl  = user->left[j + 1];
    151:         xlt = user->left[j + 2];
    152:       } else {
    153:         xl = x[row - 1];
    154:       }

    156:       if (j == 0) { /* bottom side */
    157:         xb  = user->bottom[i + 1];
    158:         xrb = user->bottom[i + 2];
    159:       } else {
    160:         xb = x[row - mx];
    161:       }

    163:       if (i + 1 == mx) { /* right side */
    164:         xr  = user->right[j + 1];
    165:         xrb = user->right[j];
    166:       } else {
    167:         xr = x[row + 1];
    168:       }

    170:       if (j + 1 == 0 + my) { /* top side */
    171:         xt  = user->top[i + 1];
    172:         xlt = user->top[i];
    173:       } else {
    174:         xt = x[row + mx];
    175:       }

    177:       if (i > 0 && j + 1 < my) xlt = x[row - 1 + mx];
    178:       if (j > 0 && i + 1 < mx) xrb = x[row + 1 - mx];

    180:       d1 = (xc - xl);
    181:       d2 = (xc - xr);
    182:       d3 = (xc - xt);
    183:       d4 = (xc - xb);
    184:       d5 = (xr - xrb);
    185:       d6 = (xrb - xb);
    186:       d7 = (xlt - xl);
    187:       d8 = (xt - xlt);

    189:       df1dxc = d1 * hydhx;
    190:       df2dxc = (d1 * hydhx + d4 * hxdhy);
    191:       df3dxc = d3 * hxdhy;
    192:       df4dxc = (d2 * hydhx + d3 * hxdhy);
    193:       df5dxc = d2 * hydhx;
    194:       df6dxc = d4 * hxdhy;

    196:       d1 *= rhx;
    197:       d2 *= rhx;
    198:       d3 *= rhy;
    199:       d4 *= rhy;
    200:       d5 *= rhy;
    201:       d6 *= rhx;
    202:       d7 *= rhy;
    203:       d8 *= rhx;

    205:       f1 = PetscSqrtScalar(1.0 + d1 * d1 + d7 * d7);
    206:       f2 = PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
    207:       f3 = PetscSqrtScalar(1.0 + d3 * d3 + d8 * d8);
    208:       f4 = PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
    209:       f5 = PetscSqrtScalar(1.0 + d2 * d2 + d5 * d5);
    210:       f6 = PetscSqrtScalar(1.0 + d4 * d4 + d6 * d6);

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

    214:       df1dxc /= f1;
    215:       df2dxc /= f2;
    216:       df3dxc /= f3;
    217:       df4dxc /= f4;
    218:       df5dxc /= f5;
    219:       df6dxc /= f6;

    221:       g[row] = (df1dxc + df2dxc + df3dxc + df4dxc + df5dxc + df6dxc) / 2.0;
    222:     }
    223:   }

    225:   for (j = 0; j < my; j++) { /* left side */
    226:     d3 = (user->left[j + 1] - user->left[j + 2]) * rhy;
    227:     d2 = (user->left[j + 1] - x[j * mx]) * rhx;
    228:     ft = ft + PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
    229:   }

    231:   for (i = 0; i < mx; i++) { /* bottom */
    232:     d2 = (user->bottom[i + 1] - user->bottom[i + 2]) * rhx;
    233:     d3 = (user->bottom[i + 1] - x[i]) * rhy;
    234:     ft = ft + PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
    235:   }

    237:   for (j = 0; j < my; j++) { /* right side */
    238:     d1 = (x[(j + 1) * mx - 1] - user->right[j + 1]) * rhx;
    239:     d4 = (user->right[j] - user->right[j + 1]) * rhy;
    240:     ft = ft + PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
    241:   }

    243:   for (i = 0; i < mx; i++) { /* top side */
    244:     d1 = (x[(my - 1) * mx + i] - user->top[i + 1]) * rhy;
    245:     d4 = (user->top[i + 1] - user->top[i]) * rhx;
    246:     ft = ft + PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
    247:   }

    249:   /* Bottom left corner */
    250:   d1 = (user->left[0] - user->left[1]) * rhy;
    251:   d2 = (user->bottom[0] - user->bottom[1]) * rhx;
    252:   ft += PetscSqrtScalar(1.0 + d1 * d1 + d2 * d2);

    254:   /* Top right corner */
    255:   d1 = (user->right[my + 1] - user->right[my]) * rhy;
    256:   d2 = (user->top[mx + 1] - user->top[mx]) * rhx;
    257:   ft += PetscSqrtScalar(1.0 + d1 * d1 + d2 * d2);

    259:   (*fcn) = ft * area;

    261:   /* Restore vectors */
    262:   PetscCall(VecRestoreArrayRead(X, &x));
    263:   PetscCall(VecRestoreArray(G, &g));
    264:   PetscCall(PetscLogFlops(67.0 * mx * my));
    265:   PetscFunctionReturn(PETSC_SUCCESS);
    266: }

    268: /* ------------------------------------------------------------------- */
    269: /*
    270:    FormHessian - Evaluates the Hessian matrix.

    272:    Input Parameters:
    273: .  tao  - the Tao context
    274: .  x    - input vector
    275: .  ptr  - optional user-defined context, as set by TaoSetHessian()

    277:    Output Parameters:
    278: .  H    - Hessian matrix
    279: .  Hpre - optionally different matrix used to compute the preconditioner

    281: */
    282: PetscErrorCode FormHessian(Tao tao, Vec X, Mat H, Mat Hpre, void *ptr)
    283: {
    284:   AppCtx *user = (AppCtx *)ptr;

    286:   PetscFunctionBeginUser;
    287:   /* Evaluate the Hessian entries*/
    288:   PetscCall(QuadraticH(user, X, H));
    289:   PetscFunctionReturn(PETSC_SUCCESS);
    290: }

    292: /* ------------------------------------------------------------------- */
    293: /*
    294:    QuadraticH - Evaluates the Hessian matrix.

    296:    Input Parameters:
    297: .  user - user-defined context, as set by TaoSetHessian()
    298: .  X    - input vector

    300:    Output Parameter:
    301: .  H    - Hessian matrix
    302: */
    303: PetscErrorCode QuadraticH(AppCtx *user, Vec X, Mat Hessian)
    304: {
    305:   PetscInt           i, j, k, row;
    306:   PetscInt           mx = user->mx, my = user->my;
    307:   PetscInt           col[7];
    308:   PetscReal          hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), hydhx = hy / hx, hxdhy = hx / hy;
    309:   PetscReal          rhx = mx + 1, rhy = my + 1;
    310:   PetscReal          f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb;
    311:   PetscReal          hl, hr, ht, hb, hc, htl, hbr;
    312:   const PetscScalar *x;
    313:   PetscReal          v[7];

    315:   PetscFunctionBeginUser;
    316:   /* Get pointers to vector data */
    317:   PetscCall(VecGetArrayRead(X, &x));

    319:   /* Initialize matrix entries to zero */
    320:   PetscCall(MatZeroEntries(Hessian));

    322:   /* Set various matrix options */
    323:   PetscCall(MatSetOption(Hessian, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE));

    325:   /* Compute Hessian over the locally owned part of the mesh */
    326:   for (i = 0; i < mx; i++) {
    327:     for (j = 0; j < my; j++) {
    328:       row = (j)*mx + (i);

    330:       xc  = x[row];
    331:       xlt = xrb = xl = xr = xb = xt = xc;

    333:       /* Left side */
    334:       if (i == 0) {
    335:         xl  = user->left[j + 1];
    336:         xlt = user->left[j + 2];
    337:       } else {
    338:         xl = x[row - 1];
    339:       }

    341:       if (j == 0) {
    342:         xb  = user->bottom[i + 1];
    343:         xrb = user->bottom[i + 2];
    344:       } else {
    345:         xb = x[row - mx];
    346:       }

    348:       if (i + 1 == mx) {
    349:         xr  = user->right[j + 1];
    350:         xrb = user->right[j];
    351:       } else {
    352:         xr = x[row + 1];
    353:       }

    355:       if (j + 1 == my) {
    356:         xt  = user->top[i + 1];
    357:         xlt = user->top[i];
    358:       } else {
    359:         xt = x[row + mx];
    360:       }

    362:       if (i > 0 && j + 1 < my) xlt = x[row - 1 + mx];
    363:       if (j > 0 && i + 1 < mx) xrb = x[row + 1 - mx];

    365:       d1 = (xc - xl) * rhx;
    366:       d2 = (xc - xr) * rhx;
    367:       d3 = (xc - xt) * rhy;
    368:       d4 = (xc - xb) * rhy;
    369:       d5 = (xrb - xr) * rhy;
    370:       d6 = (xrb - xb) * rhx;
    371:       d7 = (xlt - xl) * rhy;
    372:       d8 = (xlt - xt) * rhx;

    374:       f1 = PetscSqrtScalar(1.0 + d1 * d1 + d7 * d7);
    375:       f2 = PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
    376:       f3 = PetscSqrtScalar(1.0 + d3 * d3 + d8 * d8);
    377:       f4 = PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
    378:       f5 = PetscSqrtScalar(1.0 + d2 * d2 + d5 * d5);
    379:       f6 = PetscSqrtScalar(1.0 + d4 * d4 + d6 * d6);

    381:       hl = (-hydhx * (1.0 + d7 * d7) + d1 * d7) / (f1 * f1 * f1) + (-hydhx * (1.0 + d4 * d4) + d1 * d4) / (f2 * f2 * f2);
    382:       hr = (-hydhx * (1.0 + d5 * d5) + d2 * d5) / (f5 * f5 * f5) + (-hydhx * (1.0 + d3 * d3) + d2 * d3) / (f4 * f4 * f4);
    383:       ht = (-hxdhy * (1.0 + d8 * d8) + d3 * d8) / (f3 * f3 * f3) + (-hxdhy * (1.0 + d2 * d2) + d2 * d3) / (f4 * f4 * f4);
    384:       hb = (-hxdhy * (1.0 + d6 * d6) + d4 * d6) / (f6 * f6 * f6) + (-hxdhy * (1.0 + d1 * d1) + d1 * d4) / (f2 * f2 * f2);

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

    389:       hc = hydhx * (1.0 + d7 * d7) / (f1 * f1 * f1) + hxdhy * (1.0 + d8 * d8) / (f3 * f3 * f3) + hydhx * (1.0 + d5 * d5) / (f5 * f5 * f5) + hxdhy * (1.0 + d6 * d6) / (f6 * f6 * f6) + (hxdhy * (1.0 + d1 * d1) + hydhx * (1.0 + d4 * d4) - 2 * d1 * d4) / (f2 * f2 * f2) + (hxdhy * (1.0 + d2 * d2) + hydhx * (1.0 + d3 * d3) - 2 * d2 * d3) / (f4 * f4 * f4);

    391:       hl *= 0.5;
    392:       hr *= 0.5;
    393:       ht *= 0.5;
    394:       hb *= 0.5;
    395:       hbr *= 0.5;
    396:       htl *= 0.5;
    397:       hc *= 0.5;

    399:       k = 0;
    400:       if (j > 0) {
    401:         v[k]   = hb;
    402:         col[k] = row - mx;
    403:         k++;
    404:       }

    406:       if (j > 0 && i < mx - 1) {
    407:         v[k]   = hbr;
    408:         col[k] = row - mx + 1;
    409:         k++;
    410:       }

    412:       if (i > 0) {
    413:         v[k]   = hl;
    414:         col[k] = row - 1;
    415:         k++;
    416:       }

    418:       v[k]   = hc;
    419:       col[k] = row;
    420:       k++;

    422:       if (i < mx - 1) {
    423:         v[k]   = hr;
    424:         col[k] = row + 1;
    425:         k++;
    426:       }

    428:       if (i > 0 && j < my - 1) {
    429:         v[k]   = htl;
    430:         col[k] = row + mx - 1;
    431:         k++;
    432:       }

    434:       if (j < my - 1) {
    435:         v[k]   = ht;
    436:         col[k] = row + mx;
    437:         k++;
    438:       }

    440:       /*
    441:          Set matrix values using local numbering, which was defined
    442:          earlier, in the main routine.
    443:       */
    444:       PetscCall(MatSetValues(Hessian, 1, &row, k, col, v, INSERT_VALUES));
    445:     }
    446:   }

    448:   /* Restore vectors */
    449:   PetscCall(VecRestoreArrayRead(X, &x));

    451:   /* Assemble the matrix */
    452:   PetscCall(MatAssemblyBegin(Hessian, MAT_FINAL_ASSEMBLY));
    453:   PetscCall(MatAssemblyEnd(Hessian, MAT_FINAL_ASSEMBLY));

    455:   PetscCall(PetscLogFlops(199.0 * mx * my));
    456:   PetscFunctionReturn(PETSC_SUCCESS);
    457: }

    459: /* ------------------------------------------------------------------- */
    460: /*
    461:    MSA_BoundaryConditions -  Calculates the boundary conditions for
    462:    the region.

    464:    Input Parameter:
    465: .  user - user-defined application context

    467:    Output Parameter:
    468: .  user - user-defined application context
    469: */
    470: static PetscErrorCode MSA_BoundaryConditions(AppCtx *user)
    471: {
    472:   PetscInt   i, j, k, limit = 0;
    473:   PetscInt   maxits = 5;
    474:   PetscInt   mx = user->mx, my = user->my;
    475:   PetscInt   bsize = 0, lsize = 0, tsize = 0, rsize = 0;
    476:   PetscReal  one = 1.0, two = 2.0, three = 3.0, tol = 1e-10;
    477:   PetscReal  fnorm, det, hx, hy, xt = 0, yt = 0;
    478:   PetscReal  u1, u2, nf1, nf2, njac11, njac12, njac21, njac22;
    479:   PetscReal  b = -0.5, t = 0.5, l = -0.5, r = 0.5;
    480:   PetscReal *boundary;

    482:   PetscFunctionBeginUser;
    483:   bsize = mx + 2;
    484:   lsize = my + 2;
    485:   rsize = my + 2;
    486:   tsize = mx + 2;

    488:   PetscCall(PetscMalloc1(bsize, &user->bottom));
    489:   PetscCall(PetscMalloc1(tsize, &user->top));
    490:   PetscCall(PetscMalloc1(lsize, &user->left));
    491:   PetscCall(PetscMalloc1(rsize, &user->right));

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

    496:   for (j = 0; j < 4; j++) {
    497:     if (j == 0) {
    498:       yt       = b;
    499:       xt       = l;
    500:       limit    = bsize;
    501:       boundary = user->bottom;
    502:     } else if (j == 1) {
    503:       yt       = t;
    504:       xt       = l;
    505:       limit    = tsize;
    506:       boundary = user->top;
    507:     } else if (j == 2) {
    508:       yt       = b;
    509:       xt       = l;
    510:       limit    = lsize;
    511:       boundary = user->left;
    512:     } else { /*  if (j==3) */
    513:       yt       = b;
    514:       xt       = r;
    515:       limit    = rsize;
    516:       boundary = user->right;
    517:     }

    519:     for (i = 0; i < limit; i++) {
    520:       u1 = xt;
    521:       u2 = -yt;
    522:       for (k = 0; k < maxits; k++) {
    523:         nf1   = u1 + u1 * u2 * u2 - u1 * u1 * u1 / three - xt;
    524:         nf2   = -u2 - u1 * u1 * u2 + u2 * u2 * u2 / three - yt;
    525:         fnorm = PetscSqrtScalar(nf1 * nf1 + nf2 * nf2);
    526:         if (fnorm <= tol) break;
    527:         njac11 = one + u2 * u2 - u1 * u1;
    528:         njac12 = two * u1 * u2;
    529:         njac21 = -two * u1 * u2;
    530:         njac22 = -one - u1 * u1 + u2 * u2;
    531:         det    = njac11 * njac22 - njac21 * njac12;
    532:         u1     = u1 - (njac22 * nf1 - njac12 * nf2) / det;
    533:         u2     = u2 - (njac11 * nf2 - njac21 * nf1) / det;
    534:       }

    536:       boundary[i] = u1 * u1 - u2 * u2;
    537:       if (j == 0 || j == 1) {
    538:         xt = xt + hx;
    539:       } else { /*  if (j==2 || j==3) */
    540:         yt = yt + hy;
    541:       }
    542:     }
    543:   }
    544:   PetscFunctionReturn(PETSC_SUCCESS);
    545: }

    547: /* ------------------------------------------------------------------- */
    548: /*
    549:    MSA_InitialPoint - Calculates the initial guess in one of three ways.

    551:    Input Parameters:
    552: .  user - user-defined application context
    553: .  X - vector for initial guess

    555:    Output Parameters:
    556: .  X - newly computed initial guess
    557: */
    558: static PetscErrorCode MSA_InitialPoint(AppCtx *user, Vec X)
    559: {
    560:   PetscInt  start = -1, i, j;
    561:   PetscReal zero  = 0.0;
    562:   PetscBool flg;

    564:   PetscFunctionBeginUser;
    565:   PetscCall(VecSet(X, zero));
    566:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-start", &start, &flg));

    568:   if (flg && start == 0) { /* The zero vector is reasonable */
    569:     PetscCall(VecSet(X, zero));
    570:   } else { /* Take an average of the boundary conditions */
    571:     PetscInt   row;
    572:     PetscInt   mx = user->mx, my = user->my;
    573:     PetscReal *x;

    575:     /* Get pointers to vector data */
    576:     PetscCall(VecGetArray(X, &x));
    577:     /* Perform local computations */
    578:     for (j = 0; j < my; j++) {
    579:       for (i = 0; i < mx; i++) {
    580:         row    = (j)*mx + (i);
    581:         x[row] = (((j + 1) * user->bottom[i + 1] + (my - j + 1) * user->top[i + 1]) / (my + 2) + ((i + 1) * user->left[j + 1] + (mx - i + 1) * user->right[j + 1]) / (mx + 2)) / 2.0;
    582:       }
    583:     }
    584:     /* Restore vectors */
    585:     PetscCall(VecRestoreArray(X, &x));
    586:   }
    587:   PetscFunctionReturn(PETSC_SUCCESS);
    588: }

    590: /*TEST

    592:    build:
    593:       requires: !complex

    595:    test:
    596:       args: -tao_monitor_short -tao_type nls -mx 10 -my 8 -tao_gatol 1.e-4
    597:       requires: !single

    599:    test:
    600:       suffix: 2
    601:       args: -tao_monitor_short -tao_type bmrm -mx 10 -my 8 -tao_gatol 1.e-3
    602:       requires: !single

    604:    test:
    605:       suffix: 3
    606:       args: -tao_monitor_short -tao_type lmvm -mx 10 -my 8 -tao_gatol 1.e-3
    607:       requires: !single

    609:    test:
    610:       suffix: 4
    611:       args: -tao_monitor_short -tao_type bntr -mx 10 -my 8 -tao_gatol 1.e-4
    612:       requires: !single

    614:    test:
    615:       suffix: 4_ew
    616:       args: -tao_monitor_short -tao_type bntr -tao_ksp_ew -mx 10 -my 8 -tao_gatol 1.e-4
    617:       requires: !single

    619:    test:
    620:       suffix: 5
    621:       args: -tao_monitor_short -tao_type bntl -mx 10 -my 8 -tao_gatol 1.e-4
    622:       requires: !single

    624:    test:
    625:       suffix: 5_ew
    626:       args: -tao_monitor_short -tao_type bntl -tao_ksp_ew -mx 10 -my 8 -tao_gatol 1.e-4
    627:       requires: !single

    629:    test:
    630:       suffix: 6
    631:       args: -tao_monitor_short -tao_type bnls -mx 10 -my 8 -tao_gatol 1.e-4
    632:       requires: !single

    634:    test:
    635:       suffix: 6_ew
    636:       args: -tao_monitor_short -tao_type bnls -tao_ksp_ew -tao_bnk_ksp_ew_version 3 -mx 10 -my 8 -tao_gatol 1.e-4
    637:       requires: !single

    639:    test:
    640:       suffix: 7
    641:       args: -tao_monitor_short -tao_type bntr -mx 10 -my 8 -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4 -tao_bnk_cg_tao_monitor_short
    642:       requires: !single

    644:    test:
    645:       suffix: 8
    646:       args: -tao_monitor_short -tao_type bntl -mx 10 -my 8 -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4 -tao_bnk_cg_tao_monitor_short
    647:       requires: !single

    649:    test:
    650:       suffix: 9
    651:       args: -tao_monitor_short -tao_type bnls -mx 10 -my 8 -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4 -tao_bnk_cg_tao_monitor_short
    652:       requires: !single

    654:    test:
    655:       suffix: 10
    656:       args: -tao_monitor_short -tao_type bnls -mx 10 -my 8 -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4 -tao_mf_hessian -tao_bnk_cg_tao_monitor_short

    658:    test:
    659:       suffix: 11
    660:       args: -tao_monitor_short -tao_type bntr -mx 10 -my 8 -tao_gatol 1.e-4 -tao_mf_hessian
    661:       requires: !single

    663:    test:
    664:       suffix: 12
    665:       args: -tao_monitor_short -tao_type bntl -mx 10 -my 8 -tao_gatol 1.e-4 -tao_mf_hessian
    666:       requires: !single

    668: TEST*/