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 preconditioning matrix

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
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
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
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

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*/