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, (char *)0, 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_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, 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
280: .  flg  - flag indicating matrix structure

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

576:     /* Get pointers to vector data */
577:     PetscCall(VecGetArray(X, &x));
578:     /* Perform local computations */
579:     for (j = 0; j < my; j++) {
580:       for (i = 0; i < mx; i++) {
581:         row    = (j)*mx + (i);
582:         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;
583:       }
584:     }
585:     /* Restore vectors */
586:     PetscCall(VecRestoreArray(X, &x));
587:   }
588:   PetscFunctionReturn(PETSC_SUCCESS);
589: }

591: /*TEST

593:    build:
594:       requires: !complex

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

669: TEST*/