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, (char *)0, 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 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: PetscCall(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: PetscCall(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
479: . flg - flag indicating matrix structure
481: */
482: PetscErrorCode FormHessian(Tao tao, Vec X, Mat H, Mat Hpre, void *ptr)
483: {
484: AppCtx *user = (AppCtx *)ptr;
486: PetscFunctionBegin;
487: /* Evaluate the Hessian entries*/
488: PetscCall(QuadraticH(user, X, H));
489: PetscFunctionReturn(PETSC_SUCCESS);
490: }
492: /* ------------------------------------------------------------------- */
493: /*
494: QuadraticH - Evaluates Hessian matrix.
496: Input Parameters:
497: . user - user-defined context, as set by TaoSetHessian()
498: . X - input vector
500: Output Parameter:
501: . H - Hessian matrix
502: */
503: PetscErrorCode QuadraticH(AppCtx *user, Vec X, Mat Hessian)
504: {
505: PetscInt i, j, k;
506: PetscInt mx = user->mx, my = user->my;
507: PetscInt xs, xm, gxs, gxm, ys, ym, gys, gym;
508: PetscReal hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), hydhx = hy / hx, hxdhy = hx / hy;
509: PetscReal f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb;
510: PetscReal hl, hr, ht, hb, hc, htl, hbr;
511: PetscReal **x, v[7];
512: MatStencil col[7], row;
513: Vec localX;
514: PetscBool assembled;
516: PetscFunctionBegin;
517: /* Get local mesh boundaries */
518: PetscCall(DMGetLocalVector(user->dm, &localX));
520: PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
521: PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL));
523: /* Scatter ghost points to local vector */
524: PetscCall(DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, localX));
525: PetscCall(DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, localX));
527: /* Get pointers to vector data */
528: PetscCall(DMDAVecGetArray(user->dm, localX, (void **)&x));
530: /* Initialize matrix entries to zero */
531: PetscCall(MatAssembled(Hessian, &assembled));
532: if (assembled) PetscCall(MatZeroEntries(Hessian));
534: /* Set various matrix options */
535: PetscCall(MatSetOption(Hessian, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE));
537: /* Compute Hessian over the locally owned part of the mesh */
538: for (j = ys; j < ys + ym; j++) {
539: for (i = xs; i < xs + xm; i++) {
540: xc = x[j][i];
541: xlt = xrb = xl = xr = xb = xt = xc;
543: /* Left side */
544: if (i == 0) {
545: xl = user->left[j - ys + 1];
546: xlt = user->left[j - ys + 2];
547: } else {
548: xl = x[j][i - 1];
549: }
551: if (j == 0) {
552: xb = user->bottom[i - xs + 1];
553: xrb = user->bottom[i - xs + 2];
554: } else {
555: xb = x[j - 1][i];
556: }
558: if (i + 1 == mx) {
559: xr = user->right[j - ys + 1];
560: xrb = user->right[j - ys];
561: } else {
562: xr = x[j][i + 1];
563: }
565: if (j + 1 == my) {
566: xt = user->top[i - xs + 1];
567: xlt = user->top[i - xs];
568: } else {
569: xt = x[j + 1][i];
570: }
572: if (i > 0 && j + 1 < my) xlt = x[j + 1][i - 1];
573: if (j > 0 && i + 1 < mx) xrb = x[j - 1][i + 1];
575: d1 = (xc - xl) / hx;
576: d2 = (xc - xr) / hx;
577: d3 = (xc - xt) / hy;
578: d4 = (xc - xb) / hy;
579: d5 = (xrb - xr) / hy;
580: d6 = (xrb - xb) / hx;
581: d7 = (xlt - xl) / hy;
582: d8 = (xlt - xt) / hx;
584: f1 = PetscSqrtReal(1.0 + d1 * d1 + d7 * d7);
585: f2 = PetscSqrtReal(1.0 + d1 * d1 + d4 * d4);
586: f3 = PetscSqrtReal(1.0 + d3 * d3 + d8 * d8);
587: f4 = PetscSqrtReal(1.0 + d3 * d3 + d2 * d2);
588: f5 = PetscSqrtReal(1.0 + d2 * d2 + d5 * d5);
589: f6 = PetscSqrtReal(1.0 + d4 * d4 + d6 * d6);
591: hl = (-hydhx * (1.0 + d7 * d7) + d1 * d7) / (f1 * f1 * f1) + (-hydhx * (1.0 + d4 * d4) + d1 * d4) / (f2 * f2 * f2);
592: hr = (-hydhx * (1.0 + d5 * d5) + d2 * d5) / (f5 * f5 * f5) + (-hydhx * (1.0 + d3 * d3) + d2 * d3) / (f4 * f4 * f4);
593: ht = (-hxdhy * (1.0 + d8 * d8) + d3 * d8) / (f3 * f3 * f3) + (-hxdhy * (1.0 + d2 * d2) + d2 * d3) / (f4 * f4 * f4);
594: hb = (-hxdhy * (1.0 + d6 * d6) + d4 * d6) / (f6 * f6 * f6) + (-hxdhy * (1.0 + d1 * d1) + d1 * d4) / (f2 * f2 * f2);
596: hbr = -d2 * d5 / (f5 * f5 * f5) - d4 * d6 / (f6 * f6 * f6);
597: htl = -d1 * d7 / (f1 * f1 * f1) - d3 * d8 / (f3 * f3 * f3);
599: 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);
601: hl /= 2.0;
602: hr /= 2.0;
603: ht /= 2.0;
604: hb /= 2.0;
605: hbr /= 2.0;
606: htl /= 2.0;
607: hc /= 2.0;
609: row.j = j;
610: row.i = i;
611: k = 0;
612: if (j > 0) {
613: v[k] = hb;
614: col[k].j = j - 1;
615: col[k].i = i;
616: k++;
617: }
619: if (j > 0 && i < mx - 1) {
620: v[k] = hbr;
621: col[k].j = j - 1;
622: col[k].i = i + 1;
623: k++;
624: }
626: if (i > 0) {
627: v[k] = hl;
628: col[k].j = j;
629: col[k].i = i - 1;
630: k++;
631: }
633: v[k] = hc;
634: col[k].j = j;
635: col[k].i = i;
636: k++;
638: if (i < mx - 1) {
639: v[k] = hr;
640: col[k].j = j;
641: col[k].i = i + 1;
642: k++;
643: }
645: if (i > 0 && j < my - 1) {
646: v[k] = htl;
647: col[k].j = j + 1;
648: col[k].i = i - 1;
649: k++;
650: }
652: if (j < my - 1) {
653: v[k] = ht;
654: col[k].j = j + 1;
655: col[k].i = i;
656: k++;
657: }
659: /*
660: Set matrix values using local numbering, which was defined
661: earlier, in the main routine.
662: */
663: PetscCall(MatSetValuesStencil(Hessian, 1, &row, k, col, v, INSERT_VALUES));
664: }
665: }
667: PetscCall(DMDAVecRestoreArray(user->dm, localX, (void **)&x));
668: PetscCall(DMRestoreLocalVector(user->dm, &localX));
670: PetscCall(MatAssemblyBegin(Hessian, MAT_FINAL_ASSEMBLY));
671: PetscCall(MatAssemblyEnd(Hessian, MAT_FINAL_ASSEMBLY));
673: PetscCall(PetscLogFlops(199.0 * xm * ym));
674: PetscFunctionReturn(PETSC_SUCCESS);
675: }
677: /* ------------------------------------------------------------------- */
678: /*
679: MSA_BoundaryConditions - Calculates the boundary conditions for
680: the region.
682: Input Parameter:
683: . user - user-defined application context
685: Output Parameter:
686: . user - user-defined application context
687: */
688: static PetscErrorCode MSA_BoundaryConditions(AppCtx *user)
689: {
690: PetscInt i, j, k, limit = 0, maxits = 5;
691: PetscInt xs, ys, xm, ym, gxs, gys, gxm, gym;
692: PetscInt mx = user->mx, my = user->my;
693: PetscInt bsize = 0, lsize = 0, tsize = 0, rsize = 0;
694: PetscReal one = 1.0, two = 2.0, three = 3.0, tol = 1e-10;
695: PetscReal fnorm, det, hx, hy, xt = 0, yt = 0;
696: PetscReal u1, u2, nf1, nf2, njac11, njac12, njac21, njac22;
697: PetscReal b = -0.5, t = 0.5, l = -0.5, r = 0.5;
698: PetscReal *boundary;
699: PetscBool flg;
701: PetscFunctionBegin;
702: /* Get local mesh boundaries */
703: PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
704: PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL));
706: bsize = xm + 2;
707: lsize = ym + 2;
708: rsize = ym + 2;
709: tsize = xm + 2;
711: PetscCall(PetscMalloc1(bsize, &user->bottom));
712: PetscCall(PetscMalloc1(tsize, &user->top));
713: PetscCall(PetscMalloc1(lsize, &user->left));
714: PetscCall(PetscMalloc1(rsize, &user->right));
716: hx = (r - l) / (mx + 1);
717: hy = (t - b) / (my + 1);
719: for (j = 0; j < 4; j++) {
720: if (j == 0) {
721: yt = b;
722: xt = l + hx * xs;
723: limit = bsize;
724: boundary = user->bottom;
725: } else if (j == 1) {
726: yt = t;
727: xt = l + hx * xs;
728: limit = tsize;
729: boundary = user->top;
730: } else if (j == 2) {
731: yt = b + hy * ys;
732: xt = l;
733: limit = lsize;
734: boundary = user->left;
735: } else { /* if (j==3) */
736: yt = b + hy * ys;
737: xt = r;
738: limit = rsize;
739: boundary = user->right;
740: }
742: for (i = 0; i < limit; i++) {
743: u1 = xt;
744: u2 = -yt;
745: for (k = 0; k < maxits; k++) {
746: nf1 = u1 + u1 * u2 * u2 - u1 * u1 * u1 / three - xt;
747: nf2 = -u2 - u1 * u1 * u2 + u2 * u2 * u2 / three - yt;
748: fnorm = PetscSqrtReal(nf1 * nf1 + nf2 * nf2);
749: if (fnorm <= tol) break;
750: njac11 = one + u2 * u2 - u1 * u1;
751: njac12 = two * u1 * u2;
752: njac21 = -two * u1 * u2;
753: njac22 = -one - u1 * u1 + u2 * u2;
754: det = njac11 * njac22 - njac21 * njac12;
755: u1 = u1 - (njac22 * nf1 - njac12 * nf2) / det;
756: u2 = u2 - (njac11 * nf2 - njac21 * nf1) / det;
757: }
759: boundary[i] = u1 * u1 - u2 * u2;
760: if (j == 0 || j == 1) {
761: xt = xt + hx;
762: } else { /* if (j==2 || j==3) */
763: yt = yt + hy;
764: }
765: }
766: }
768: /* Scale the boundary if desired */
769: if (1 == 1) {
770: PetscReal scl = 1.0;
772: PetscCall(PetscOptionsGetReal(NULL, NULL, "-bottom", &scl, &flg));
773: if (flg) {
774: for (i = 0; i < bsize; i++) user->bottom[i] *= scl;
775: }
777: PetscCall(PetscOptionsGetReal(NULL, NULL, "-top", &scl, &flg));
778: if (flg) {
779: for (i = 0; i < tsize; i++) user->top[i] *= scl;
780: }
782: PetscCall(PetscOptionsGetReal(NULL, NULL, "-right", &scl, &flg));
783: if (flg) {
784: for (i = 0; i < rsize; i++) user->right[i] *= scl;
785: }
787: PetscCall(PetscOptionsGetReal(NULL, NULL, "-left", &scl, &flg));
788: if (flg) {
789: for (i = 0; i < lsize; i++) user->left[i] *= scl;
790: }
791: }
792: PetscFunctionReturn(PETSC_SUCCESS);
793: }
795: /* ------------------------------------------------------------------- */
796: /*
797: MSA_InitialPoint - Calculates the initial guess in one of three ways.
799: Input Parameters:
800: . user - user-defined application context
801: . X - vector for initial guess
803: Output Parameters:
804: . X - newly computed initial guess
805: */
806: static PetscErrorCode MSA_InitialPoint(AppCtx *user, Vec X)
807: {
808: PetscInt start2 = -1, i, j;
809: PetscReal start1 = 0;
810: PetscBool flg1, flg2;
812: PetscFunctionBegin;
813: PetscCall(PetscOptionsGetReal(NULL, NULL, "-start", &start1, &flg1));
814: PetscCall(PetscOptionsGetInt(NULL, NULL, "-random", &start2, &flg2));
816: if (flg1) { /* The zero vector is reasonable */
817: PetscCall(VecSet(X, start1));
818: } else if (flg2 && start2 > 0) { /* Try a random start between -0.5 and 0.5 */
819: PetscRandom rctx;
820: PetscReal np5 = -0.5;
822: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx));
823: PetscCall(VecSetRandom(X, rctx));
824: PetscCall(PetscRandomDestroy(&rctx));
825: PetscCall(VecShift(X, np5));
826: } else { /* Take an average of the boundary conditions */
827: PetscInt xs, xm, ys, ym;
828: PetscInt mx = user->mx, my = user->my;
829: PetscReal **x;
831: /* Get local mesh boundaries */
832: PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
834: /* Get pointers to vector data */
835: PetscCall(DMDAVecGetArray(user->dm, X, (void **)&x));
837: /* Perform local computations */
838: for (j = ys; j < ys + ym; j++) {
839: 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;
840: }
841: PetscCall(DMDAVecRestoreArray(user->dm, X, (void **)&x));
842: PetscCall(PetscLogFlops(9.0 * xm * ym));
843: }
844: PetscFunctionReturn(PETSC_SUCCESS);
845: }
847: /*-----------------------------------------------------------------------*/
848: PetscErrorCode My_Monitor(Tao tao, void *ctx)
849: {
850: Vec X;
852: PetscFunctionBegin;
853: PetscCall(TaoGetSolution(tao, &X));
854: PetscCall(VecView(X, PETSC_VIEWER_STDOUT_WORLD));
855: PetscFunctionReturn(PETSC_SUCCESS);
856: }
858: /*TEST
860: build:
861: requires: !complex
863: test:
864: args: -tao_monitor_short -tao_type lmvm -da_grid_x 10 -da_grid_y 8 -tao_gatol 1.e-3
865: requires: !single
867: test:
868: suffix: 2
869: nsize: 2
870: args: -tao_monitor_short -tao_type nls -tao_nls_ksp_max_it 15 -tao_gatol 1.e-4
871: filter: grep -v "nls ksp"
872: requires: !single
874: test:
875: suffix: 2_snes
876: nsize: 2
877: args: -tao_monitor_short -tao_type snes -ksp_converged_maxits -ksp_max_it 15 -snes_atol 1.e-4
878: filter: grep -v "nls ksp"
879: requires: !single
881: test:
882: suffix: 3
883: nsize: 3
884: args: -tao_monitor_short -tao_type cg -tao_cg_type fr -da_grid_x 10 -da_grid_y 10 -tao_gatol 1.e-3
885: requires: !single
887: test:
888: suffix: 3_snes
889: nsize: 3
890: 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
891: requires: !single
893: test:
894: suffix: 4_snes_ngmres
895: 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}
896: requires: !single
898: test:
899: suffix: 5
900: nsize: 2
901: args: -tao_monitor_short -tao_type bmrm -da_grid_x 10 -da_grid_y 8 -tao_gatol 1.e-3
902: requires: !single
904: TEST*/