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(TaoSetDM(tao, user.dm));
71: PetscCall(TaoSetType(tao, TAOCG));
73: /*
74: Extract global vector from DA for the vector of variables -- PETSc routine
75: Compute the initial solution -- application specific, see below
76: Set this vector for use by TAO -- TAO routine
77: */
78: PetscCall(DMCreateGlobalVector(user.dm, &x));
79: PetscCall(MSA_BoundaryConditions(&user));
80: PetscCall(MSA_InitialPoint(&user, x));
81: PetscCall(TaoSetSolution(tao, x));
83: /*
84: Initialize the Application context for use in function evaluations -- application specific, see below.
85: Set routines for function and gradient evaluation
86: */
87: PetscCall(TaoSetObjective(tao, FormFunction, (void *)&user));
88: PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user));
90: /*
91: Given the command line arguments, calculate the hessian with either the user-
92: provided function FormHessian, or the default finite-difference driven Hessian
93: functions
94: */
95: PetscCall(PetscOptionsHasName(NULL, NULL, "-tao_fddefault", &fddefault));
96: PetscCall(PetscOptionsHasName(NULL, NULL, "-tao_fdcoloring", &fdcoloring));
98: /*
99: Create a matrix data structure to store the Hessian and set
100: the Hessian evaluation routine.
101: Set the matrix nonzero structure to be used for Hessian evaluations
102: */
103: PetscCall(DMCreateMatrix(user.dm, &user.H));
104: PetscCall(MatSetOption(user.H, MAT_SYMMETRIC, PETSC_TRUE));
106: if (fdcoloring) {
107: PetscCall(DMCreateColoring(user.dm, IS_COLORING_GLOBAL, &iscoloring));
108: PetscCall(MatFDColoringCreate(user.H, iscoloring, &matfdcoloring));
109: PetscCall(ISColoringDestroy(&iscoloring));
110: PetscCall(MatFDColoringSetFunction(matfdcoloring, (MatFDColoringFn *)FormGradient, (void *)&user));
111: PetscCall(MatFDColoringSetFromOptions(matfdcoloring));
112: PetscCall(TaoSetHessian(tao, user.H, user.H, TaoDefaultComputeHessianColor, (void *)matfdcoloring));
113: } else if (fddefault) {
114: PetscCall(TaoSetHessian(tao, user.H, user.H, TaoDefaultComputeHessian, (void *)NULL));
115: } else {
116: PetscCall(TaoSetHessian(tao, user.H, user.H, FormHessian, (void *)&user));
117: }
119: /*
120: If my_monitor option is in command line, then use the user-provided
121: monitoring function
122: */
123: PetscCall(PetscOptionsHasName(NULL, NULL, "-my_monitor", &viewmat));
124: if (viewmat) PetscCall(TaoMonitorSet(tao, My_Monitor, NULL, NULL));
126: /* Check for any tao command line options */
127: PetscCall(TaoSetFromOptions(tao));
129: /* SOLVE THE APPLICATION */
130: PetscCall(TaoSolve(tao));
132: PetscCall(TaoView(tao, PETSC_VIEWER_STDOUT_WORLD));
134: /* Free TAO data structures */
135: PetscCall(TaoDestroy(&tao));
137: /* Free PETSc data structures */
138: PetscCall(VecDestroy(&x));
139: PetscCall(MatDestroy(&user.H));
140: if (fdcoloring) PetscCall(MatFDColoringDestroy(&matfdcoloring));
141: PetscCall(PetscFree(user.bottom));
142: PetscCall(PetscFree(user.top));
143: PetscCall(PetscFree(user.left));
144: PetscCall(PetscFree(user.right));
145: PetscCall(DMDestroy(&user.dm));
146: PetscCall(PetscFinalize());
147: return 0;
148: }
150: /* -------------------------------------------------------------------- */
151: /*
152: FormFunction - Evaluates the objective function.
154: Input Parameters:
155: . tao - the Tao context
156: . X - input vector
157: . userCtx - optional user-defined context, as set by TaoSetObjective()
159: Output Parameters:
160: . fcn - the newly evaluated function
161: */
162: PetscErrorCode FormFunction(Tao tao, Vec X, PetscReal *fcn, void *userCtx)
163: {
164: AppCtx *user = (AppCtx *)userCtx;
165: PetscInt i, j;
166: PetscInt mx = user->mx, my = user->my;
167: PetscInt xs, xm, gxs, gxm, ys, ym, gys, gym;
168: PetscReal ft = 0.0;
169: PetscReal hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), area = 0.5 * hx * hy;
170: PetscReal rhx = mx + 1, rhy = my + 1;
171: PetscReal f2, f4, d1, d2, d3, d4, xc, xl, xr, xt, xb;
172: PetscReal **x;
173: Vec localX;
175: PetscFunctionBegin;
176: /* Get local mesh boundaries */
177: PetscCall(DMGetLocalVector(user->dm, &localX));
178: PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
179: PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL));
181: /* Scatter ghost points to local vector */
182: PetscCall(DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, localX));
183: PetscCall(DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, localX));
185: /* Get pointers to vector data */
186: PetscCall(DMDAVecGetArray(user->dm, localX, (void **)&x));
188: /* Compute function and gradient over the locally owned part of the mesh */
189: for (j = ys; j < ys + ym; j++) {
190: for (i = xs; i < xs + xm; i++) {
191: xc = x[j][i];
193: if (i == 0) { /* left side */
194: xl = user->left[j - ys + 1];
195: } else {
196: xl = x[j][i - 1];
197: }
199: if (j == 0) { /* bottom side */
200: xb = user->bottom[i - xs + 1];
201: } else {
202: xb = x[j - 1][i];
203: }
205: if (i + 1 == gxs + gxm) { /* right side */
206: xr = user->right[j - ys + 1];
207: } else {
208: xr = x[j][i + 1];
209: }
211: if (j + 1 == gys + gym) { /* top side */
212: xt = user->top[i - xs + 1];
213: } else {
214: xt = x[j + 1][i];
215: }
217: d1 = (xc - xl);
218: d2 = (xc - xr);
219: d3 = (xc - xt);
220: d4 = (xc - xb);
222: d1 *= rhx;
223: d2 *= rhx;
224: d3 *= rhy;
225: d4 *= rhy;
227: f2 = PetscSqrtReal(1.0 + d1 * d1 + d4 * d4);
228: f4 = PetscSqrtReal(1.0 + d3 * d3 + d2 * d2);
230: ft = ft + (f2 + f4);
231: }
232: }
234: /* Compute triangular areas along the border of the domain. */
235: if (xs == 0) { /* left side */
236: for (j = ys; j < ys + ym; j++) {
237: d3 = (user->left[j - ys + 1] - user->left[j - ys + 2]) * rhy;
238: d2 = (user->left[j - ys + 1] - x[j][0]) * rhx;
239: ft = ft + PetscSqrtReal(1.0 + d3 * d3 + d2 * d2);
240: }
241: }
242: if (ys == 0) { /* bottom side */
243: for (i = xs; i < xs + xm; i++) {
244: d2 = (user->bottom[i + 1 - xs] - user->bottom[i - xs + 2]) * rhx;
245: d3 = (user->bottom[i - xs + 1] - x[0][i]) * rhy;
246: ft = ft + PetscSqrtReal(1.0 + d3 * d3 + d2 * d2);
247: }
248: }
249: if (xs + xm == mx) { /* right side */
250: for (j = ys; j < ys + ym; j++) {
251: d1 = (x[j][mx - 1] - user->right[j - ys + 1]) * rhx;
252: d4 = (user->right[j - ys] - user->right[j - ys + 1]) * rhy;
253: ft = ft + PetscSqrtReal(1.0 + d1 * d1 + d4 * d4);
254: }
255: }
256: if (ys + ym == my) { /* top side */
257: for (i = xs; i < xs + xm; i++) {
258: d1 = (x[my - 1][i] - user->top[i - xs + 1]) * rhy;
259: d4 = (user->top[i - xs + 1] - user->top[i - xs]) * rhx;
260: ft = ft + PetscSqrtReal(1.0 + d1 * d1 + d4 * d4);
261: }
262: }
263: if (ys == 0 && xs == 0) {
264: d1 = (user->left[0] - user->left[1]) / hy;
265: d2 = (user->bottom[0] - user->bottom[1]) * rhx;
266: ft += PetscSqrtReal(1.0 + d1 * d1 + d2 * d2);
267: }
268: if (ys + ym == my && xs + xm == mx) {
269: d1 = (user->right[ym + 1] - user->right[ym]) * rhy;
270: d2 = (user->top[xm + 1] - user->top[xm]) * rhx;
271: ft += PetscSqrtReal(1.0 + d1 * d1 + d2 * d2);
272: }
274: ft = ft * area;
275: PetscCallMPI(MPIU_Allreduce(&ft, fcn, 1, MPIU_REAL, MPIU_SUM, MPI_COMM_WORLD));
277: /* Restore vectors */
278: PetscCall(DMDAVecRestoreArray(user->dm, localX, (void **)&x));
279: PetscCall(DMRestoreLocalVector(user->dm, &localX));
280: PetscFunctionReturn(PETSC_SUCCESS);
281: }
283: /* -------------------------------------------------------------------- */
284: /*
285: FormFunctionGradient - Evaluates the function and corresponding gradient.
287: Input Parameters:
288: . tao - the Tao context
289: . X - input vector
290: . userCtx - optional user-defined context, as set by TaoSetObjectiveAndGradient()
292: Output Parameters:
293: . fcn - the newly evaluated function
294: . G - vector containing the newly evaluated gradient
295: */
296: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn, Vec G, void *userCtx)
297: {
298: AppCtx *user = (AppCtx *)userCtx;
299: PetscInt i, j;
300: PetscInt mx = user->mx, my = user->my;
301: PetscInt xs, xm, gxs, gxm, ys, ym, gys, gym;
302: PetscReal ft = 0.0;
303: PetscReal hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), hydhx = hy / hx, hxdhy = hx / hy, area = 0.5 * hx * hy;
304: PetscReal rhx = mx + 1, rhy = my + 1;
305: PetscReal f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb;
306: PetscReal df1dxc, df2dxc, df3dxc, df4dxc, df5dxc, df6dxc;
307: PetscReal **g, **x;
308: Vec localX;
310: PetscFunctionBegin;
311: /* Get local mesh boundaries */
312: PetscCall(DMGetLocalVector(user->dm, &localX));
313: PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
314: PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL));
316: /* Scatter ghost points to local vector */
317: PetscCall(DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, localX));
318: PetscCall(DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, localX));
320: /* Get pointers to vector data */
321: PetscCall(DMDAVecGetArray(user->dm, localX, (void **)&x));
322: PetscCall(DMDAVecGetArray(user->dm, G, (void **)&g));
324: /* Compute function and gradient over the locally owned part of the mesh */
325: for (j = ys; j < ys + ym; j++) {
326: for (i = xs; i < xs + xm; i++) {
327: xc = x[j][i];
328: xlt = xrb = xl = xr = xb = xt = xc;
330: if (i == 0) { /* left side */
331: xl = user->left[j - ys + 1];
332: xlt = user->left[j - ys + 2];
333: } else {
334: xl = x[j][i - 1];
335: }
337: if (j == 0) { /* bottom side */
338: xb = user->bottom[i - xs + 1];
339: xrb = user->bottom[i - xs + 2];
340: } else {
341: xb = x[j - 1][i];
342: }
344: if (i + 1 == gxs + gxm) { /* right side */
345: xr = user->right[j - ys + 1];
346: xrb = user->right[j - ys];
347: } else {
348: xr = x[j][i + 1];
349: }
351: if (j + 1 == gys + gym) { /* top side */
352: xt = user->top[i - xs + 1];
353: xlt = user->top[i - xs];
354: } else {
355: xt = x[j + 1][i];
356: }
358: if (i > gxs && j + 1 < gys + gym) xlt = x[j + 1][i - 1];
359: if (j > gys && i + 1 < gxs + gxm) xrb = x[j - 1][i + 1];
361: d1 = (xc - xl);
362: d2 = (xc - xr);
363: d3 = (xc - xt);
364: d4 = (xc - xb);
365: d5 = (xr - xrb);
366: d6 = (xrb - xb);
367: d7 = (xlt - xl);
368: d8 = (xt - xlt);
370: df1dxc = d1 * hydhx;
371: df2dxc = (d1 * hydhx + d4 * hxdhy);
372: df3dxc = d3 * hxdhy;
373: df4dxc = (d2 * hydhx + d3 * hxdhy);
374: df5dxc = d2 * hydhx;
375: df6dxc = d4 * hxdhy;
377: d1 *= rhx;
378: d2 *= rhx;
379: d3 *= rhy;
380: d4 *= rhy;
381: d5 *= rhy;
382: d6 *= rhx;
383: d7 *= rhy;
384: d8 *= rhx;
386: f1 = PetscSqrtReal(1.0 + d1 * d1 + d7 * d7);
387: f2 = PetscSqrtReal(1.0 + d1 * d1 + d4 * d4);
388: f3 = PetscSqrtReal(1.0 + d3 * d3 + d8 * d8);
389: f4 = PetscSqrtReal(1.0 + d3 * d3 + d2 * d2);
390: f5 = PetscSqrtReal(1.0 + d2 * d2 + d5 * d5);
391: f6 = PetscSqrtReal(1.0 + d4 * d4 + d6 * d6);
393: ft = ft + (f2 + f4);
395: df1dxc /= f1;
396: df2dxc /= f2;
397: df3dxc /= f3;
398: df4dxc /= f4;
399: df5dxc /= f5;
400: df6dxc /= f6;
402: g[j][i] = (df1dxc + df2dxc + df3dxc + df4dxc + df5dxc + df6dxc) * 0.5;
403: }
404: }
406: /* Compute triangular areas along the border of the domain. */
407: if (xs == 0) { /* left side */
408: for (j = ys; j < ys + ym; j++) {
409: d3 = (user->left[j - ys + 1] - user->left[j - ys + 2]) * rhy;
410: d2 = (user->left[j - ys + 1] - x[j][0]) * rhx;
411: ft = ft + PetscSqrtReal(1.0 + d3 * d3 + d2 * d2);
412: }
413: }
414: if (ys == 0) { /* bottom side */
415: for (i = xs; i < xs + xm; i++) {
416: d2 = (user->bottom[i + 1 - xs] - user->bottom[i - xs + 2]) * rhx;
417: d3 = (user->bottom[i - xs + 1] - x[0][i]) * rhy;
418: ft = ft + PetscSqrtReal(1.0 + d3 * d3 + d2 * d2);
419: }
420: }
422: if (xs + xm == mx) { /* right side */
423: for (j = ys; j < ys + ym; j++) {
424: d1 = (x[j][mx - 1] - user->right[j - ys + 1]) * rhx;
425: d4 = (user->right[j - ys] - user->right[j - ys + 1]) * rhy;
426: ft = ft + PetscSqrtReal(1.0 + d1 * d1 + d4 * d4);
427: }
428: }
429: if (ys + ym == my) { /* top side */
430: for (i = xs; i < xs + xm; i++) {
431: d1 = (x[my - 1][i] - user->top[i - xs + 1]) * rhy;
432: d4 = (user->top[i - xs + 1] - user->top[i - xs]) * rhx;
433: ft = ft + PetscSqrtReal(1.0 + d1 * d1 + d4 * d4);
434: }
435: }
437: if (ys == 0 && xs == 0) {
438: d1 = (user->left[0] - user->left[1]) / hy;
439: d2 = (user->bottom[0] - user->bottom[1]) * rhx;
440: ft += PetscSqrtReal(1.0 + d1 * d1 + d2 * d2);
441: }
442: if (ys + ym == my && xs + xm == mx) {
443: d1 = (user->right[ym + 1] - user->right[ym]) * rhy;
444: d2 = (user->top[xm + 1] - user->top[xm]) * rhx;
445: ft += PetscSqrtReal(1.0 + d1 * d1 + d2 * d2);
446: }
448: ft = ft * area;
449: PetscCallMPI(MPIU_Allreduce(&ft, fcn, 1, MPIU_REAL, MPIU_SUM, MPI_COMM_WORLD));
451: /* Restore vectors */
452: PetscCall(DMDAVecRestoreArray(user->dm, localX, (void **)&x));
453: PetscCall(DMDAVecRestoreArray(user->dm, G, (void **)&g));
454: PetscCall(DMRestoreLocalVector(user->dm, &localX));
455: PetscCall(PetscLogFlops(67.0 * xm * ym));
456: PetscFunctionReturn(PETSC_SUCCESS);
457: }
459: PetscErrorCode FormGradient(Tao tao, Vec X, Vec G, void *userCtx)
460: {
461: PetscReal fcn;
463: PetscFunctionBegin;
464: PetscCall(FormFunctionGradient(tao, X, &fcn, G, userCtx));
465: PetscFunctionReturn(PETSC_SUCCESS);
466: }
468: /* ------------------------------------------------------------------- */
469: /*
470: FormHessian - Evaluates Hessian matrix.
472: Input Parameters:
473: . tao - the Tao context
474: . x - input vector
475: . ptr - optional user-defined context, as set by TaoSetHessian()
477: Output Parameters:
478: . H - Hessian matrix
479: . Hpre - optionally different matrix used to compute the preconditioner
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, PetscCtx 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*/