Actual source code: minsurf2.c

  1: /* Program usage: mpiexec -n <proc> minsurf2 [-help] [all TAO options] */

  3: /*
  4:   Include "petsctao.h" so we can use TAO solvers.
  5:   petscdmda.h for distributed array
  6: */
  7: #include <petsctao.h>
  8: #include <petscdmda.h>

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

 21: /*
 22:    User-defined application context - contains data needed by the
 23:    application-provided call-back routines, FormFunctionGradient()
 24:    and FormHessian().
 25: */
 26: typedef struct {
 27:   PetscInt   mx, my;                      /* discretization in x, y directions */
 28:   PetscReal *bottom, *top, *left, *right; /* boundary values */
 29:   DM         dm;                          /* distributed array data structure */
 30:   Mat        H;                           /* Hessian */
 31: } AppCtx;

 33: /* -------- User-defined Routines --------- */

 35: static PetscErrorCode MSA_BoundaryConditions(AppCtx *);
 36: static PetscErrorCode MSA_InitialPoint(AppCtx *, Vec);
 37: PetscErrorCode        QuadraticH(AppCtx *, Vec, Mat);
 38: PetscErrorCode        FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);
 39: PetscErrorCode        FormGradient(Tao, Vec, Vec, void *);
 40: PetscErrorCode        FormHessian(Tao, Vec, Mat, Mat, void *);
 41: PetscErrorCode        My_Monitor(Tao, void *);

 43: int main(int argc, char **argv)
 44: {
 45:   PetscInt      Nx, Ny;                /* number of processors in x- and y- directions */
 46:   Vec           x;                     /* solution, gradient vectors */
 47:   PetscBool     flg, viewmat;          /* flags */
 48:   PetscBool     fddefault, fdcoloring; /* flags */
 49:   Tao           tao;                   /* TAO solver context */
 50:   AppCtx        user;                  /* user-defined work context */
 51:   ISColoring    iscoloring;
 52:   MatFDColoring matfdcoloring;

 54:   /* Initialize TAO */
 56:   PetscInitialize(&argc, &argv, (char *)0, help);

 58:   /* Specify dimension of the problem */
 59:   user.mx = 10;
 60:   user.my = 10;

 62:   /* Check for any command line arguments that override defaults */
 63:   PetscOptionsGetInt(NULL, NULL, "-mx", &user.mx, &flg);
 64:   PetscOptionsGetInt(NULL, NULL, "-my", &user.my, &flg);

 66:   PetscPrintf(MPI_COMM_WORLD, "\n---- Minimum Surface Area Problem -----\n");
 67:   PetscPrintf(MPI_COMM_WORLD, "mx: %" PetscInt_FMT "     my: %" PetscInt_FMT "   \n\n", user.mx, user.my);

 69:   /* Let PETSc determine the vector distribution */
 70:   Nx = PETSC_DECIDE;
 71:   Ny = PETSC_DECIDE;

 73:   /* Create distributed array (DM) to manage parallel grid and vectors  */
 74:   DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, user.mx, user.my, Nx, Ny, 1, 1, NULL, NULL, &user.dm);
 75:   DMSetFromOptions(user.dm);
 76:   DMSetUp(user.dm);

 78:   /* Create TAO solver and set desired solution method.*/
 79:   TaoCreate(PETSC_COMM_WORLD, &tao);
 80:   TaoSetType(tao, TAOCG);

 82:   /*
 83:      Extract global vector from DA for the vector of variables --  PETSC routine
 84:      Compute the initial solution                              --  application specific, see below
 85:      Set this vector for use by TAO                            --  TAO routine
 86:   */
 87:   DMCreateGlobalVector(user.dm, &x);
 88:   MSA_BoundaryConditions(&user);
 89:   MSA_InitialPoint(&user, x);
 90:   TaoSetSolution(tao, x);

 92:   /*
 93:      Initialize the Application context for use in function evaluations  --  application specific, see below.
 94:      Set routines for function and gradient evaluation
 95:   */
 96:   TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user);

 98:   /*
 99:      Given the command line arguments, calculate the hessian with either the user-
100:      provided function FormHessian, or the default finite-difference driven Hessian
101:      functions
102:   */
103:   PetscOptionsHasName(NULL, NULL, "-tao_fddefault", &fddefault);
104:   PetscOptionsHasName(NULL, NULL, "-tao_fdcoloring", &fdcoloring);

106:   /*
107:      Create a matrix data structure to store the Hessian and set
108:      the Hessian evaluation routine.
109:      Set the matrix structure to be used for Hessian evaluations
110:   */
111:   DMCreateMatrix(user.dm, &user.H);
112:   MatSetOption(user.H, MAT_SYMMETRIC, PETSC_TRUE);

114:   if (fdcoloring) {
115:     DMCreateColoring(user.dm, IS_COLORING_GLOBAL, &iscoloring);
116:     MatFDColoringCreate(user.H, iscoloring, &matfdcoloring);
117:     ISColoringDestroy(&iscoloring);
118:     MatFDColoringSetFunction(matfdcoloring, (PetscErrorCode(*)(void))FormGradient, (void *)&user);
119:     MatFDColoringSetFromOptions(matfdcoloring);
120:     TaoSetHessian(tao, user.H, user.H, TaoDefaultComputeHessianColor, (void *)matfdcoloring);
121:   } else if (fddefault) {
122:     TaoSetHessian(tao, user.H, user.H, TaoDefaultComputeHessian, (void *)NULL);
123:   } else {
124:     TaoSetHessian(tao, user.H, user.H, FormHessian, (void *)&user);
125:   }

127:   /*
128:      If my_monitor option is in command line, then use the user-provided
129:      monitoring function
130:   */
131:   PetscOptionsHasName(NULL, NULL, "-my_monitor", &viewmat);
132:   if (viewmat) TaoSetMonitor(tao, My_Monitor, NULL, NULL);

134:   /* Check for any tao command line options */
135:   TaoSetFromOptions(tao);

137:   /* SOLVE THE APPLICATION */
138:   TaoSolve(tao);

140:   TaoView(tao, PETSC_VIEWER_STDOUT_WORLD);

142:   /* Free TAO data structures */
143:   TaoDestroy(&tao);

145:   /* Free PETSc data structures */
146:   VecDestroy(&x);
147:   MatDestroy(&user.H);
148:   if (fdcoloring) MatFDColoringDestroy(&matfdcoloring);
149:   PetscFree(user.bottom);
150:   PetscFree(user.top);
151:   PetscFree(user.left);
152:   PetscFree(user.right);
153:   DMDestroy(&user.dm);
154:   PetscFinalize();
155:   return 0;
156: }

158: PetscErrorCode FormGradient(Tao tao, Vec X, Vec G, void *userCtx)
159: {
160:   PetscReal fcn;

162:   FormFunctionGradient(tao, X, &fcn, G, userCtx);
163:   return 0;
164: }

166: /* -------------------------------------------------------------------- */
167: /*  FormFunctionGradient - Evaluates the function and corresponding gradient.

169:     Input Parameters:
170: .   tao     - the Tao context
171: .   XX      - input vector
172: .   userCtx - optional user-defined context, as set by TaoSetObjectiveAndGradient()

174:     Output Parameters:
175: .   fcn     - the newly evaluated function
176: .   GG       - vector containing the newly evaluated gradient
177: */
178: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn, Vec G, void *userCtx)
179: {
180:   AppCtx     *user = (AppCtx *)userCtx;
181:   PetscInt    i, j;
182:   PetscInt    mx = user->mx, my = user->my;
183:   PetscInt    xs, xm, gxs, gxm, ys, ym, gys, gym;
184:   PetscReal   ft = 0.0;
185:   PetscReal   hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), hydhx = hy / hx, hxdhy = hx / hy, area = 0.5 * hx * hy;
186:   PetscReal   rhx = mx + 1, rhy = my + 1;
187:   PetscReal   f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb;
188:   PetscReal   df1dxc, df2dxc, df3dxc, df4dxc, df5dxc, df6dxc;
189:   PetscReal **g, **x;
190:   Vec         localX;

192:   /* Get local mesh boundaries */
193:   DMGetLocalVector(user->dm, &localX);
194:   DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL);
195:   DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL);

197:   /* Scatter ghost points to local vector */
198:   DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, localX);
199:   DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, localX);

201:   /* Get pointers to vector data */
202:   DMDAVecGetArray(user->dm, localX, (void **)&x);
203:   DMDAVecGetArray(user->dm, G, (void **)&g);

205:   /* Compute function and gradient over the locally owned part of the mesh */
206:   for (j = ys; j < ys + ym; j++) {
207:     for (i = xs; i < xs + xm; i++) {
208:       xc  = x[j][i];
209:       xlt = xrb = xl = xr = xb = xt = xc;

211:       if (i == 0) { /* left side */
212:         xl  = user->left[j - ys + 1];
213:         xlt = user->left[j - ys + 2];
214:       } else {
215:         xl = x[j][i - 1];
216:       }

218:       if (j == 0) { /* bottom side */
219:         xb  = user->bottom[i - xs + 1];
220:         xrb = user->bottom[i - xs + 2];
221:       } else {
222:         xb = x[j - 1][i];
223:       }

225:       if (i + 1 == gxs + gxm) { /* right side */
226:         xr  = user->right[j - ys + 1];
227:         xrb = user->right[j - ys];
228:       } else {
229:         xr = x[j][i + 1];
230:       }

232:       if (j + 1 == gys + gym) { /* top side */
233:         xt  = user->top[i - xs + 1];
234:         xlt = user->top[i - xs];
235:       } else {
236:         xt = x[j + 1][i];
237:       }

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

242:       d1 = (xc - xl);
243:       d2 = (xc - xr);
244:       d3 = (xc - xt);
245:       d4 = (xc - xb);
246:       d5 = (xr - xrb);
247:       d6 = (xrb - xb);
248:       d7 = (xlt - xl);
249:       d8 = (xt - xlt);

251:       df1dxc = d1 * hydhx;
252:       df2dxc = (d1 * hydhx + d4 * hxdhy);
253:       df3dxc = d3 * hxdhy;
254:       df4dxc = (d2 * hydhx + d3 * hxdhy);
255:       df5dxc = d2 * hydhx;
256:       df6dxc = d4 * hxdhy;

258:       d1 *= rhx;
259:       d2 *= rhx;
260:       d3 *= rhy;
261:       d4 *= rhy;
262:       d5 *= rhy;
263:       d6 *= rhx;
264:       d7 *= rhy;
265:       d8 *= rhx;

267:       f1 = PetscSqrtReal(1.0 + d1 * d1 + d7 * d7);
268:       f2 = PetscSqrtReal(1.0 + d1 * d1 + d4 * d4);
269:       f3 = PetscSqrtReal(1.0 + d3 * d3 + d8 * d8);
270:       f4 = PetscSqrtReal(1.0 + d3 * d3 + d2 * d2);
271:       f5 = PetscSqrtReal(1.0 + d2 * d2 + d5 * d5);
272:       f6 = PetscSqrtReal(1.0 + d4 * d4 + d6 * d6);

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

276:       df1dxc /= f1;
277:       df2dxc /= f2;
278:       df3dxc /= f3;
279:       df4dxc /= f4;
280:       df5dxc /= f5;
281:       df6dxc /= f6;

283:       g[j][i] = (df1dxc + df2dxc + df3dxc + df4dxc + df5dxc + df6dxc) * 0.5;
284:     }
285:   }

287:   /* Compute triangular areas along the border of the domain. */
288:   if (xs == 0) { /* left side */
289:     for (j = ys; j < ys + ym; j++) {
290:       d3 = (user->left[j - ys + 1] - user->left[j - ys + 2]) * rhy;
291:       d2 = (user->left[j - ys + 1] - x[j][0]) * rhx;
292:       ft = ft + PetscSqrtReal(1.0 + d3 * d3 + d2 * d2);
293:     }
294:   }
295:   if (ys == 0) { /* bottom side */
296:     for (i = xs; i < xs + xm; i++) {
297:       d2 = (user->bottom[i + 1 - xs] - user->bottom[i - xs + 2]) * rhx;
298:       d3 = (user->bottom[i - xs + 1] - x[0][i]) * rhy;
299:       ft = ft + PetscSqrtReal(1.0 + d3 * d3 + d2 * d2);
300:     }
301:   }

303:   if (xs + xm == mx) { /* right side */
304:     for (j = ys; j < ys + ym; j++) {
305:       d1 = (x[j][mx - 1] - user->right[j - ys + 1]) * rhx;
306:       d4 = (user->right[j - ys] - user->right[j - ys + 1]) * rhy;
307:       ft = ft + PetscSqrtReal(1.0 + d1 * d1 + d4 * d4);
308:     }
309:   }
310:   if (ys + ym == my) { /* top side */
311:     for (i = xs; i < xs + xm; i++) {
312:       d1 = (x[my - 1][i] - user->top[i - xs + 1]) * rhy;
313:       d4 = (user->top[i - xs + 1] - user->top[i - xs]) * rhx;
314:       ft = ft + PetscSqrtReal(1.0 + d1 * d1 + d4 * d4);
315:     }
316:   }

318:   if (ys == 0 && xs == 0) {
319:     d1 = (user->left[0] - user->left[1]) / hy;
320:     d2 = (user->bottom[0] - user->bottom[1]) * rhx;
321:     ft += PetscSqrtReal(1.0 + d1 * d1 + d2 * d2);
322:   }
323:   if (ys + ym == my && xs + xm == mx) {
324:     d1 = (user->right[ym + 1] - user->right[ym]) * rhy;
325:     d2 = (user->top[xm + 1] - user->top[xm]) * rhx;
326:     ft += PetscSqrtReal(1.0 + d1 * d1 + d2 * d2);
327:   }

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

332:   /* Restore vectors */
333:   DMDAVecRestoreArray(user->dm, localX, (void **)&x);
334:   DMDAVecRestoreArray(user->dm, G, (void **)&g);

336:   /* Scatter values to global vector */
337:   DMRestoreLocalVector(user->dm, &localX);
338:   PetscLogFlops(67.0 * xm * ym);
339:   return 0;
340: }

342: /* ------------------------------------------------------------------- */
343: /*
344:    FormHessian - Evaluates Hessian matrix.

346:    Input Parameters:
347: .  tao  - the Tao context
348: .  x    - input vector
349: .  ptr  - optional user-defined context, as set by TaoSetHessian()

351:    Output Parameters:
352: .  H    - Hessian matrix
353: .  Hpre - optionally different preconditioning matrix
354: .  flg  - flag indicating matrix structure

356: */
357: PetscErrorCode FormHessian(Tao tao, Vec X, Mat H, Mat Hpre, void *ptr)
358: {
359:   AppCtx *user = (AppCtx *)ptr;

361:   /* Evaluate the Hessian entries*/
362:   QuadraticH(user, X, H);
363:   return 0;
364: }

366: /* ------------------------------------------------------------------- */
367: /*
368:    QuadraticH - Evaluates Hessian matrix.

370:    Input Parameters:
371: .  user - user-defined context, as set by TaoSetHessian()
372: .  X    - input vector

374:    Output Parameter:
375: .  H    - Hessian matrix
376: */
377: PetscErrorCode QuadraticH(AppCtx *user, Vec X, Mat Hessian)
378: {
379:   PetscInt    i, j, k;
380:   PetscInt    mx = user->mx, my = user->my;
381:   PetscInt    xs, xm, gxs, gxm, ys, ym, gys, gym;
382:   PetscReal   hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), hydhx = hy / hx, hxdhy = hx / hy;
383:   PetscReal   f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb;
384:   PetscReal   hl, hr, ht, hb, hc, htl, hbr;
385:   PetscReal **x, v[7];
386:   MatStencil  col[7], row;
387:   Vec         localX;
388:   PetscBool   assembled;

390:   /* Get local mesh boundaries */
391:   DMGetLocalVector(user->dm, &localX);

393:   DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL);
394:   DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL);

396:   /* Scatter ghost points to local vector */
397:   DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, localX);
398:   DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, localX);

400:   /* Get pointers to vector data */
401:   DMDAVecGetArray(user->dm, localX, (void **)&x);

403:   /* Initialize matrix entries to zero */
404:   MatAssembled(Hessian, &assembled);
405:   if (assembled) MatZeroEntries(Hessian);

407:   /* Set various matrix options */
408:   MatSetOption(Hessian, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE);

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

412:   for (j = ys; j < ys + ym; j++) {
413:     for (i = xs; i < xs + xm; i++) {
414:       xc  = x[j][i];
415:       xlt = xrb = xl = xr = xb = xt = xc;

417:       /* Left side */
418:       if (i == 0) {
419:         xl  = user->left[j - ys + 1];
420:         xlt = user->left[j - ys + 2];
421:       } else {
422:         xl = x[j][i - 1];
423:       }

425:       if (j == 0) {
426:         xb  = user->bottom[i - xs + 1];
427:         xrb = user->bottom[i - xs + 2];
428:       } else {
429:         xb = x[j - 1][i];
430:       }

432:       if (i + 1 == mx) {
433:         xr  = user->right[j - ys + 1];
434:         xrb = user->right[j - ys];
435:       } else {
436:         xr = x[j][i + 1];
437:       }

439:       if (j + 1 == my) {
440:         xt  = user->top[i - xs + 1];
441:         xlt = user->top[i - xs];
442:       } else {
443:         xt = x[j + 1][i];
444:       }

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

449:       d1 = (xc - xl) / hx;
450:       d2 = (xc - xr) / hx;
451:       d3 = (xc - xt) / hy;
452:       d4 = (xc - xb) / hy;
453:       d5 = (xrb - xr) / hy;
454:       d6 = (xrb - xb) / hx;
455:       d7 = (xlt - xl) / hy;
456:       d8 = (xlt - xt) / hx;

458:       f1 = PetscSqrtReal(1.0 + d1 * d1 + d7 * d7);
459:       f2 = PetscSqrtReal(1.0 + d1 * d1 + d4 * d4);
460:       f3 = PetscSqrtReal(1.0 + d3 * d3 + d8 * d8);
461:       f4 = PetscSqrtReal(1.0 + d3 * d3 + d2 * d2);
462:       f5 = PetscSqrtReal(1.0 + d2 * d2 + d5 * d5);
463:       f6 = PetscSqrtReal(1.0 + d4 * d4 + d6 * d6);

465:       hl = (-hydhx * (1.0 + d7 * d7) + d1 * d7) / (f1 * f1 * f1) + (-hydhx * (1.0 + d4 * d4) + d1 * d4) / (f2 * f2 * f2);
466:       hr = (-hydhx * (1.0 + d5 * d5) + d2 * d5) / (f5 * f5 * f5) + (-hydhx * (1.0 + d3 * d3) + d2 * d3) / (f4 * f4 * f4);
467:       ht = (-hxdhy * (1.0 + d8 * d8) + d3 * d8) / (f3 * f3 * f3) + (-hxdhy * (1.0 + d2 * d2) + d2 * d3) / (f4 * f4 * f4);
468:       hb = (-hxdhy * (1.0 + d6 * d6) + d4 * d6) / (f6 * f6 * f6) + (-hxdhy * (1.0 + d1 * d1) + d1 * d4) / (f2 * f2 * f2);

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

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

475:       hl /= 2.0;
476:       hr /= 2.0;
477:       ht /= 2.0;
478:       hb /= 2.0;
479:       hbr /= 2.0;
480:       htl /= 2.0;
481:       hc /= 2.0;

483:       row.j = j;
484:       row.i = i;
485:       k     = 0;
486:       if (j > 0) {
487:         v[k]     = hb;
488:         col[k].j = j - 1;
489:         col[k].i = i;
490:         k++;
491:       }

493:       if (j > 0 && i < mx - 1) {
494:         v[k]     = hbr;
495:         col[k].j = j - 1;
496:         col[k].i = i + 1;
497:         k++;
498:       }

500:       if (i > 0) {
501:         v[k]     = hl;
502:         col[k].j = j;
503:         col[k].i = i - 1;
504:         k++;
505:       }

507:       v[k]     = hc;
508:       col[k].j = j;
509:       col[k].i = i;
510:       k++;

512:       if (i < mx - 1) {
513:         v[k]     = hr;
514:         col[k].j = j;
515:         col[k].i = i + 1;
516:         k++;
517:       }

519:       if (i > 0 && j < my - 1) {
520:         v[k]     = htl;
521:         col[k].j = j + 1;
522:         col[k].i = i - 1;
523:         k++;
524:       }

526:       if (j < my - 1) {
527:         v[k]     = ht;
528:         col[k].j = j + 1;
529:         col[k].i = i;
530:         k++;
531:       }

533:       /*
534:          Set matrix values using local numbering, which was defined
535:          earlier, in the main routine.
536:       */
537:       MatSetValuesStencil(Hessian, 1, &row, k, col, v, INSERT_VALUES);
538:     }
539:   }

541:   DMDAVecRestoreArray(user->dm, localX, (void **)&x);
542:   DMRestoreLocalVector(user->dm, &localX);

544:   MatAssemblyBegin(Hessian, MAT_FINAL_ASSEMBLY);
545:   MatAssemblyEnd(Hessian, MAT_FINAL_ASSEMBLY);

547:   PetscLogFlops(199.0 * xm * ym);
548:   return 0;
549: }

551: /* ------------------------------------------------------------------- */
552: /*
553:    MSA_BoundaryConditions -  Calculates the boundary conditions for
554:    the region.

556:    Input Parameter:
557: .  user - user-defined application context

559:    Output Parameter:
560: .  user - user-defined application context
561: */
562: static PetscErrorCode MSA_BoundaryConditions(AppCtx *user)
563: {
564:   PetscInt   i, j, k, limit = 0, maxits = 5;
565:   PetscInt   xs, ys, xm, ym, gxs, gys, gxm, gym;
566:   PetscInt   mx = user->mx, my = user->my;
567:   PetscInt   bsize = 0, lsize = 0, tsize = 0, rsize = 0;
568:   PetscReal  one = 1.0, two = 2.0, three = 3.0, tol = 1e-10;
569:   PetscReal  fnorm, det, hx, hy, xt = 0, yt = 0;
570:   PetscReal  u1, u2, nf1, nf2, njac11, njac12, njac21, njac22;
571:   PetscReal  b = -0.5, t = 0.5, l = -0.5, r = 0.5;
572:   PetscReal *boundary;
573:   PetscBool  flg;

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

579:   bsize = xm + 2;
580:   lsize = ym + 2;
581:   rsize = ym + 2;
582:   tsize = xm + 2;

584:   PetscMalloc1(bsize, &user->bottom);
585:   PetscMalloc1(tsize, &user->top);
586:   PetscMalloc1(lsize, &user->left);
587:   PetscMalloc1(rsize, &user->right);

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

592:   for (j = 0; j < 4; j++) {
593:     if (j == 0) {
594:       yt       = b;
595:       xt       = l + hx * xs;
596:       limit    = bsize;
597:       boundary = user->bottom;
598:     } else if (j == 1) {
599:       yt       = t;
600:       xt       = l + hx * xs;
601:       limit    = tsize;
602:       boundary = user->top;
603:     } else if (j == 2) {
604:       yt       = b + hy * ys;
605:       xt       = l;
606:       limit    = lsize;
607:       boundary = user->left;
608:     } else { /* if (j==3) */
609:       yt       = b + hy * ys;
610:       xt       = r;
611:       limit    = rsize;
612:       boundary = user->right;
613:     }

615:     for (i = 0; i < limit; i++) {
616:       u1 = xt;
617:       u2 = -yt;
618:       for (k = 0; k < maxits; k++) {
619:         nf1   = u1 + u1 * u2 * u2 - u1 * u1 * u1 / three - xt;
620:         nf2   = -u2 - u1 * u1 * u2 + u2 * u2 * u2 / three - yt;
621:         fnorm = PetscSqrtReal(nf1 * nf1 + nf2 * nf2);
622:         if (fnorm <= tol) break;
623:         njac11 = one + u2 * u2 - u1 * u1;
624:         njac12 = two * u1 * u2;
625:         njac21 = -two * u1 * u2;
626:         njac22 = -one - u1 * u1 + u2 * u2;
627:         det    = njac11 * njac22 - njac21 * njac12;
628:         u1     = u1 - (njac22 * nf1 - njac12 * nf2) / det;
629:         u2     = u2 - (njac11 * nf2 - njac21 * nf1) / det;
630:       }

632:       boundary[i] = u1 * u1 - u2 * u2;
633:       if (j == 0 || j == 1) {
634:         xt = xt + hx;
635:       } else { /*  if (j==2 || j==3) */
636:         yt = yt + hy;
637:       }
638:     }
639:   }

641:   /* Scale the boundary if desired */
642:   if (1 == 1) {
643:     PetscReal scl = 1.0;

645:     PetscOptionsGetReal(NULL, NULL, "-bottom", &scl, &flg);
646:     if (flg) {
647:       for (i = 0; i < bsize; i++) user->bottom[i] *= scl;
648:     }

650:     PetscOptionsGetReal(NULL, NULL, "-top", &scl, &flg);
651:     if (flg) {
652:       for (i = 0; i < tsize; i++) user->top[i] *= scl;
653:     }

655:     PetscOptionsGetReal(NULL, NULL, "-right", &scl, &flg);
656:     if (flg) {
657:       for (i = 0; i < rsize; i++) user->right[i] *= scl;
658:     }

660:     PetscOptionsGetReal(NULL, NULL, "-left", &scl, &flg);
661:     if (flg) {
662:       for (i = 0; i < lsize; i++) user->left[i] *= scl;
663:     }
664:   }
665:   return 0;
666: }

668: /* ------------------------------------------------------------------- */
669: /*
670:    MSA_InitialPoint - Calculates the initial guess in one of three ways.

672:    Input Parameters:
673: .  user - user-defined application context
674: .  X - vector for initial guess

676:    Output Parameters:
677: .  X - newly computed initial guess
678: */
679: static PetscErrorCode MSA_InitialPoint(AppCtx *user, Vec X)
680: {
681:   PetscInt  start2 = -1, i, j;
682:   PetscReal start1 = 0;
683:   PetscBool flg1, flg2;

685:   PetscOptionsGetReal(NULL, NULL, "-start", &start1, &flg1);
686:   PetscOptionsGetInt(NULL, NULL, "-random", &start2, &flg2);

688:   if (flg1) { /* The zero vector is reasonable */

690:     VecSet(X, start1);

692:   } else if (flg2 && start2 > 0) { /* Try a random start between -0.5 and 0.5 */
693:     PetscRandom rctx;
694:     PetscReal   np5 = -0.5;

696:     PetscRandomCreate(PETSC_COMM_WORLD, &rctx);
697:     VecSetRandom(X, rctx);
698:     PetscRandomDestroy(&rctx);
699:     VecShift(X, np5);

701:   } else { /* Take an average of the boundary conditions */
702:     PetscInt    xs, xm, ys, ym;
703:     PetscInt    mx = user->mx, my = user->my;
704:     PetscReal **x;

706:     /* Get local mesh boundaries */
707:     DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL);

709:     /* Get pointers to vector data */
710:     DMDAVecGetArray(user->dm, X, (void **)&x);

712:     /* Perform local computations */
713:     for (j = ys; j < ys + ym; j++) {
714:       for (i = xs; i < xs + xm; i++) x[j][i] = (((j + 1) * user->bottom[i - xs + 1] + (my - j + 1) * user->top[i - xs + 1]) / (my + 2) + ((i + 1) * user->left[j - ys + 1] + (mx - i + 1) * user->right[j - ys + 1]) / (mx + 2)) / 2.0;
715:     }
716:     DMDAVecRestoreArray(user->dm, X, (void **)&x);
717:     PetscLogFlops(9.0 * xm * ym);
718:   }
719:   return 0;
720: }

722: /*-----------------------------------------------------------------------*/
723: PetscErrorCode My_Monitor(Tao tao, void *ctx)
724: {
725:   Vec X;

727:   TaoGetSolution(tao, &X);
728:   VecView(X, PETSC_VIEWER_STDOUT_WORLD);
729:   return 0;
730: }

732: /*TEST

734:    build:
735:       requires: !complex

737:    test:
738:       args: -tao_smonitor -tao_type lmvm -mx 10 -my 8 -tao_gatol 1.e-3
739:       requires: !single

741:    test:
742:       suffix: 2
743:       nsize: 2
744:       args: -tao_smonitor -tao_type nls -tao_nls_ksp_max_it 15 -tao_gatol 1.e-4
745:       filter: grep -v "nls ksp"
746:       requires: !single

748:    test:
749:       suffix: 3
750:       nsize: 3
751:       args: -tao_smonitor -tao_type cg -tao_cg_type fr -mx 10 -my 10 -tao_gatol 1.e-3
752:       requires: !single

754:    test:
755:       suffix: 5
756:       nsize: 2
757:       args: -tao_smonitor -tao_type bmrm -mx 10 -my 8 -tao_gatol 1.e-3
758:       requires: !single

760: TEST*/