Actual source code: plate2.c
1: #include <petscdmda.h>
2: #include <petsctao.h>
4: static char help[] = "This example demonstrates use of the TAO package to \n\
5: solve an unconstrained minimization problem. This example is based on a \n\
6: problem from the MINPACK-2 test suite. Given a rectangular 2-D domain, \n\
7: boundary values along the edges of the domain, and a plate represented by \n\
8: lower boundary conditions, the objective is to find the\n\
9: surface with the minimal area that satisfies the boundary conditions.\n\
10: The command line options are:\n\
11: -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
12: -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
13: -bmx <bxg>, where <bxg> = number of grid points under plate in 1st direction\n\
14: -bmy <byg>, where <byg> = number of grid points under plate in 2nd direction\n\
15: -bheight <ht>, where <ht> = height of the plate\n\
16: -start <st>, where <st> =0 for zero vector, <st> >0 for random start, and <st> <0 \n\
17: for an average of the boundary conditions\n\n";
19: /*
20: User-defined application context - contains data needed by the
21: application-provided call-back routines, FormFunctionGradient(),
22: FormHessian().
23: */
24: typedef struct {
25: /* problem parameters */
26: PetscReal bheight; /* Height of plate under the surface */
27: PetscInt mx, my; /* discretization in x, y directions */
28: PetscInt bmx, bmy; /* Size of plate under the surface */
29: Vec Bottom, Top, Left, Right; /* boundary values */
31: /* Working space */
32: Vec localX, localV; /* ghosted local vector */
33: DM dm; /* distributed array data structure */
34: Mat H;
35: } AppCtx;
37: /* -------- User-defined Routines --------- */
39: static PetscErrorCode MSA_BoundaryConditions(AppCtx *);
40: static PetscErrorCode MSA_InitialPoint(AppCtx *, Vec);
41: static PetscErrorCode MSA_Plate(Vec, Vec, void *);
42: PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);
43: PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void *);
45: /* For testing matrix-free submatrices */
46: PetscErrorCode MatrixFreeHessian(Tao, Vec, Mat, Mat, void *);
47: PetscErrorCode MyMatMult(Mat, Vec, Vec);
49: int main(int argc, char **argv)
50: {
51: PetscInt Nx, Ny; /* number of processors in x- and y- directions */
52: PetscInt m, N; /* number of local and global elements in vectors */
53: Vec x, xl, xu; /* solution vector and bounds*/
54: PetscBool flg; /* A return variable when checking for user options */
55: Tao tao; /* Tao solver context */
56: ISLocalToGlobalMapping isltog; /* local-to-global mapping object */
57: Mat H_shell; /* to test matrix-free submatrices */
58: AppCtx user; /* user-defined work context */
60: /* Initialize PETSc, TAO */
61: PetscFunctionBeginUser;
62: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
64: /* Specify default dimension of the problem */
65: user.mx = 10;
66: user.my = 10;
67: user.bheight = 0.1;
69: /* Check for any command line arguments that override defaults */
70: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &user.mx, &flg));
71: PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &user.my, &flg));
72: PetscCall(PetscOptionsGetReal(NULL, NULL, "-bheight", &user.bheight, &flg));
74: user.bmx = user.mx / 2;
75: user.bmy = user.my / 2;
76: PetscCall(PetscOptionsGetInt(NULL, NULL, "-bmx", &user.bmx, &flg));
77: PetscCall(PetscOptionsGetInt(NULL, NULL, "-bmy", &user.bmy, &flg));
79: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n---- Minimum Surface Area With Plate Problem -----\n"));
80: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "mx:%" PetscInt_FMT ", my:%" PetscInt_FMT ", bmx:%" PetscInt_FMT ", bmy:%" PetscInt_FMT ", height:%g\n", user.mx, user.my, user.bmx, user.bmy, (double)user.bheight));
82: /* Calculate any derived values from parameters */
83: N = user.mx * user.my;
85: /* Let Petsc determine the dimensions of the local vectors */
86: Nx = PETSC_DECIDE;
87: Ny = PETSC_DECIDE;
89: /*
90: A two dimensional distributed array will help define this problem,
91: which derives from an elliptic PDE on two dimensional domain. From
92: the distributed array, Create the vectors.
93: */
94: PetscCall(DMDACreate2d(MPI_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, user.mx, user.my, Nx, Ny, 1, 1, NULL, NULL, &user.dm));
95: PetscCall(DMSetFromOptions(user.dm));
96: PetscCall(DMSetUp(user.dm));
97: /*
98: Extract global and local vectors from DM; The local vectors are
99: used solely as work space for the evaluation of the function,
100: gradient, and Hessian. Duplicate for remaining vectors that are
101: the same types.
102: */
103: PetscCall(DMCreateGlobalVector(user.dm, &x)); /* Solution */
104: PetscCall(DMCreateLocalVector(user.dm, &user.localX));
105: PetscCall(VecDuplicate(user.localX, &user.localV));
107: PetscCall(VecDuplicate(x, &xl));
108: PetscCall(VecDuplicate(x, &xu));
110: /* The TAO code begins here */
112: /*
113: Create TAO solver and set desired solution method
114: The method must either be TAOTRON or TAOBLMVM
115: If TAOBLMVM is used, then hessian function is not called.
116: */
117: PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao));
118: PetscCall(TaoSetType(tao, TAOBLMVM));
120: /* Set initial solution guess; */
121: PetscCall(MSA_BoundaryConditions(&user));
122: PetscCall(MSA_InitialPoint(&user, x));
123: PetscCall(TaoSetSolution(tao, x));
125: /* Set routines for function, gradient and hessian evaluation */
126: PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user));
128: PetscCall(VecGetLocalSize(x, &m));
129: PetscCall(DMCreateMatrix(user.dm, &user.H));
130: PetscCall(MatSetOption(user.H, MAT_SYMMETRIC, PETSC_TRUE));
132: PetscCall(DMGetLocalToGlobalMapping(user.dm, &isltog));
133: PetscCall(MatSetLocalToGlobalMapping(user.H, isltog, isltog));
134: PetscCall(PetscOptionsHasName(NULL, NULL, "-matrixfree", &flg));
135: if (flg) {
136: PetscCall(MatCreateShell(PETSC_COMM_WORLD, m, m, N, N, (void *)&user, &H_shell));
137: PetscCall(MatShellSetOperation(H_shell, MATOP_MULT, (void (*)(void))MyMatMult));
138: PetscCall(MatSetOption(H_shell, MAT_SYMMETRIC, PETSC_TRUE));
139: PetscCall(TaoSetHessian(tao, H_shell, H_shell, MatrixFreeHessian, (void *)&user));
140: } else {
141: PetscCall(TaoSetHessian(tao, user.H, user.H, FormHessian, (void *)&user));
142: }
144: /* Set Variable bounds */
145: PetscCall(MSA_Plate(xl, xu, (void *)&user));
146: PetscCall(TaoSetVariableBounds(tao, xl, xu));
148: /* Check for any tao command line options */
149: PetscCall(TaoSetFromOptions(tao));
151: /* SOLVE THE APPLICATION */
152: PetscCall(TaoSolve(tao));
154: PetscCall(TaoView(tao, PETSC_VIEWER_STDOUT_WORLD));
156: /* Free TAO data structures */
157: PetscCall(TaoDestroy(&tao));
159: /* Free PETSc data structures */
160: PetscCall(VecDestroy(&x));
161: PetscCall(VecDestroy(&xl));
162: PetscCall(VecDestroy(&xu));
163: PetscCall(MatDestroy(&user.H));
164: PetscCall(VecDestroy(&user.localX));
165: PetscCall(VecDestroy(&user.localV));
166: PetscCall(VecDestroy(&user.Bottom));
167: PetscCall(VecDestroy(&user.Top));
168: PetscCall(VecDestroy(&user.Left));
169: PetscCall(VecDestroy(&user.Right));
170: PetscCall(DMDestroy(&user.dm));
171: if (flg) PetscCall(MatDestroy(&H_shell));
172: PetscCall(PetscFinalize());
173: return 0;
174: }
176: /* FormFunctionGradient - Evaluates f(x) and gradient g(x).
178: Input Parameters:
179: . tao - the Tao context
180: . X - input vector
181: . userCtx - optional user-defined context, as set by TaoSetObjectiveAndGradient()
183: Output Parameters:
184: . fcn - the function value
185: . G - vector containing the newly evaluated gradient
187: Notes:
188: In this case, we discretize the domain and Create triangles. The
189: surface of each triangle is planar, whose surface area can be easily
190: computed. The total surface area is found by sweeping through the grid
191: and computing the surface area of the two triangles that have their
192: right angle at the grid point. The diagonal line segments on the
193: grid that define the triangles run from top left to lower right.
194: The numbering of points starts at the lower left and runs left to
195: right, then bottom to top.
196: */
197: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn, Vec G, void *userCtx)
198: {
199: AppCtx *user = (AppCtx *)userCtx;
200: PetscInt i, j, row;
201: PetscInt mx = user->mx, my = user->my;
202: PetscInt xs, xm, gxs, gxm, ys, ym, gys, gym;
203: PetscReal ft = 0;
204: PetscReal zero = 0.0;
205: PetscReal hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), hydhx = hy / hx, hxdhy = hx / hy, area = 0.5 * hx * hy;
206: PetscReal rhx = mx + 1, rhy = my + 1;
207: PetscReal f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb;
208: PetscReal df1dxc, df2dxc, df3dxc, df4dxc, df5dxc, df6dxc;
209: PetscReal *g, *x, *left, *right, *bottom, *top;
210: Vec localX = user->localX, localG = user->localV;
212: PetscFunctionBeginUser;
213: /* Get local mesh boundaries */
214: PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
215: PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL));
217: /* Scatter ghost points to local vector */
218: PetscCall(DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, localX));
219: PetscCall(DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, localX));
221: /* Initialize vector to zero */
222: PetscCall(VecSet(localG, zero));
224: /* Get pointers to vector data */
225: PetscCall(VecGetArray(localX, &x));
226: PetscCall(VecGetArray(localG, &g));
227: PetscCall(VecGetArray(user->Top, &top));
228: PetscCall(VecGetArray(user->Bottom, &bottom));
229: PetscCall(VecGetArray(user->Left, &left));
230: PetscCall(VecGetArray(user->Right, &right));
232: /* Compute function over the locally owned part of the mesh */
233: for (j = ys; j < ys + ym; j++) {
234: for (i = xs; i < xs + xm; i++) {
235: row = (j - gys) * gxm + (i - gxs);
237: xc = x[row];
238: xlt = xrb = xl = xr = xb = xt = xc;
240: if (i == 0) { /* left side */
241: xl = left[j - ys + 1];
242: xlt = left[j - ys + 2];
243: } else {
244: xl = x[row - 1];
245: }
247: if (j == 0) { /* bottom side */
248: xb = bottom[i - xs + 1];
249: xrb = bottom[i - xs + 2];
250: } else {
251: xb = x[row - gxm];
252: }
254: if (i + 1 == gxs + gxm) { /* right side */
255: xr = right[j - ys + 1];
256: xrb = right[j - ys];
257: } else {
258: xr = x[row + 1];
259: }
261: if (j + 1 == gys + gym) { /* top side */
262: xt = top[i - xs + 1];
263: xlt = top[i - xs];
264: } else {
265: xt = x[row + gxm];
266: }
268: if (i > gxs && j + 1 < gys + gym) xlt = x[row - 1 + gxm];
269: if (j > gys && i + 1 < gxs + gxm) xrb = x[row + 1 - gxm];
271: d1 = (xc - xl);
272: d2 = (xc - xr);
273: d3 = (xc - xt);
274: d4 = (xc - xb);
275: d5 = (xr - xrb);
276: d6 = (xrb - xb);
277: d7 = (xlt - xl);
278: d8 = (xt - xlt);
280: df1dxc = d1 * hydhx;
281: df2dxc = (d1 * hydhx + d4 * hxdhy);
282: df3dxc = d3 * hxdhy;
283: df4dxc = (d2 * hydhx + d3 * hxdhy);
284: df5dxc = d2 * hydhx;
285: df6dxc = d4 * hxdhy;
287: d1 *= rhx;
288: d2 *= rhx;
289: d3 *= rhy;
290: d4 *= rhy;
291: d5 *= rhy;
292: d6 *= rhx;
293: d7 *= rhy;
294: d8 *= rhx;
296: f1 = PetscSqrtScalar(1.0 + d1 * d1 + d7 * d7);
297: f2 = PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
298: f3 = PetscSqrtScalar(1.0 + d3 * d3 + d8 * d8);
299: f4 = PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
300: f5 = PetscSqrtScalar(1.0 + d2 * d2 + d5 * d5);
301: f6 = PetscSqrtScalar(1.0 + d4 * d4 + d6 * d6);
303: ft = ft + (f2 + f4);
305: df1dxc /= f1;
306: df2dxc /= f2;
307: df3dxc /= f3;
308: df4dxc /= f4;
309: df5dxc /= f5;
310: df6dxc /= f6;
312: g[row] = (df1dxc + df2dxc + df3dxc + df4dxc + df5dxc + df6dxc) * 0.5;
313: }
314: }
316: /* Compute triangular areas along the border of the domain. */
317: if (xs == 0) { /* left side */
318: for (j = ys; j < ys + ym; j++) {
319: d3 = (left[j - ys + 1] - left[j - ys + 2]) * rhy;
320: d2 = (left[j - ys + 1] - x[(j - gys) * gxm]) * rhx;
321: ft = ft + PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
322: }
323: }
324: if (ys == 0) { /* bottom side */
325: for (i = xs; i < xs + xm; i++) {
326: d2 = (bottom[i + 1 - xs] - bottom[i - xs + 2]) * rhx;
327: d3 = (bottom[i - xs + 1] - x[i - gxs]) * rhy;
328: ft = ft + PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
329: }
330: }
332: if (xs + xm == mx) { /* right side */
333: for (j = ys; j < ys + ym; j++) {
334: d1 = (x[(j + 1 - gys) * gxm - 1] - right[j - ys + 1]) * rhx;
335: d4 = (right[j - ys] - right[j - ys + 1]) * rhy;
336: ft = ft + PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
337: }
338: }
339: if (ys + ym == my) { /* top side */
340: for (i = xs; i < xs + xm; i++) {
341: d1 = (x[(gym - 1) * gxm + i - gxs] - top[i - xs + 1]) * rhy;
342: d4 = (top[i - xs + 1] - top[i - xs]) * rhx;
343: ft = ft + PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
344: }
345: }
347: if (ys == 0 && xs == 0) {
348: d1 = (left[0] - left[1]) * rhy;
349: d2 = (bottom[0] - bottom[1]) * rhx;
350: ft += PetscSqrtScalar(1.0 + d1 * d1 + d2 * d2);
351: }
352: if (ys + ym == my && xs + xm == mx) {
353: d1 = (right[ym + 1] - right[ym]) * rhy;
354: d2 = (top[xm + 1] - top[xm]) * rhx;
355: ft += PetscSqrtScalar(1.0 + d1 * d1 + d2 * d2);
356: }
358: ft = ft * area;
359: PetscCall(MPIU_Allreduce(&ft, fcn, 1, MPIU_REAL, MPIU_SUM, MPI_COMM_WORLD));
361: /* Restore vectors */
362: PetscCall(VecRestoreArray(localX, &x));
363: PetscCall(VecRestoreArray(localG, &g));
364: PetscCall(VecRestoreArray(user->Left, &left));
365: PetscCall(VecRestoreArray(user->Top, &top));
366: PetscCall(VecRestoreArray(user->Bottom, &bottom));
367: PetscCall(VecRestoreArray(user->Right, &right));
369: /* Scatter values to global vector */
370: PetscCall(DMLocalToGlobalBegin(user->dm, localG, INSERT_VALUES, G));
371: PetscCall(DMLocalToGlobalEnd(user->dm, localG, INSERT_VALUES, G));
373: PetscCall(PetscLogFlops(70.0 * xm * ym));
375: PetscFunctionReturn(PETSC_SUCCESS);
376: }
378: /* ------------------------------------------------------------------- */
379: /*
380: FormHessian - Evaluates Hessian matrix.
382: Input Parameters:
383: . tao - the Tao context
384: . x - input vector
385: . ptr - optional user-defined context, as set by TaoSetHessian()
387: Output Parameters:
388: . A - Hessian matrix
389: . B - optionally different preconditioning matrix
391: Notes:
392: Due to mesh point reordering with DMs, we must always work
393: with the local mesh points, and then transform them to the new
394: global numbering with the local-to-global mapping. We cannot work
395: directly with the global numbers for the original uniprocessor mesh!
397: Two methods are available for imposing this transformation
398: when setting matrix entries:
399: (A) MatSetValuesLocal(), using the local ordering (including
400: ghost points!)
401: - Do the following two steps once, before calling TaoSolve()
402: - Use DMGetISLocalToGlobalMapping() to extract the
403: local-to-global map from the DM
404: - Associate this map with the matrix by calling
405: MatSetLocalToGlobalMapping()
406: - Then set matrix entries using the local ordering
407: by calling MatSetValuesLocal()
408: (B) MatSetValues(), using the global ordering
409: - Use DMGetGlobalIndices() to extract the local-to-global map
410: - Then apply this map explicitly yourself
411: - Set matrix entries using the global ordering by calling
412: MatSetValues()
413: Option (A) seems cleaner/easier in many cases, and is the procedure
414: used in this example.
415: */
416: PetscErrorCode FormHessian(Tao tao, Vec X, Mat Hptr, Mat Hessian, void *ptr)
417: {
418: AppCtx *user = (AppCtx *)ptr;
419: PetscInt i, j, k, row;
420: PetscInt mx = user->mx, my = user->my;
421: PetscInt xs, xm, gxs, gxm, ys, ym, gys, gym, col[7];
422: PetscReal hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), hydhx = hy / hx, hxdhy = hx / hy;
423: PetscReal rhx = mx + 1, rhy = my + 1;
424: PetscReal f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb;
425: PetscReal hl, hr, ht, hb, hc, htl, hbr;
426: PetscReal *x, *left, *right, *bottom, *top;
427: PetscReal v[7];
428: Vec localX = user->localX;
429: PetscBool assembled;
431: PetscFunctionBeginUser;
432: /* Set various matrix options */
433: PetscCall(MatSetOption(Hessian, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE));
435: /* Initialize matrix entries to zero */
436: PetscCall(MatAssembled(Hessian, &assembled));
437: if (assembled) PetscCall(MatZeroEntries(Hessian));
439: /* Get local mesh boundaries */
440: PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
441: PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL));
443: /* Scatter ghost points to local vector */
444: PetscCall(DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, localX));
445: PetscCall(DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, localX));
447: /* Get pointers to vector data */
448: PetscCall(VecGetArray(localX, &x));
449: PetscCall(VecGetArray(user->Top, &top));
450: PetscCall(VecGetArray(user->Bottom, &bottom));
451: PetscCall(VecGetArray(user->Left, &left));
452: PetscCall(VecGetArray(user->Right, &right));
454: /* Compute Hessian over the locally owned part of the mesh */
456: for (i = xs; i < xs + xm; i++) {
457: for (j = ys; j < ys + ym; j++) {
458: row = (j - gys) * gxm + (i - gxs);
460: xc = x[row];
461: xlt = xrb = xl = xr = xb = xt = xc;
463: /* Left side */
464: if (i == gxs) {
465: xl = left[j - ys + 1];
466: xlt = left[j - ys + 2];
467: } else {
468: xl = x[row - 1];
469: }
471: if (j == gys) {
472: xb = bottom[i - xs + 1];
473: xrb = bottom[i - xs + 2];
474: } else {
475: xb = x[row - gxm];
476: }
478: if (i + 1 == gxs + gxm) {
479: xr = right[j - ys + 1];
480: xrb = right[j - ys];
481: } else {
482: xr = x[row + 1];
483: }
485: if (j + 1 == gys + gym) {
486: xt = top[i - xs + 1];
487: xlt = top[i - xs];
488: } else {
489: xt = x[row + gxm];
490: }
492: if (i > gxs && j + 1 < gys + gym) xlt = x[row - 1 + gxm];
493: if (j > gys && i + 1 < gxs + gxm) xrb = x[row + 1 - gxm];
495: d1 = (xc - xl) * rhx;
496: d2 = (xc - xr) * rhx;
497: d3 = (xc - xt) * rhy;
498: d4 = (xc - xb) * rhy;
499: d5 = (xrb - xr) * rhy;
500: d6 = (xrb - xb) * rhx;
501: d7 = (xlt - xl) * rhy;
502: d8 = (xlt - xt) * rhx;
504: f1 = PetscSqrtScalar(1.0 + d1 * d1 + d7 * d7);
505: f2 = PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
506: f3 = PetscSqrtScalar(1.0 + d3 * d3 + d8 * d8);
507: f4 = PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
508: f5 = PetscSqrtScalar(1.0 + d2 * d2 + d5 * d5);
509: f6 = PetscSqrtScalar(1.0 + d4 * d4 + d6 * d6);
511: hl = (-hydhx * (1.0 + d7 * d7) + d1 * d7) / (f1 * f1 * f1) + (-hydhx * (1.0 + d4 * d4) + d1 * d4) / (f2 * f2 * f2);
512: hr = (-hydhx * (1.0 + d5 * d5) + d2 * d5) / (f5 * f5 * f5) + (-hydhx * (1.0 + d3 * d3) + d2 * d3) / (f4 * f4 * f4);
513: ht = (-hxdhy * (1.0 + d8 * d8) + d3 * d8) / (f3 * f3 * f3) + (-hxdhy * (1.0 + d2 * d2) + d2 * d3) / (f4 * f4 * f4);
514: hb = (-hxdhy * (1.0 + d6 * d6) + d4 * d6) / (f6 * f6 * f6) + (-hxdhy * (1.0 + d1 * d1) + d1 * d4) / (f2 * f2 * f2);
516: hbr = -d2 * d5 / (f5 * f5 * f5) - d4 * d6 / (f6 * f6 * f6);
517: htl = -d1 * d7 / (f1 * f1 * f1) - d3 * d8 / (f3 * f3 * f3);
519: 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);
521: hl *= 0.5;
522: hr *= 0.5;
523: ht *= 0.5;
524: hb *= 0.5;
525: hbr *= 0.5;
526: htl *= 0.5;
527: hc *= 0.5;
529: k = 0;
530: if (j > 0) {
531: v[k] = hb;
532: col[k] = row - gxm;
533: k++;
534: }
536: if (j > 0 && i < mx - 1) {
537: v[k] = hbr;
538: col[k] = row - gxm + 1;
539: k++;
540: }
542: if (i > 0) {
543: v[k] = hl;
544: col[k] = row - 1;
545: k++;
546: }
548: v[k] = hc;
549: col[k] = row;
550: k++;
552: if (i < mx - 1) {
553: v[k] = hr;
554: col[k] = row + 1;
555: k++;
556: }
558: if (i > 0 && j < my - 1) {
559: v[k] = htl;
560: col[k] = row + gxm - 1;
561: k++;
562: }
564: if (j < my - 1) {
565: v[k] = ht;
566: col[k] = row + gxm;
567: k++;
568: }
570: /*
571: Set matrix values using local numbering, which was defined
572: earlier, in the main routine.
573: */
574: PetscCall(MatSetValuesLocal(Hessian, 1, &row, k, col, v, INSERT_VALUES));
575: }
576: }
578: /* Restore vectors */
579: PetscCall(VecRestoreArray(localX, &x));
580: PetscCall(VecRestoreArray(user->Left, &left));
581: PetscCall(VecRestoreArray(user->Top, &top));
582: PetscCall(VecRestoreArray(user->Bottom, &bottom));
583: PetscCall(VecRestoreArray(user->Right, &right));
585: /* Assemble the matrix */
586: PetscCall(MatAssemblyBegin(Hessian, MAT_FINAL_ASSEMBLY));
587: PetscCall(MatAssemblyEnd(Hessian, MAT_FINAL_ASSEMBLY));
589: PetscCall(PetscLogFlops(199.0 * xm * ym));
590: PetscFunctionReturn(PETSC_SUCCESS);
591: }
593: /* ------------------------------------------------------------------- */
594: /*
595: MSA_BoundaryConditions - Calculates the boundary conditions for
596: the region.
598: Input Parameter:
599: . user - user-defined application context
601: Output Parameter:
602: . user - user-defined application context
603: */
604: static PetscErrorCode MSA_BoundaryConditions(AppCtx *user)
605: {
606: PetscInt i, j, k, maxits = 5, limit = 0;
607: PetscInt xs, ys, xm, ym, gxs, gys, gxm, gym;
608: PetscInt mx = user->mx, my = user->my;
609: PetscInt bsize = 0, lsize = 0, tsize = 0, rsize = 0;
610: PetscReal one = 1.0, two = 2.0, three = 3.0, scl = 1.0, tol = 1e-10;
611: PetscReal fnorm, det, hx, hy, xt = 0, yt = 0;
612: PetscReal u1, u2, nf1, nf2, njac11, njac12, njac21, njac22;
613: PetscReal b = -0.5, t = 0.5, l = -0.5, r = 0.5;
614: PetscReal *boundary;
615: PetscBool flg;
616: Vec Bottom, Top, Right, Left;
618: PetscFunctionBeginUser;
619: /* Get local mesh boundaries */
620: PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
621: PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL));
623: bsize = xm + 2;
624: lsize = ym + 2;
625: rsize = ym + 2;
626: tsize = xm + 2;
628: PetscCall(VecCreateMPI(MPI_COMM_WORLD, bsize, PETSC_DECIDE, &Bottom));
629: PetscCall(VecCreateMPI(MPI_COMM_WORLD, tsize, PETSC_DECIDE, &Top));
630: PetscCall(VecCreateMPI(MPI_COMM_WORLD, lsize, PETSC_DECIDE, &Left));
631: PetscCall(VecCreateMPI(MPI_COMM_WORLD, rsize, PETSC_DECIDE, &Right));
633: user->Top = Top;
634: user->Left = Left;
635: user->Bottom = Bottom;
636: user->Right = Right;
638: hx = (r - l) / (mx + 1);
639: hy = (t - b) / (my + 1);
641: for (j = 0; j < 4; j++) {
642: if (j == 0) {
643: yt = b;
644: xt = l + hx * xs;
645: limit = bsize;
646: PetscCall(VecGetArray(Bottom, &boundary));
647: } else if (j == 1) {
648: yt = t;
649: xt = l + hx * xs;
650: limit = tsize;
651: PetscCall(VecGetArray(Top, &boundary));
652: } else if (j == 2) {
653: yt = b + hy * ys;
654: xt = l;
655: limit = lsize;
656: PetscCall(VecGetArray(Left, &boundary));
657: } else if (j == 3) {
658: yt = b + hy * ys;
659: xt = r;
660: limit = rsize;
661: PetscCall(VecGetArray(Right, &boundary));
662: }
664: for (i = 0; i < limit; i++) {
665: u1 = xt;
666: u2 = -yt;
667: for (k = 0; k < maxits; k++) {
668: nf1 = u1 + u1 * u2 * u2 - u1 * u1 * u1 / three - xt;
669: nf2 = -u2 - u1 * u1 * u2 + u2 * u2 * u2 / three - yt;
670: fnorm = PetscSqrtScalar(nf1 * nf1 + nf2 * nf2);
671: if (fnorm <= tol) break;
672: njac11 = one + u2 * u2 - u1 * u1;
673: njac12 = two * u1 * u2;
674: njac21 = -two * u1 * u2;
675: njac22 = -one - u1 * u1 + u2 * u2;
676: det = njac11 * njac22 - njac21 * njac12;
677: u1 = u1 - (njac22 * nf1 - njac12 * nf2) / det;
678: u2 = u2 - (njac11 * nf2 - njac21 * nf1) / det;
679: }
681: boundary[i] = u1 * u1 - u2 * u2;
682: if (j == 0 || j == 1) {
683: xt = xt + hx;
684: } else if (j == 2 || j == 3) {
685: yt = yt + hy;
686: }
687: }
688: if (j == 0) {
689: PetscCall(VecRestoreArray(Bottom, &boundary));
690: } else if (j == 1) {
691: PetscCall(VecRestoreArray(Top, &boundary));
692: } else if (j == 2) {
693: PetscCall(VecRestoreArray(Left, &boundary));
694: } else if (j == 3) {
695: PetscCall(VecRestoreArray(Right, &boundary));
696: }
697: }
699: /* Scale the boundary if desired */
701: PetscCall(PetscOptionsGetReal(NULL, NULL, "-bottom", &scl, &flg));
702: if (flg) PetscCall(VecScale(Bottom, scl));
703: PetscCall(PetscOptionsGetReal(NULL, NULL, "-top", &scl, &flg));
704: if (flg) PetscCall(VecScale(Top, scl));
705: PetscCall(PetscOptionsGetReal(NULL, NULL, "-right", &scl, &flg));
706: if (flg) PetscCall(VecScale(Right, scl));
708: PetscCall(PetscOptionsGetReal(NULL, NULL, "-left", &scl, &flg));
709: if (flg) PetscCall(VecScale(Left, scl));
710: PetscFunctionReturn(PETSC_SUCCESS);
711: }
713: /* ------------------------------------------------------------------- */
714: /*
715: MSA_Plate - Calculates an obstacle for surface to stretch over.
717: Input Parameter:
718: . user - user-defined application context
720: Output Parameter:
721: . user - user-defined application context
722: */
723: static PetscErrorCode MSA_Plate(Vec XL, Vec XU, void *ctx)
724: {
725: AppCtx *user = (AppCtx *)ctx;
726: PetscInt i, j, row;
727: PetscInt xs, ys, xm, ym;
728: PetscInt mx = user->mx, my = user->my, bmy, bmx;
729: PetscReal t1, t2, t3;
730: PetscReal *xl, lb = PETSC_NINFINITY, ub = PETSC_INFINITY;
731: PetscBool cylinder;
733: PetscFunctionBeginUser;
734: user->bmy = PetscMax(0, user->bmy);
735: user->bmy = PetscMin(my, user->bmy);
736: user->bmx = PetscMax(0, user->bmx);
737: user->bmx = PetscMin(mx, user->bmx);
738: bmy = user->bmy;
739: bmx = user->bmx;
741: PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
743: PetscCall(VecSet(XL, lb));
744: PetscCall(VecSet(XU, ub));
746: PetscCall(VecGetArray(XL, &xl));
748: PetscCall(PetscOptionsHasName(NULL, NULL, "-cylinder", &cylinder));
749: /* Compute the optional lower box */
750: if (cylinder) {
751: for (i = xs; i < xs + xm; i++) {
752: for (j = ys; j < ys + ym; j++) {
753: row = (j - ys) * xm + (i - xs);
754: t1 = (2.0 * i - mx) * bmy;
755: t2 = (2.0 * j - my) * bmx;
756: t3 = bmx * bmx * bmy * bmy;
757: if (t1 * t1 + t2 * t2 <= t3) xl[row] = user->bheight;
758: }
759: }
760: } else {
761: /* Compute the optional lower box */
762: for (i = xs; i < xs + xm; i++) {
763: for (j = ys; j < ys + ym; j++) {
764: row = (j - ys) * xm + (i - xs);
765: if (i >= (mx - bmx) / 2 && i < mx - (mx - bmx) / 2 && j >= (my - bmy) / 2 && j < my - (my - bmy) / 2) xl[row] = user->bheight;
766: }
767: }
768: }
769: PetscCall(VecRestoreArray(XL, &xl));
771: PetscFunctionReturn(PETSC_SUCCESS);
772: }
774: /* ------------------------------------------------------------------- */
775: /*
776: MSA_InitialPoint - Calculates the initial guess in one of three ways.
778: Input Parameters:
779: . user - user-defined application context
780: . X - vector for initial guess
782: Output Parameters:
783: . X - newly computed initial guess
784: */
785: static PetscErrorCode MSA_InitialPoint(AppCtx *user, Vec X)
786: {
787: PetscInt start = -1, i, j;
788: PetscReal zero = 0.0;
789: PetscBool flg;
791: PetscFunctionBeginUser;
792: PetscCall(PetscOptionsGetInt(NULL, NULL, "-start", &start, &flg));
793: if (flg && start == 0) { /* The zero vector is reasonable */
794: PetscCall(VecSet(X, zero));
795: } else if (flg && start > 0) { /* Try a random start between -0.5 and 0.5 */
796: PetscRandom rctx;
797: PetscReal np5 = -0.5;
799: PetscCall(PetscRandomCreate(MPI_COMM_WORLD, &rctx));
800: for (i = 0; i < start; i++) PetscCall(VecSetRandom(X, rctx));
801: PetscCall(PetscRandomDestroy(&rctx));
802: PetscCall(VecShift(X, np5));
804: } else { /* Take an average of the boundary conditions */
806: PetscInt row, xs, xm, gxs, gxm, ys, ym, gys, gym;
807: PetscInt mx = user->mx, my = user->my;
808: PetscReal *x, *left, *right, *bottom, *top;
809: Vec localX = user->localX;
811: /* Get local mesh boundaries */
812: PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
813: PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL));
815: /* Get pointers to vector data */
816: PetscCall(VecGetArray(user->Top, &top));
817: PetscCall(VecGetArray(user->Bottom, &bottom));
818: PetscCall(VecGetArray(user->Left, &left));
819: PetscCall(VecGetArray(user->Right, &right));
821: PetscCall(VecGetArray(localX, &x));
822: /* Perform local computations */
823: for (j = ys; j < ys + ym; j++) {
824: for (i = xs; i < xs + xm; i++) {
825: row = (j - gys) * gxm + (i - gxs);
826: x[row] = ((j + 1) * bottom[i - xs + 1] / my + (my - j + 1) * top[i - xs + 1] / (my + 2) + (i + 1) * left[j - ys + 1] / mx + (mx - i + 1) * right[j - ys + 1] / (mx + 2)) / 2.0;
827: }
828: }
830: /* Restore vectors */
831: PetscCall(VecRestoreArray(localX, &x));
833: PetscCall(VecRestoreArray(user->Left, &left));
834: PetscCall(VecRestoreArray(user->Top, &top));
835: PetscCall(VecRestoreArray(user->Bottom, &bottom));
836: PetscCall(VecRestoreArray(user->Right, &right));
838: /* Scatter values into global vector */
839: PetscCall(DMLocalToGlobalBegin(user->dm, localX, INSERT_VALUES, X));
840: PetscCall(DMLocalToGlobalEnd(user->dm, localX, INSERT_VALUES, X));
841: }
842: PetscFunctionReturn(PETSC_SUCCESS);
843: }
845: /* For testing matrix-free submatrices */
846: PetscErrorCode MatrixFreeHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ptr)
847: {
848: AppCtx *user = (AppCtx *)ptr;
850: PetscFunctionBegin;
851: PetscCall(FormHessian(tao, x, user->H, user->H, ptr));
852: PetscFunctionReturn(PETSC_SUCCESS);
853: }
855: PetscErrorCode MyMatMult(Mat H_shell, Vec X, Vec Y)
856: {
857: void *ptr;
858: AppCtx *user;
860: PetscFunctionBegin;
861: PetscCall(MatShellGetContext(H_shell, &ptr));
862: user = (AppCtx *)ptr;
863: PetscCall(MatMult(user->H, X, Y));
864: PetscFunctionReturn(PETSC_SUCCESS);
865: }
867: /*TEST
869: build:
870: requires: !complex
872: test:
873: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type tron -tao_gatol 1.e-5
874: requires: !single
876: test:
877: suffix: 2
878: nsize: 2
879: args: -tao_smonitor -mx 8 -my 8 -bmx 2 -bmy 5 -bheight 0.3 -tao_type blmvm -tao_gatol 1.e-4
880: requires: !single
882: test:
883: suffix: 3
884: nsize: 3
885: args: -tao_smonitor -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_type tron -tao_gatol 1.e-5
886: requires: !single
888: test:
889: suffix: 4
890: nsize: 3
891: args: -tao_smonitor -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_subset_type mask -tao_type tron -tao_gatol 1.e-5
892: requires: !single
894: test:
895: suffix: 5
896: nsize: 3
897: args: -tao_smonitor -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_subset_type matrixfree -matrixfree -pc_type none -tao_type tron -tao_gatol 1.e-5
898: requires: !single
900: test:
901: suffix: 6
902: nsize: 3
903: args: -tao_smonitor -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_subset_type matrixfree -matrixfree -tao_type blmvm -tao_gatol 1.e-4
904: requires: !single
906: test:
907: suffix: 7
908: nsize: 3
909: args: -tao_smonitor -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_subset_type matrixfree -pc_type none -tao_type gpcg -tao_gatol 1.e-5
910: requires: !single
912: test:
913: suffix: 8
914: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bncg -tao_bncg_type gd -tao_gatol 1e-4
915: requires: !single
917: test:
918: suffix: 9
919: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bncg -tao_gatol 1e-4
920: requires: !single
922: test:
923: suffix: 10
924: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bnls -tao_gatol 1e-5
925: requires: !single
927: test:
928: suffix: 11
929: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntr -tao_gatol 1e-5
930: requires: !single
932: test:
933: suffix: 12
934: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntl -tao_gatol 1e-5
935: requires: !single
937: test:
938: suffix: 13
939: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bnls -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
940: requires: !single
942: test:
943: suffix: 14
944: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntr -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
945: requires: !single
947: test:
948: suffix: 15
949: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntl -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
950: requires: !single
952: test:
953: suffix: 16
954: args: -tao_smonitor -mx 8 -my 8 -bmx 2 -bmy 5 -bheight 0.3 -tao_gatol 1e-4 -tao_type bqnls
955: requires: !single
957: test:
958: suffix: 17
959: args: -tao_smonitor -mx 8 -my 8 -bmx 2 -bmy 5 -bheight 0.3 -tao_gatol 1e-4 -tao_type bqnkls -tao_bqnk_mat_type lmvmbfgs
960: requires: !single
962: test:
963: suffix: 18
964: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bnls -tao_gatol 1e-5 -tao_mf_hessian
965: requires: !single
967: test:
968: suffix: 19
969: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntr -tao_gatol 1e-5 -tao_mf_hessian
970: requires: !single
972: test:
973: suffix: 20
974: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntl -tao_gatol 1e-5 -tao_mf_hessian
976: TEST*/