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: -da_grid_x <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
17: -da_grid_y <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, FormFunction(),
24: FormFunctionGradient(), 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: static PetscErrorCode QuadraticH(AppCtx *, Vec, Mat);
38: static PetscErrorCode FormFunction(Tao, Vec, PetscReal *, void *);
39: static PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);
40: static PetscErrorCode FormGradient(Tao, Vec, Vec, void *);
41: static PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void *);
42: static PetscErrorCode My_Monitor(Tao, void *);
44: int main(int argc, char **argv)
45: {
46: Vec x; /* solution, gradient vectors */
47: PetscBool 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 */
55: PetscFunctionBeginUser;
56: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
58: /* Create distributed array (DM) to manage parallel grid and vectors */
59: PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 10, 10, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &user.dm));
60: PetscCall(DMSetFromOptions(user.dm));
61: PetscCall(DMSetUp(user.dm));
62: PetscCall(DMDASetUniformCoordinates(user.dm, -0.5, 0.5, -0.5, 0.5, PETSC_DECIDE, PETSC_DECIDE));
63: PetscCall(DMDAGetInfo(user.dm, PETSC_IGNORE, &user.mx, &user.my, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE));
65: PetscCall(PetscPrintf(MPI_COMM_WORLD, "\n---- Minimum Surface Area Problem -----\n"));
66: PetscCall(PetscPrintf(MPI_COMM_WORLD, "mx: %" PetscInt_FMT " my: %" PetscInt_FMT " \n\n", user.mx, user.my));
68: /* Create TAO solver and set desired solution method.*/
69: PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao));
70: PetscCall(TaoSetType(tao, TAOCG));
72: /*
73: Extract global vector from DA for the vector of variables -- PETSC routine
74: Compute the initial solution -- application specific, see below
75: Set this vector for use by TAO -- TAO routine
76: */
77: PetscCall(DMCreateGlobalVector(user.dm, &x));
78: PetscCall(MSA_BoundaryConditions(&user));
79: PetscCall(MSA_InitialPoint(&user, x));
80: PetscCall(TaoSetSolution(tao, x));
82: /*
83: Initialize the Application context for use in function evaluations -- application specific, see below.
84: Set routines for function and gradient evaluation
85: */
86: PetscCall(TaoSetObjective(tao, FormFunction, (void *)&user));
87: PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user));
89: /*
90: Given the command line arguments, calculate the hessian with either the user-
91: provided function FormHessian, or the default finite-difference driven Hessian
92: functions
93: */
94: PetscCall(PetscOptionsHasName(NULL, NULL, "-tao_fddefault", &fddefault));
95: PetscCall(PetscOptionsHasName(NULL, NULL, "-tao_fdcoloring", &fdcoloring));
97: /*
98: Create a matrix data structure to store the Hessian and set
99: the Hessian evaluation routine.
100: Set the matrix nonzero structure to be used for Hessian evaluations
101: */
102: PetscCall(DMCreateMatrix(user.dm, &user.H));
103: PetscCall(MatSetOption(user.H, MAT_SYMMETRIC, PETSC_TRUE));
105: if (fdcoloring) {
106: PetscCall(DMCreateColoring(user.dm, IS_COLORING_GLOBAL, &iscoloring));
107: PetscCall(MatFDColoringCreate(user.H, iscoloring, &matfdcoloring));
108: PetscCall(ISColoringDestroy(&iscoloring));
109: PetscCall(MatFDColoringSetFunction(matfdcoloring, (PetscErrorCode (*)(void))FormGradient, (void *)&user));
110: PetscCall(MatFDColoringSetFromOptions(matfdcoloring));
111: PetscCall(TaoSetHessian(tao, user.H, user.H, TaoDefaultComputeHessianColor, (void *)matfdcoloring));
112: } else if (fddefault) {
113: PetscCall(TaoSetHessian(tao, user.H, user.H, TaoDefaultComputeHessian, (void *)NULL));
114: } else {
115: PetscCall(TaoSetHessian(tao, user.H, user.H, FormHessian, (void *)&user));
116: }
118: /*
119: If my_monitor option is in command line, then use the user-provided
120: monitoring function
121: */
122: PetscCall(PetscOptionsHasName(NULL, NULL, "-my_monitor", &viewmat));
123: if (viewmat) PetscCall(TaoMonitorSet(tao, My_Monitor, NULL, NULL));
125: /* Check for any tao command line options */
126: PetscCall(TaoSetFromOptions(tao));
128: /* SOLVE THE APPLICATION */
129: PetscCall(TaoSolve(tao));
131: PetscCall(TaoView(tao, PETSC_VIEWER_STDOUT_WORLD));
133: /* Free TAO data structures */
134: PetscCall(TaoDestroy(&tao));
136: /* Free PETSc data structures */
137: PetscCall(VecDestroy(&x));
138: PetscCall(MatDestroy(&user.H));
139: if (fdcoloring) PetscCall(MatFDColoringDestroy(&matfdcoloring));
140: PetscCall(PetscFree(user.bottom));
141: PetscCall(PetscFree(user.top));
142: PetscCall(PetscFree(user.left));
143: PetscCall(PetscFree(user.right));
144: PetscCall(DMDestroy(&user.dm));
145: PetscCall(PetscFinalize());
146: return 0;
147: }
149: /* -------------------------------------------------------------------- */
150: /*
151: FormFunction - Evaluates the objective function.
153: Input Parameters:
154: . tao - the Tao context
155: . X - input vector
156: . userCtx - optional user-defined context, as set by TaoSetObjective()
158: Output Parameters:
159: . fcn - the newly evaluated function
160: */
161: PetscErrorCode FormFunction(Tao tao, Vec X, PetscReal *fcn, void *userCtx)
162: {
163: AppCtx *user = (AppCtx *)userCtx;
164: PetscInt i, j;
165: PetscInt mx = user->mx, my = user->my;
166: PetscInt xs, xm, gxs, gxm, ys, ym, gys, gym;
167: PetscReal ft = 0.0;
168: PetscReal hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), area = 0.5 * hx * hy;
169: PetscReal rhx = mx + 1, rhy = my + 1;
170: PetscReal f2, f4, d1, d2, d3, d4, xc, xl, xr, xt, xb;
171: PetscReal **x;
172: Vec localX;
174: PetscFunctionBegin;
175: /* Get local mesh boundaries */
176: PetscCall(DMGetLocalVector(user->dm, &localX));
177: PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
178: PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL));
180: /* Scatter ghost points to local vector */
181: PetscCall(DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, localX));
182: PetscCall(DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, localX));
184: /* Get pointers to vector data */
185: PetscCall(DMDAVecGetArray(user->dm, localX, (void **)&x));
187: /* Compute function and gradient over the locally owned part of the mesh */
188: for (j = ys; j < ys + ym; j++) {
189: for (i = xs; i < xs + xm; i++) {
190: xc = x[j][i];
192: if (i == 0) { /* left side */
193: xl = user->left[j - ys + 1];
194: } else {
195: xl = x[j][i - 1];
196: }
198: if (j == 0) { /* bottom side */
199: xb = user->bottom[i - xs + 1];
200: } else {
201: xb = x[j - 1][i];
202: }
204: if (i + 1 == gxs + gxm) { /* right side */
205: xr = user->right[j - ys + 1];
206: } else {
207: xr = x[j][i + 1];
208: }
210: if (j + 1 == gys + gym) { /* top side */
211: xt = user->top[i - xs + 1];
212: } else {
213: xt = x[j + 1][i];
214: }
216: d1 = (xc - xl);
217: d2 = (xc - xr);
218: d3 = (xc - xt);
219: d4 = (xc - xb);
221: d1 *= rhx;
222: d2 *= rhx;
223: d3 *= rhy;
224: d4 *= rhy;
226: f2 = PetscSqrtReal(1.0 + d1 * d1 + d4 * d4);
227: f4 = PetscSqrtReal(1.0 + d3 * d3 + d2 * d2);
229: ft = ft + (f2 + f4);
230: }
231: }
233: /* Compute triangular areas along the border of the domain. */
234: if (xs == 0) { /* left side */
235: for (j = ys; j < ys + ym; j++) {
236: d3 = (user->left[j - ys + 1] - user->left[j - ys + 2]) * rhy;
237: d2 = (user->left[j - ys + 1] - x[j][0]) * rhx;
238: ft = ft + PetscSqrtReal(1.0 + d3 * d3 + d2 * d2);
239: }
240: }
241: if (ys == 0) { /* bottom side */
242: for (i = xs; i < xs + xm; i++) {
243: d2 = (user->bottom[i + 1 - xs] - user->bottom[i - xs + 2]) * rhx;
244: d3 = (user->bottom[i - xs + 1] - x[0][i]) * rhy;
245: ft = ft + PetscSqrtReal(1.0 + d3 * d3 + d2 * d2);
246: }
247: }
248: if (xs + xm == mx) { /* right side */
249: for (j = ys; j < ys + ym; j++) {
250: d1 = (x[j][mx - 1] - user->right[j - ys + 1]) * rhx;
251: d4 = (user->right[j - ys] - user->right[j - ys + 1]) * rhy;
252: ft = ft + PetscSqrtReal(1.0 + d1 * d1 + d4 * d4);
253: }
254: }
255: if (ys + ym == my) { /* top side */
256: for (i = xs; i < xs + xm; i++) {
257: d1 = (x[my - 1][i] - user->top[i - xs + 1]) * rhy;
258: d4 = (user->top[i - xs + 1] - user->top[i - xs]) * rhx;
259: ft = ft + PetscSqrtReal(1.0 + d1 * d1 + d4 * d4);
260: }
261: }
262: if (ys == 0 && xs == 0) {
263: d1 = (user->left[0] - user->left[1]) / hy;
264: d2 = (user->bottom[0] - user->bottom[1]) * rhx;
265: ft += PetscSqrtReal(1.0 + d1 * d1 + d2 * d2);
266: }
267: if (ys + ym == my && xs + xm == mx) {
268: d1 = (user->right[ym + 1] - user->right[ym]) * rhy;
269: d2 = (user->top[xm + 1] - user->top[xm]) * rhx;
270: ft += PetscSqrtReal(1.0 + d1 * d1 + d2 * d2);
271: }
273: ft = ft * area;
274: PetscCallMPI(MPIU_Allreduce(&ft, fcn, 1, MPIU_REAL, MPIU_SUM, MPI_COMM_WORLD));
276: /* Restore vectors */
277: PetscCall(DMDAVecRestoreArray(user->dm, localX, (void **)&x));
278: PetscCall(DMRestoreLocalVector(user->dm, &localX));
279: PetscFunctionReturn(PETSC_SUCCESS);
280: }
282: /* -------------------------------------------------------------------- */
283: /*
284: FormFunctionGradient - Evaluates the function and corresponding gradient.
286: Input Parameters:
287: . tao - the Tao context
288: . X - input vector
289: . userCtx - optional user-defined context, as set by TaoSetObjectiveAndGradient()
291: Output Parameters:
292: . fcn - the newly evaluated function
293: . G - vector containing the newly evaluated gradient
294: */
295: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn, Vec G, void *userCtx)
296: {
297: AppCtx *user = (AppCtx *)userCtx;
298: PetscInt i, j;
299: PetscInt mx = user->mx, my = user->my;
300: PetscInt xs, xm, gxs, gxm, ys, ym, gys, gym;
301: PetscReal ft = 0.0;
302: PetscReal hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), hydhx = hy / hx, hxdhy = hx / hy, area = 0.5 * hx * hy;
303: PetscReal rhx = mx + 1, rhy = my + 1;
304: PetscReal f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb;
305: PetscReal df1dxc, df2dxc, df3dxc, df4dxc, df5dxc, df6dxc;
306: PetscReal **g, **x;
307: Vec localX;
309: PetscFunctionBegin;
310: /* Get local mesh boundaries */
311: PetscCall(DMGetLocalVector(user->dm, &localX));
312: PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
313: PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL));
315: /* Scatter ghost points to local vector */
316: PetscCall(DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, localX));
317: PetscCall(DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, localX));
319: /* Get pointers to vector data */
320: PetscCall(DMDAVecGetArray(user->dm, localX, (void **)&x));
321: PetscCall(DMDAVecGetArray(user->dm, G, (void **)&g));
323: /* Compute function and gradient over the locally owned part of the mesh */
324: for (j = ys; j < ys + ym; j++) {
325: for (i = xs; i < xs + xm; i++) {
326: xc = x[j][i];
327: xlt = xrb = xl = xr = xb = xt = xc;
329: if (i == 0) { /* left side */
330: xl = user->left[j - ys + 1];
331: xlt = user->left[j - ys + 2];
332: } else {
333: xl = x[j][i - 1];
334: }
336: if (j == 0) { /* bottom side */
337: xb = user->bottom[i - xs + 1];
338: xrb = user->bottom[i - xs + 2];
339: } else {
340: xb = x[j - 1][i];
341: }
343: if (i + 1 == gxs + gxm) { /* right side */
344: xr = user->right[j - ys + 1];
345: xrb = user->right[j - ys];
346: } else {
347: xr = x[j][i + 1];
348: }
350: if (j + 1 == gys + gym) { /* top side */
351: xt = user->top[i - xs + 1];
352: xlt = user->top[i - xs];
353: } else {
354: xt = x[j + 1][i];
355: }
357: if (i > gxs && j + 1 < gys + gym) xlt = x[j + 1][i - 1];
358: if (j > gys && i + 1 < gxs + gxm) xrb = x[j - 1][i + 1];
360: d1 = (xc - xl);
361: d2 = (xc - xr);
362: d3 = (xc - xt);
363: d4 = (xc - xb);
364: d5 = (xr - xrb);
365: d6 = (xrb - xb);
366: d7 = (xlt - xl);
367: d8 = (xt - xlt);
369: df1dxc = d1 * hydhx;
370: df2dxc = (d1 * hydhx + d4 * hxdhy);
371: df3dxc = d3 * hxdhy;
372: df4dxc = (d2 * hydhx + d3 * hxdhy);
373: df5dxc = d2 * hydhx;
374: df6dxc = d4 * hxdhy;
376: d1 *= rhx;
377: d2 *= rhx;
378: d3 *= rhy;
379: d4 *= rhy;
380: d5 *= rhy;
381: d6 *= rhx;
382: d7 *= rhy;
383: d8 *= rhx;
385: f1 = PetscSqrtReal(1.0 + d1 * d1 + d7 * d7);
386: f2 = PetscSqrtReal(1.0 + d1 * d1 + d4 * d4);
387: f3 = PetscSqrtReal(1.0 + d3 * d3 + d8 * d8);
388: f4 = PetscSqrtReal(1.0 + d3 * d3 + d2 * d2);
389: f5 = PetscSqrtReal(1.0 + d2 * d2 + d5 * d5);
390: f6 = PetscSqrtReal(1.0 + d4 * d4 + d6 * d6);
392: ft = ft + (f2 + f4);
394: df1dxc /= f1;
395: df2dxc /= f2;
396: df3dxc /= f3;
397: df4dxc /= f4;
398: df5dxc /= f5;
399: df6dxc /= f6;
401: g[j][i] = (df1dxc + df2dxc + df3dxc + df4dxc + df5dxc + df6dxc) * 0.5;
402: }
403: }
405: /* Compute triangular areas along the border of the domain. */
406: if (xs == 0) { /* left side */
407: for (j = ys; j < ys + ym; j++) {
408: d3 = (user->left[j - ys + 1] - user->left[j - ys + 2]) * rhy;
409: d2 = (user->left[j - ys + 1] - x[j][0]) * rhx;
410: ft = ft + PetscSqrtReal(1.0 + d3 * d3 + d2 * d2);
411: }
412: }
413: if (ys == 0) { /* bottom side */
414: for (i = xs; i < xs + xm; i++) {
415: d2 = (user->bottom[i + 1 - xs] - user->bottom[i - xs + 2]) * rhx;
416: d3 = (user->bottom[i - xs + 1] - x[0][i]) * rhy;
417: ft = ft + PetscSqrtReal(1.0 + d3 * d3 + d2 * d2);
418: }
419: }
421: if (xs + xm == mx) { /* right side */
422: for (j = ys; j < ys + ym; j++) {
423: d1 = (x[j][mx - 1] - user->right[j - ys + 1]) * rhx;
424: d4 = (user->right[j - ys] - user->right[j - ys + 1]) * rhy;
425: ft = ft + PetscSqrtReal(1.0 + d1 * d1 + d4 * d4);
426: }
427: }
428: if (ys + ym == my) { /* top side */
429: for (i = xs; i < xs + xm; i++) {
430: d1 = (x[my - 1][i] - user->top[i - xs + 1]) * rhy;
431: d4 = (user->top[i - xs + 1] - user->top[i - xs]) * rhx;
432: ft = ft + PetscSqrtReal(1.0 + d1 * d1 + d4 * d4);
433: }
434: }
436: if (ys == 0 && xs == 0) {
437: d1 = (user->left[0] - user->left[1]) / hy;
438: d2 = (user->bottom[0] - user->bottom[1]) * rhx;
439: ft += PetscSqrtReal(1.0 + d1 * d1 + d2 * d2);
440: }
441: if (ys + ym == my && xs + xm == mx) {
442: d1 = (user->right[ym + 1] - user->right[ym]) * rhy;
443: d2 = (user->top[xm + 1] - user->top[xm]) * rhx;
444: ft += PetscSqrtReal(1.0 + d1 * d1 + d2 * d2);
445: }
447: ft = ft * area;
448: PetscCallMPI(MPIU_Allreduce(&ft, fcn, 1, MPIU_REAL, MPIU_SUM, MPI_COMM_WORLD));
450: /* Restore vectors */
451: PetscCall(DMDAVecRestoreArray(user->dm, localX, (void **)&x));
452: PetscCall(DMDAVecRestoreArray(user->dm, G, (void **)&g));
453: PetscCall(DMRestoreLocalVector(user->dm, &localX));
454: PetscCall(PetscLogFlops(67.0 * xm * ym));
455: PetscFunctionReturn(PETSC_SUCCESS);
456: }
458: PetscErrorCode FormGradient(Tao tao, Vec X, Vec G, void *userCtx)
459: {
460: PetscReal fcn;
462: PetscFunctionBegin;
463: PetscCall(FormFunctionGradient(tao, X, &fcn, G, userCtx));
464: PetscFunctionReturn(PETSC_SUCCESS);
465: }
467: /* ------------------------------------------------------------------- */
468: /*
469: FormHessian - Evaluates Hessian matrix.
471: Input Parameters:
472: . tao - the Tao context
473: . x - input vector
474: . ptr - optional user-defined context, as set by TaoSetHessian()
476: Output Parameters:
477: . H - Hessian matrix
478: . Hpre - optionally different preconditioning matrix
480: */
481: PetscErrorCode FormHessian(Tao tao, Vec X, Mat H, Mat Hpre, void *ptr)
482: {
483: AppCtx *user = (AppCtx *)ptr;
485: PetscFunctionBegin;
486: /* Evaluate the Hessian entries*/
487: PetscCall(QuadraticH(user, X, H));
488: PetscFunctionReturn(PETSC_SUCCESS);
489: }
491: /* ------------------------------------------------------------------- */
492: /*
493: QuadraticH - Evaluates Hessian matrix.
495: Input Parameters:
496: . user - user-defined context, as set by TaoSetHessian()
497: . X - input vector
499: Output Parameter:
500: . H - Hessian matrix
501: */
502: PetscErrorCode QuadraticH(AppCtx *user, Vec X, Mat Hessian)
503: {
504: PetscInt i, j, k;
505: PetscInt mx = user->mx, my = user->my;
506: PetscInt xs, xm, gxs, gxm, ys, ym, gys, gym;
507: PetscReal hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), hydhx = hy / hx, hxdhy = hx / hy;
508: PetscReal f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb;
509: PetscReal hl, hr, ht, hb, hc, htl, hbr;
510: PetscReal **x, v[7];
511: MatStencil col[7], row;
512: Vec localX;
513: PetscBool assembled;
515: PetscFunctionBegin;
516: /* Get local mesh boundaries */
517: PetscCall(DMGetLocalVector(user->dm, &localX));
519: PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
520: PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL));
522: /* Scatter ghost points to local vector */
523: PetscCall(DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, localX));
524: PetscCall(DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, localX));
526: /* Get pointers to vector data */
527: PetscCall(DMDAVecGetArray(user->dm, localX, (void **)&x));
529: /* Initialize matrix entries to zero */
530: PetscCall(MatAssembled(Hessian, &assembled));
531: if (assembled) PetscCall(MatZeroEntries(Hessian));
533: /* Set various matrix options */
534: PetscCall(MatSetOption(Hessian, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE));
536: /* Compute Hessian over the locally owned part of the mesh */
537: for (j = ys; j < ys + ym; j++) {
538: for (i = xs; i < xs + xm; i++) {
539: xc = x[j][i];
540: xlt = xrb = xl = xr = xb = xt = xc;
542: /* Left side */
543: if (i == 0) {
544: xl = user->left[j - ys + 1];
545: xlt = user->left[j - ys + 2];
546: } else {
547: xl = x[j][i - 1];
548: }
550: if (j == 0) {
551: xb = user->bottom[i - xs + 1];
552: xrb = user->bottom[i - xs + 2];
553: } else {
554: xb = x[j - 1][i];
555: }
557: if (i + 1 == mx) {
558: xr = user->right[j - ys + 1];
559: xrb = user->right[j - ys];
560: } else {
561: xr = x[j][i + 1];
562: }
564: if (j + 1 == my) {
565: xt = user->top[i - xs + 1];
566: xlt = user->top[i - xs];
567: } else {
568: xt = x[j + 1][i];
569: }
571: if (i > 0 && j + 1 < my) xlt = x[j + 1][i - 1];
572: if (j > 0 && i + 1 < mx) xrb = x[j - 1][i + 1];
574: d1 = (xc - xl) / hx;
575: d2 = (xc - xr) / hx;
576: d3 = (xc - xt) / hy;
577: d4 = (xc - xb) / hy;
578: d5 = (xrb - xr) / hy;
579: d6 = (xrb - xb) / hx;
580: d7 = (xlt - xl) / hy;
581: d8 = (xlt - xt) / hx;
583: f1 = PetscSqrtReal(1.0 + d1 * d1 + d7 * d7);
584: f2 = PetscSqrtReal(1.0 + d1 * d1 + d4 * d4);
585: f3 = PetscSqrtReal(1.0 + d3 * d3 + d8 * d8);
586: f4 = PetscSqrtReal(1.0 + d3 * d3 + d2 * d2);
587: f5 = PetscSqrtReal(1.0 + d2 * d2 + d5 * d5);
588: f6 = PetscSqrtReal(1.0 + d4 * d4 + d6 * d6);
590: hl = (-hydhx * (1.0 + d7 * d7) + d1 * d7) / (f1 * f1 * f1) + (-hydhx * (1.0 + d4 * d4) + d1 * d4) / (f2 * f2 * f2);
591: hr = (-hydhx * (1.0 + d5 * d5) + d2 * d5) / (f5 * f5 * f5) + (-hydhx * (1.0 + d3 * d3) + d2 * d3) / (f4 * f4 * f4);
592: ht = (-hxdhy * (1.0 + d8 * d8) + d3 * d8) / (f3 * f3 * f3) + (-hxdhy * (1.0 + d2 * d2) + d2 * d3) / (f4 * f4 * f4);
593: hb = (-hxdhy * (1.0 + d6 * d6) + d4 * d6) / (f6 * f6 * f6) + (-hxdhy * (1.0 + d1 * d1) + d1 * d4) / (f2 * f2 * f2);
595: hbr = -d2 * d5 / (f5 * f5 * f5) - d4 * d6 / (f6 * f6 * f6);
596: htl = -d1 * d7 / (f1 * f1 * f1) - d3 * d8 / (f3 * f3 * f3);
598: 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);
600: hl /= 2.0;
601: hr /= 2.0;
602: ht /= 2.0;
603: hb /= 2.0;
604: hbr /= 2.0;
605: htl /= 2.0;
606: hc /= 2.0;
608: row.j = j;
609: row.i = i;
610: k = 0;
611: if (j > 0) {
612: v[k] = hb;
613: col[k].j = j - 1;
614: col[k].i = i;
615: k++;
616: }
618: if (j > 0 && i < mx - 1) {
619: v[k] = hbr;
620: col[k].j = j - 1;
621: col[k].i = i + 1;
622: k++;
623: }
625: if (i > 0) {
626: v[k] = hl;
627: col[k].j = j;
628: col[k].i = i - 1;
629: k++;
630: }
632: v[k] = hc;
633: col[k].j = j;
634: col[k].i = i;
635: k++;
637: if (i < mx - 1) {
638: v[k] = hr;
639: col[k].j = j;
640: col[k].i = i + 1;
641: k++;
642: }
644: if (i > 0 && j < my - 1) {
645: v[k] = htl;
646: col[k].j = j + 1;
647: col[k].i = i - 1;
648: k++;
649: }
651: if (j < my - 1) {
652: v[k] = ht;
653: col[k].j = j + 1;
654: col[k].i = i;
655: k++;
656: }
658: /*
659: Set matrix values using local numbering, which was defined
660: earlier, in the main routine.
661: */
662: PetscCall(MatSetValuesStencil(Hessian, 1, &row, k, col, v, INSERT_VALUES));
663: }
664: }
666: PetscCall(DMDAVecRestoreArray(user->dm, localX, (void **)&x));
667: PetscCall(DMRestoreLocalVector(user->dm, &localX));
669: PetscCall(MatAssemblyBegin(Hessian, MAT_FINAL_ASSEMBLY));
670: PetscCall(MatAssemblyEnd(Hessian, MAT_FINAL_ASSEMBLY));
672: PetscCall(PetscLogFlops(199.0 * xm * ym));
673: PetscFunctionReturn(PETSC_SUCCESS);
674: }
676: /* ------------------------------------------------------------------- */
677: /*
678: MSA_BoundaryConditions - Calculates the boundary conditions for
679: the region.
681: Input Parameter:
682: . user - user-defined application context
684: Output Parameter:
685: . user - user-defined application context
686: */
687: static PetscErrorCode MSA_BoundaryConditions(AppCtx *user)
688: {
689: PetscInt i, j, k, limit = 0, maxits = 5;
690: PetscInt xs, ys, xm, ym, gxs, gys, gxm, gym;
691: PetscInt mx = user->mx, my = user->my;
692: PetscInt bsize = 0, lsize = 0, tsize = 0, rsize = 0;
693: PetscReal one = 1.0, two = 2.0, three = 3.0, tol = 1e-10;
694: PetscReal fnorm, det, hx, hy, xt = 0, yt = 0;
695: PetscReal u1, u2, nf1, nf2, njac11, njac12, njac21, njac22;
696: PetscReal b = -0.5, t = 0.5, l = -0.5, r = 0.5;
697: PetscReal *boundary;
698: PetscBool flg;
700: PetscFunctionBegin;
701: /* Get local mesh boundaries */
702: PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
703: PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL));
705: bsize = xm + 2;
706: lsize = ym + 2;
707: rsize = ym + 2;
708: tsize = xm + 2;
710: PetscCall(PetscMalloc1(bsize, &user->bottom));
711: PetscCall(PetscMalloc1(tsize, &user->top));
712: PetscCall(PetscMalloc1(lsize, &user->left));
713: PetscCall(PetscMalloc1(rsize, &user->right));
715: hx = (r - l) / (mx + 1);
716: hy = (t - b) / (my + 1);
718: for (j = 0; j < 4; j++) {
719: if (j == 0) {
720: yt = b;
721: xt = l + hx * xs;
722: limit = bsize;
723: boundary = user->bottom;
724: } else if (j == 1) {
725: yt = t;
726: xt = l + hx * xs;
727: limit = tsize;
728: boundary = user->top;
729: } else if (j == 2) {
730: yt = b + hy * ys;
731: xt = l;
732: limit = lsize;
733: boundary = user->left;
734: } else { /* if (j==3) */
735: yt = b + hy * ys;
736: xt = r;
737: limit = rsize;
738: boundary = user->right;
739: }
741: for (i = 0; i < limit; i++) {
742: u1 = xt;
743: u2 = -yt;
744: for (k = 0; k < maxits; k++) {
745: nf1 = u1 + u1 * u2 * u2 - u1 * u1 * u1 / three - xt;
746: nf2 = -u2 - u1 * u1 * u2 + u2 * u2 * u2 / three - yt;
747: fnorm = PetscSqrtReal(nf1 * nf1 + nf2 * nf2);
748: if (fnorm <= tol) break;
749: njac11 = one + u2 * u2 - u1 * u1;
750: njac12 = two * u1 * u2;
751: njac21 = -two * u1 * u2;
752: njac22 = -one - u1 * u1 + u2 * u2;
753: det = njac11 * njac22 - njac21 * njac12;
754: u1 = u1 - (njac22 * nf1 - njac12 * nf2) / det;
755: u2 = u2 - (njac11 * nf2 - njac21 * nf1) / det;
756: }
758: boundary[i] = u1 * u1 - u2 * u2;
759: if (j == 0 || j == 1) {
760: xt = xt + hx;
761: } else { /* if (j==2 || j==3) */
762: yt = yt + hy;
763: }
764: }
765: }
767: /* Scale the boundary if desired */
768: if (1 == 1) {
769: PetscReal scl = 1.0;
771: PetscCall(PetscOptionsGetReal(NULL, NULL, "-bottom", &scl, &flg));
772: if (flg) {
773: for (i = 0; i < bsize; i++) user->bottom[i] *= scl;
774: }
776: PetscCall(PetscOptionsGetReal(NULL, NULL, "-top", &scl, &flg));
777: if (flg) {
778: for (i = 0; i < tsize; i++) user->top[i] *= scl;
779: }
781: PetscCall(PetscOptionsGetReal(NULL, NULL, "-right", &scl, &flg));
782: if (flg) {
783: for (i = 0; i < rsize; i++) user->right[i] *= scl;
784: }
786: PetscCall(PetscOptionsGetReal(NULL, NULL, "-left", &scl, &flg));
787: if (flg) {
788: for (i = 0; i < lsize; i++) user->left[i] *= scl;
789: }
790: }
791: PetscFunctionReturn(PETSC_SUCCESS);
792: }
794: /* ------------------------------------------------------------------- */
795: /*
796: MSA_InitialPoint - Calculates the initial guess in one of three ways.
798: Input Parameters:
799: . user - user-defined application context
800: . X - vector for initial guess
802: Output Parameters:
803: . X - newly computed initial guess
804: */
805: static PetscErrorCode MSA_InitialPoint(AppCtx *user, Vec X)
806: {
807: PetscInt start2 = -1, i, j;
808: PetscReal start1 = 0;
809: PetscBool flg1, flg2;
811: PetscFunctionBegin;
812: PetscCall(PetscOptionsGetReal(NULL, NULL, "-start", &start1, &flg1));
813: PetscCall(PetscOptionsGetInt(NULL, NULL, "-random", &start2, &flg2));
815: if (flg1) { /* The zero vector is reasonable */
816: PetscCall(VecSet(X, start1));
817: } else if (flg2 && start2 > 0) { /* Try a random start between -0.5 and 0.5 */
818: PetscRandom rctx;
819: PetscReal np5 = -0.5;
821: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx));
822: PetscCall(VecSetRandom(X, rctx));
823: PetscCall(PetscRandomDestroy(&rctx));
824: PetscCall(VecShift(X, np5));
825: } else { /* Take an average of the boundary conditions */
826: PetscInt xs, xm, ys, ym;
827: PetscInt mx = user->mx, my = user->my;
828: PetscReal **x;
830: /* Get local mesh boundaries */
831: PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
833: /* Get pointers to vector data */
834: PetscCall(DMDAVecGetArray(user->dm, X, (void **)&x));
836: /* Perform local computations */
837: for (j = ys; j < ys + ym; j++) {
838: 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;
839: }
840: PetscCall(DMDAVecRestoreArray(user->dm, X, (void **)&x));
841: PetscCall(PetscLogFlops(9.0 * xm * ym));
842: }
843: PetscFunctionReturn(PETSC_SUCCESS);
844: }
846: /*-----------------------------------------------------------------------*/
847: PetscErrorCode My_Monitor(Tao tao, void *ctx)
848: {
849: Vec X;
851: PetscFunctionBegin;
852: PetscCall(TaoGetSolution(tao, &X));
853: PetscCall(VecView(X, PETSC_VIEWER_STDOUT_WORLD));
854: PetscFunctionReturn(PETSC_SUCCESS);
855: }
857: /*TEST
859: build:
860: requires: !complex
862: test:
863: args: -tao_monitor_short -tao_type lmvm -da_grid_x 10 -da_grid_y 8 -tao_gatol 1.e-3
864: requires: !single
866: test:
867: suffix: 2
868: nsize: 2
869: args: -tao_monitor_short -tao_type nls -tao_nls_ksp_max_it 15 -tao_gatol 1.e-4
870: filter: grep -v "nls ksp"
871: requires: !single
873: test:
874: suffix: 2_snes
875: nsize: 2
876: args: -tao_monitor_short -tao_type snes -ksp_converged_maxits -ksp_max_it 15 -snes_atol 1.e-4
877: filter: grep -v "nls ksp"
878: requires: !single
880: test:
881: suffix: 3
882: nsize: 3
883: args: -tao_monitor_short -tao_type cg -tao_cg_type fr -da_grid_x 10 -da_grid_y 10 -tao_gatol 1.e-3
884: requires: !single
886: test:
887: suffix: 3_snes
888: nsize: 3
889: args: -tao_monitor_short -tao_type snes -snes_type ncg -snes_ncg_type fr -da_grid_x 10 -da_grid_y 10 -snes_atol 1.e-4
890: requires: !single
892: test:
893: suffix: 4_snes_ngmres
894: args: -tao_type snes -snes_type ngmres -npc_snes_type ncg -da_grid_x 10 -da_grid_y 10 -snes_atol 1.e-4 -npc_snes_ncg_type fr -snes_converged_reason -snes_monitor ::ascii_info_detail -snes_ngmres_monitor -snes_ngmres_select_type {{linesearch difference}separate output}
895: requires: !single
897: test:
898: suffix: 5
899: nsize: 2
900: args: -tao_monitor_short -tao_type bmrm -da_grid_x 10 -da_grid_y 8 -tao_gatol 1.e-3
901: requires: !single
903: TEST*/