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 */
62: 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: PetscOptionsGetInt(NULL, NULL, "-mx", &user.mx, &flg);
71: PetscOptionsGetInt(NULL, NULL, "-my", &user.my, &flg);
72: PetscOptionsGetReal(NULL, NULL, "-bheight", &user.bheight, &flg);
74: user.bmx = user.mx / 2;
75: user.bmy = user.my / 2;
76: PetscOptionsGetInt(NULL, NULL, "-bmx", &user.bmx, &flg);
77: PetscOptionsGetInt(NULL, NULL, "-bmy", &user.bmy, &flg);
79: PetscPrintf(PETSC_COMM_WORLD, "\n---- Minimum Surface Area With Plate Problem -----\n");
80: 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: 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: DMSetFromOptions(user.dm);
96: 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: DMCreateGlobalVector(user.dm, &x); /* Solution */
104: DMCreateLocalVector(user.dm, &user.localX);
105: VecDuplicate(user.localX, &user.localV);
107: VecDuplicate(x, &xl);
108: 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: TaoCreate(PETSC_COMM_WORLD, &tao);
118: TaoSetType(tao, TAOBLMVM);
120: /* Set initial solution guess; */
121: MSA_BoundaryConditions(&user);
122: MSA_InitialPoint(&user, x);
123: TaoSetSolution(tao, x);
125: /* Set routines for function, gradient and hessian evaluation */
126: TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user);
128: VecGetLocalSize(x, &m);
129: MatCreateAIJ(MPI_COMM_WORLD, m, m, N, N, 7, NULL, 3, NULL, &(user.H));
130: MatSetOption(user.H, MAT_SYMMETRIC, PETSC_TRUE);
132: DMGetLocalToGlobalMapping(user.dm, &isltog);
133: MatSetLocalToGlobalMapping(user.H, isltog, isltog);
134: PetscOptionsHasName(NULL, NULL, "-matrixfree", &flg);
135: if (flg) {
136: MatCreateShell(PETSC_COMM_WORLD, m, m, N, N, (void *)&user, &H_shell);
137: MatShellSetOperation(H_shell, MATOP_MULT, (void (*)(void))MyMatMult);
138: MatSetOption(H_shell, MAT_SYMMETRIC, PETSC_TRUE);
139: TaoSetHessian(tao, H_shell, H_shell, MatrixFreeHessian, (void *)&user);
140: } else {
141: TaoSetHessian(tao, user.H, user.H, FormHessian, (void *)&user);
142: }
144: /* Set Variable bounds */
145: MSA_Plate(xl, xu, (void *)&user);
146: TaoSetVariableBounds(tao, xl, xu);
148: /* Check for any tao command line options */
149: TaoSetFromOptions(tao);
151: /* SOLVE THE APPLICATION */
152: TaoSolve(tao);
154: TaoView(tao, PETSC_VIEWER_STDOUT_WORLD);
156: /* Free TAO data structures */
157: TaoDestroy(&tao);
159: /* Free PETSc data structures */
160: VecDestroy(&x);
161: VecDestroy(&xl);
162: VecDestroy(&xu);
163: MatDestroy(&user.H);
164: VecDestroy(&user.localX);
165: VecDestroy(&user.localV);
166: VecDestroy(&user.Bottom);
167: VecDestroy(&user.Top);
168: VecDestroy(&user.Left);
169: VecDestroy(&user.Right);
170: DMDestroy(&user.dm);
171: if (flg) MatDestroy(&H_shell);
172: 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: /* Get local mesh boundaries */
213: DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL);
214: DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL);
216: /* Scatter ghost points to local vector */
217: DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, localX);
218: DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, localX);
220: /* Initialize vector to zero */
221: VecSet(localG, zero);
223: /* Get pointers to vector data */
224: VecGetArray(localX, &x);
225: VecGetArray(localG, &g);
226: VecGetArray(user->Top, &top);
227: VecGetArray(user->Bottom, &bottom);
228: VecGetArray(user->Left, &left);
229: VecGetArray(user->Right, &right);
231: /* Compute function over the locally owned part of the mesh */
232: for (j = ys; j < ys + ym; j++) {
233: for (i = xs; i < xs + xm; i++) {
234: row = (j - gys) * gxm + (i - gxs);
236: xc = x[row];
237: xlt = xrb = xl = xr = xb = xt = xc;
239: if (i == 0) { /* left side */
240: xl = left[j - ys + 1];
241: xlt = left[j - ys + 2];
242: } else {
243: xl = x[row - 1];
244: }
246: if (j == 0) { /* bottom side */
247: xb = bottom[i - xs + 1];
248: xrb = bottom[i - xs + 2];
249: } else {
250: xb = x[row - gxm];
251: }
253: if (i + 1 == gxs + gxm) { /* right side */
254: xr = right[j - ys + 1];
255: xrb = right[j - ys];
256: } else {
257: xr = x[row + 1];
258: }
260: if (j + 1 == gys + gym) { /* top side */
261: xt = top[i - xs + 1];
262: xlt = top[i - xs];
263: } else {
264: xt = x[row + gxm];
265: }
267: if (i > gxs && j + 1 < gys + gym) xlt = x[row - 1 + gxm];
268: if (j > gys && i + 1 < gxs + gxm) xrb = x[row + 1 - gxm];
270: d1 = (xc - xl);
271: d2 = (xc - xr);
272: d3 = (xc - xt);
273: d4 = (xc - xb);
274: d5 = (xr - xrb);
275: d6 = (xrb - xb);
276: d7 = (xlt - xl);
277: d8 = (xt - xlt);
279: df1dxc = d1 * hydhx;
280: df2dxc = (d1 * hydhx + d4 * hxdhy);
281: df3dxc = d3 * hxdhy;
282: df4dxc = (d2 * hydhx + d3 * hxdhy);
283: df5dxc = d2 * hydhx;
284: df6dxc = d4 * hxdhy;
286: d1 *= rhx;
287: d2 *= rhx;
288: d3 *= rhy;
289: d4 *= rhy;
290: d5 *= rhy;
291: d6 *= rhx;
292: d7 *= rhy;
293: d8 *= rhx;
295: f1 = PetscSqrtScalar(1.0 + d1 * d1 + d7 * d7);
296: f2 = PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
297: f3 = PetscSqrtScalar(1.0 + d3 * d3 + d8 * d8);
298: f4 = PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
299: f5 = PetscSqrtScalar(1.0 + d2 * d2 + d5 * d5);
300: f6 = PetscSqrtScalar(1.0 + d4 * d4 + d6 * d6);
302: ft = ft + (f2 + f4);
304: df1dxc /= f1;
305: df2dxc /= f2;
306: df3dxc /= f3;
307: df4dxc /= f4;
308: df5dxc /= f5;
309: df6dxc /= f6;
311: g[row] = (df1dxc + df2dxc + df3dxc + df4dxc + df5dxc + df6dxc) * 0.5;
312: }
313: }
315: /* Compute triangular areas along the border of the domain. */
316: if (xs == 0) { /* left side */
317: for (j = ys; j < ys + ym; j++) {
318: d3 = (left[j - ys + 1] - left[j - ys + 2]) * rhy;
319: d2 = (left[j - ys + 1] - x[(j - gys) * gxm]) * rhx;
320: ft = ft + PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
321: }
322: }
323: if (ys == 0) { /* bottom side */
324: for (i = xs; i < xs + xm; i++) {
325: d2 = (bottom[i + 1 - xs] - bottom[i - xs + 2]) * rhx;
326: d3 = (bottom[i - xs + 1] - x[i - gxs]) * rhy;
327: ft = ft + PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
328: }
329: }
331: if (xs + xm == mx) { /* right side */
332: for (j = ys; j < ys + ym; j++) {
333: d1 = (x[(j + 1 - gys) * gxm - 1] - right[j - ys + 1]) * rhx;
334: d4 = (right[j - ys] - right[j - ys + 1]) * rhy;
335: ft = ft + PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
336: }
337: }
338: if (ys + ym == my) { /* top side */
339: for (i = xs; i < xs + xm; i++) {
340: d1 = (x[(gym - 1) * gxm + i - gxs] - top[i - xs + 1]) * rhy;
341: d4 = (top[i - xs + 1] - top[i - xs]) * rhx;
342: ft = ft + PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
343: }
344: }
346: if (ys == 0 && xs == 0) {
347: d1 = (left[0] - left[1]) * rhy;
348: d2 = (bottom[0] - bottom[1]) * rhx;
349: ft += PetscSqrtScalar(1.0 + d1 * d1 + d2 * d2);
350: }
351: if (ys + ym == my && xs + xm == mx) {
352: d1 = (right[ym + 1] - right[ym]) * rhy;
353: d2 = (top[xm + 1] - top[xm]) * rhx;
354: ft += PetscSqrtScalar(1.0 + d1 * d1 + d2 * d2);
355: }
357: ft = ft * area;
358: MPI_Allreduce(&ft, fcn, 1, MPIU_REAL, MPIU_SUM, MPI_COMM_WORLD);
360: /* Restore vectors */
361: VecRestoreArray(localX, &x);
362: VecRestoreArray(localG, &g);
363: VecRestoreArray(user->Left, &left);
364: VecRestoreArray(user->Top, &top);
365: VecRestoreArray(user->Bottom, &bottom);
366: VecRestoreArray(user->Right, &right);
368: /* Scatter values to global vector */
369: DMLocalToGlobalBegin(user->dm, localG, INSERT_VALUES, G);
370: DMLocalToGlobalEnd(user->dm, localG, INSERT_VALUES, G);
372: PetscLogFlops(70.0 * xm * ym);
374: return 0;
375: }
377: /* ------------------------------------------------------------------- */
378: /*
379: FormHessian - Evaluates Hessian matrix.
381: Input Parameters:
382: . tao - the Tao context
383: . x - input vector
384: . ptr - optional user-defined context, as set by TaoSetHessian()
386: Output Parameters:
387: . A - Hessian matrix
388: . B - optionally different preconditioning matrix
390: Notes:
391: Due to mesh point reordering with DMs, we must always work
392: with the local mesh points, and then transform them to the new
393: global numbering with the local-to-global mapping. We cannot work
394: directly with the global numbers for the original uniprocessor mesh!
396: Two methods are available for imposing this transformation
397: when setting matrix entries:
398: (A) MatSetValuesLocal(), using the local ordering (including
399: ghost points!)
400: - Do the following two steps once, before calling TaoSolve()
401: - Use DMGetISLocalToGlobalMapping() to extract the
402: local-to-global map from the DM
403: - Associate this map with the matrix by calling
404: MatSetLocalToGlobalMapping()
405: - Then set matrix entries using the local ordering
406: by calling MatSetValuesLocal()
407: (B) MatSetValues(), using the global ordering
408: - Use DMGetGlobalIndices() to extract the local-to-global map
409: - Then apply this map explicitly yourself
410: - Set matrix entries using the global ordering by calling
411: MatSetValues()
412: Option (A) seems cleaner/easier in many cases, and is the procedure
413: used in this example.
414: */
415: PetscErrorCode FormHessian(Tao tao, Vec X, Mat Hptr, Mat Hessian, void *ptr)
416: {
417: AppCtx *user = (AppCtx *)ptr;
418: PetscInt i, j, k, row;
419: PetscInt mx = user->mx, my = user->my;
420: PetscInt xs, xm, gxs, gxm, ys, ym, gys, gym, col[7];
421: PetscReal hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), hydhx = hy / hx, hxdhy = hx / hy;
422: PetscReal rhx = mx + 1, rhy = my + 1;
423: PetscReal f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb;
424: PetscReal hl, hr, ht, hb, hc, htl, hbr;
425: PetscReal *x, *left, *right, *bottom, *top;
426: PetscReal v[7];
427: Vec localX = user->localX;
428: PetscBool assembled;
430: /* Set various matrix options */
431: MatSetOption(Hessian, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE);
433: /* Initialize matrix entries to zero */
434: MatAssembled(Hessian, &assembled);
435: if (assembled) MatZeroEntries(Hessian);
437: /* Get local mesh boundaries */
438: DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL);
439: DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL);
441: /* Scatter ghost points to local vector */
442: DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, localX);
443: DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, localX);
445: /* Get pointers to vector data */
446: VecGetArray(localX, &x);
447: VecGetArray(user->Top, &top);
448: VecGetArray(user->Bottom, &bottom);
449: VecGetArray(user->Left, &left);
450: VecGetArray(user->Right, &right);
452: /* Compute Hessian over the locally owned part of the mesh */
454: for (i = xs; i < xs + xm; i++) {
455: for (j = ys; j < ys + ym; j++) {
456: row = (j - gys) * gxm + (i - gxs);
458: xc = x[row];
459: xlt = xrb = xl = xr = xb = xt = xc;
461: /* Left side */
462: if (i == gxs) {
463: xl = left[j - ys + 1];
464: xlt = left[j - ys + 2];
465: } else {
466: xl = x[row - 1];
467: }
469: if (j == gys) {
470: xb = bottom[i - xs + 1];
471: xrb = bottom[i - xs + 2];
472: } else {
473: xb = x[row - gxm];
474: }
476: if (i + 1 == gxs + gxm) {
477: xr = right[j - ys + 1];
478: xrb = right[j - ys];
479: } else {
480: xr = x[row + 1];
481: }
483: if (j + 1 == gys + gym) {
484: xt = top[i - xs + 1];
485: xlt = top[i - xs];
486: } else {
487: xt = x[row + gxm];
488: }
490: if (i > gxs && j + 1 < gys + gym) xlt = x[row - 1 + gxm];
491: if (j > gys && i + 1 < gxs + gxm) xrb = x[row + 1 - gxm];
493: d1 = (xc - xl) * rhx;
494: d2 = (xc - xr) * rhx;
495: d3 = (xc - xt) * rhy;
496: d4 = (xc - xb) * rhy;
497: d5 = (xrb - xr) * rhy;
498: d6 = (xrb - xb) * rhx;
499: d7 = (xlt - xl) * rhy;
500: d8 = (xlt - xt) * rhx;
502: f1 = PetscSqrtScalar(1.0 + d1 * d1 + d7 * d7);
503: f2 = PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
504: f3 = PetscSqrtScalar(1.0 + d3 * d3 + d8 * d8);
505: f4 = PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
506: f5 = PetscSqrtScalar(1.0 + d2 * d2 + d5 * d5);
507: f6 = PetscSqrtScalar(1.0 + d4 * d4 + d6 * d6);
509: hl = (-hydhx * (1.0 + d7 * d7) + d1 * d7) / (f1 * f1 * f1) + (-hydhx * (1.0 + d4 * d4) + d1 * d4) / (f2 * f2 * f2);
510: hr = (-hydhx * (1.0 + d5 * d5) + d2 * d5) / (f5 * f5 * f5) + (-hydhx * (1.0 + d3 * d3) + d2 * d3) / (f4 * f4 * f4);
511: ht = (-hxdhy * (1.0 + d8 * d8) + d3 * d8) / (f3 * f3 * f3) + (-hxdhy * (1.0 + d2 * d2) + d2 * d3) / (f4 * f4 * f4);
512: hb = (-hxdhy * (1.0 + d6 * d6) + d4 * d6) / (f6 * f6 * f6) + (-hxdhy * (1.0 + d1 * d1) + d1 * d4) / (f2 * f2 * f2);
514: hbr = -d2 * d5 / (f5 * f5 * f5) - d4 * d6 / (f6 * f6 * f6);
515: htl = -d1 * d7 / (f1 * f1 * f1) - d3 * d8 / (f3 * f3 * f3);
517: 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);
519: hl *= 0.5;
520: hr *= 0.5;
521: ht *= 0.5;
522: hb *= 0.5;
523: hbr *= 0.5;
524: htl *= 0.5;
525: hc *= 0.5;
527: k = 0;
528: if (j > 0) {
529: v[k] = hb;
530: col[k] = row - gxm;
531: k++;
532: }
534: if (j > 0 && i < mx - 1) {
535: v[k] = hbr;
536: col[k] = row - gxm + 1;
537: k++;
538: }
540: if (i > 0) {
541: v[k] = hl;
542: col[k] = row - 1;
543: k++;
544: }
546: v[k] = hc;
547: col[k] = row;
548: k++;
550: if (i < mx - 1) {
551: v[k] = hr;
552: col[k] = row + 1;
553: k++;
554: }
556: if (i > 0 && j < my - 1) {
557: v[k] = htl;
558: col[k] = row + gxm - 1;
559: k++;
560: }
562: if (j < my - 1) {
563: v[k] = ht;
564: col[k] = row + gxm;
565: k++;
566: }
568: /*
569: Set matrix values using local numbering, which was defined
570: earlier, in the main routine.
571: */
572: MatSetValuesLocal(Hessian, 1, &row, k, col, v, INSERT_VALUES);
573: }
574: }
576: /* Restore vectors */
577: VecRestoreArray(localX, &x);
578: VecRestoreArray(user->Left, &left);
579: VecRestoreArray(user->Top, &top);
580: VecRestoreArray(user->Bottom, &bottom);
581: VecRestoreArray(user->Right, &right);
583: /* Assemble the matrix */
584: MatAssemblyBegin(Hessian, MAT_FINAL_ASSEMBLY);
585: MatAssemblyEnd(Hessian, MAT_FINAL_ASSEMBLY);
587: PetscLogFlops(199.0 * xm * ym);
588: return 0;
589: }
591: /* ------------------------------------------------------------------- */
592: /*
593: MSA_BoundaryConditions - Calculates the boundary conditions for
594: the region.
596: Input Parameter:
597: . user - user-defined application context
599: Output Parameter:
600: . user - user-defined application context
601: */
602: static PetscErrorCode MSA_BoundaryConditions(AppCtx *user)
603: {
604: PetscInt i, j, k, maxits = 5, limit = 0;
605: PetscInt xs, ys, xm, ym, gxs, gys, gxm, gym;
606: PetscInt mx = user->mx, my = user->my;
607: PetscInt bsize = 0, lsize = 0, tsize = 0, rsize = 0;
608: PetscReal one = 1.0, two = 2.0, three = 3.0, scl = 1.0, tol = 1e-10;
609: PetscReal fnorm, det, hx, hy, xt = 0, yt = 0;
610: PetscReal u1, u2, nf1, nf2, njac11, njac12, njac21, njac22;
611: PetscReal b = -0.5, t = 0.5, l = -0.5, r = 0.5;
612: PetscReal *boundary;
613: PetscBool flg;
614: Vec Bottom, Top, Right, Left;
616: /* Get local mesh boundaries */
617: DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL);
618: DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL);
620: bsize = xm + 2;
621: lsize = ym + 2;
622: rsize = ym + 2;
623: tsize = xm + 2;
625: VecCreateMPI(MPI_COMM_WORLD, bsize, PETSC_DECIDE, &Bottom);
626: VecCreateMPI(MPI_COMM_WORLD, tsize, PETSC_DECIDE, &Top);
627: VecCreateMPI(MPI_COMM_WORLD, lsize, PETSC_DECIDE, &Left);
628: VecCreateMPI(MPI_COMM_WORLD, rsize, PETSC_DECIDE, &Right);
630: user->Top = Top;
631: user->Left = Left;
632: user->Bottom = Bottom;
633: user->Right = Right;
635: hx = (r - l) / (mx + 1);
636: hy = (t - b) / (my + 1);
638: for (j = 0; j < 4; j++) {
639: if (j == 0) {
640: yt = b;
641: xt = l + hx * xs;
642: limit = bsize;
643: VecGetArray(Bottom, &boundary);
644: } else if (j == 1) {
645: yt = t;
646: xt = l + hx * xs;
647: limit = tsize;
648: VecGetArray(Top, &boundary);
649: } else if (j == 2) {
650: yt = b + hy * ys;
651: xt = l;
652: limit = lsize;
653: VecGetArray(Left, &boundary);
654: } else if (j == 3) {
655: yt = b + hy * ys;
656: xt = r;
657: limit = rsize;
658: VecGetArray(Right, &boundary);
659: }
661: for (i = 0; i < limit; i++) {
662: u1 = xt;
663: u2 = -yt;
664: for (k = 0; k < maxits; k++) {
665: nf1 = u1 + u1 * u2 * u2 - u1 * u1 * u1 / three - xt;
666: nf2 = -u2 - u1 * u1 * u2 + u2 * u2 * u2 / three - yt;
667: fnorm = PetscSqrtScalar(nf1 * nf1 + nf2 * nf2);
668: if (fnorm <= tol) break;
669: njac11 = one + u2 * u2 - u1 * u1;
670: njac12 = two * u1 * u2;
671: njac21 = -two * u1 * u2;
672: njac22 = -one - u1 * u1 + u2 * u2;
673: det = njac11 * njac22 - njac21 * njac12;
674: u1 = u1 - (njac22 * nf1 - njac12 * nf2) / det;
675: u2 = u2 - (njac11 * nf2 - njac21 * nf1) / det;
676: }
678: boundary[i] = u1 * u1 - u2 * u2;
679: if (j == 0 || j == 1) {
680: xt = xt + hx;
681: } else if (j == 2 || j == 3) {
682: yt = yt + hy;
683: }
684: }
685: if (j == 0) {
686: VecRestoreArray(Bottom, &boundary);
687: } else if (j == 1) {
688: VecRestoreArray(Top, &boundary);
689: } else if (j == 2) {
690: VecRestoreArray(Left, &boundary);
691: } else if (j == 3) {
692: VecRestoreArray(Right, &boundary);
693: }
694: }
696: /* Scale the boundary if desired */
698: PetscOptionsGetReal(NULL, NULL, "-bottom", &scl, &flg);
699: if (flg) VecScale(Bottom, scl);
700: PetscOptionsGetReal(NULL, NULL, "-top", &scl, &flg);
701: if (flg) VecScale(Top, scl);
702: PetscOptionsGetReal(NULL, NULL, "-right", &scl, &flg);
703: if (flg) VecScale(Right, scl);
705: PetscOptionsGetReal(NULL, NULL, "-left", &scl, &flg);
706: if (flg) VecScale(Left, scl);
707: return 0;
708: }
710: /* ------------------------------------------------------------------- */
711: /*
712: MSA_Plate - Calculates an obstacle for surface to stretch over.
714: Input Parameter:
715: . user - user-defined application context
717: Output Parameter:
718: . user - user-defined application context
719: */
720: static PetscErrorCode MSA_Plate(Vec XL, Vec XU, void *ctx)
721: {
722: AppCtx *user = (AppCtx *)ctx;
723: PetscInt i, j, row;
724: PetscInt xs, ys, xm, ym;
725: PetscInt mx = user->mx, my = user->my, bmy, bmx;
726: PetscReal t1, t2, t3;
727: PetscReal *xl, lb = PETSC_NINFINITY, ub = PETSC_INFINITY;
728: PetscBool cylinder;
730: user->bmy = PetscMax(0, user->bmy);
731: user->bmy = PetscMin(my, user->bmy);
732: user->bmx = PetscMax(0, user->bmx);
733: user->bmx = PetscMin(mx, user->bmx);
734: bmy = user->bmy;
735: bmx = user->bmx;
737: DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL);
739: VecSet(XL, lb);
740: VecSet(XU, ub);
742: VecGetArray(XL, &xl);
744: PetscOptionsHasName(NULL, NULL, "-cylinder", &cylinder);
745: /* Compute the optional lower box */
746: if (cylinder) {
747: for (i = xs; i < xs + xm; i++) {
748: for (j = ys; j < ys + ym; j++) {
749: row = (j - ys) * xm + (i - xs);
750: t1 = (2.0 * i - mx) * bmy;
751: t2 = (2.0 * j - my) * bmx;
752: t3 = bmx * bmx * bmy * bmy;
753: if (t1 * t1 + t2 * t2 <= t3) xl[row] = user->bheight;
754: }
755: }
756: } else {
757: /* Compute the optional lower box */
758: for (i = xs; i < xs + xm; i++) {
759: for (j = ys; j < ys + ym; j++) {
760: row = (j - ys) * xm + (i - xs);
761: if (i >= (mx - bmx) / 2 && i < mx - (mx - bmx) / 2 && j >= (my - bmy) / 2 && j < my - (my - bmy) / 2) xl[row] = user->bheight;
762: }
763: }
764: }
765: VecRestoreArray(XL, &xl);
767: return 0;
768: }
770: /* ------------------------------------------------------------------- */
771: /*
772: MSA_InitialPoint - Calculates the initial guess in one of three ways.
774: Input Parameters:
775: . user - user-defined application context
776: . X - vector for initial guess
778: Output Parameters:
779: . X - newly computed initial guess
780: */
781: static PetscErrorCode MSA_InitialPoint(AppCtx *user, Vec X)
782: {
783: PetscInt start = -1, i, j;
784: PetscReal zero = 0.0;
785: PetscBool flg;
787: PetscOptionsGetInt(NULL, NULL, "-start", &start, &flg);
788: if (flg && start == 0) { /* The zero vector is reasonable */
789: VecSet(X, zero);
790: } else if (flg && start > 0) { /* Try a random start between -0.5 and 0.5 */
791: PetscRandom rctx;
792: PetscReal np5 = -0.5;
794: PetscRandomCreate(MPI_COMM_WORLD, &rctx);
795: for (i = 0; i < start; i++) VecSetRandom(X, rctx);
796: PetscRandomDestroy(&rctx);
797: VecShift(X, np5);
799: } else { /* Take an average of the boundary conditions */
801: PetscInt row, xs, xm, gxs, gxm, ys, ym, gys, gym;
802: PetscInt mx = user->mx, my = user->my;
803: PetscReal *x, *left, *right, *bottom, *top;
804: Vec localX = user->localX;
806: /* Get local mesh boundaries */
807: DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL);
808: DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL);
810: /* Get pointers to vector data */
811: VecGetArray(user->Top, &top);
812: VecGetArray(user->Bottom, &bottom);
813: VecGetArray(user->Left, &left);
814: VecGetArray(user->Right, &right);
816: VecGetArray(localX, &x);
817: /* Perform local computations */
818: for (j = ys; j < ys + ym; j++) {
819: for (i = xs; i < xs + xm; i++) {
820: row = (j - gys) * gxm + (i - gxs);
821: 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;
822: }
823: }
825: /* Restore vectors */
826: VecRestoreArray(localX, &x);
828: VecRestoreArray(user->Left, &left);
829: VecRestoreArray(user->Top, &top);
830: VecRestoreArray(user->Bottom, &bottom);
831: VecRestoreArray(user->Right, &right);
833: /* Scatter values into global vector */
834: DMLocalToGlobalBegin(user->dm, localX, INSERT_VALUES, X);
835: DMLocalToGlobalEnd(user->dm, localX, INSERT_VALUES, X);
836: }
837: return 0;
838: }
840: /* For testing matrix free submatrices */
841: PetscErrorCode MatrixFreeHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ptr)
842: {
843: AppCtx *user = (AppCtx *)ptr;
844: FormHessian(tao, x, user->H, user->H, ptr);
845: return 0;
846: }
847: PetscErrorCode MyMatMult(Mat H_shell, Vec X, Vec Y)
848: {
849: void *ptr;
850: AppCtx *user;
851: MatShellGetContext(H_shell, &ptr);
852: user = (AppCtx *)ptr;
853: MatMult(user->H, X, Y);
854: return 0;
855: }
857: /*TEST
859: build:
860: requires: !complex
862: test:
863: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type tron -tao_gatol 1.e-5
864: requires: !single
866: test:
867: suffix: 2
868: nsize: 2
869: args: -tao_smonitor -mx 8 -my 8 -bmx 2 -bmy 5 -bheight 0.3 -tao_type blmvm -tao_gatol 1.e-4
870: requires: !single
872: test:
873: suffix: 3
874: nsize: 3
875: args: -tao_smonitor -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_type tron -tao_gatol 1.e-5
876: requires: !single
878: test:
879: suffix: 4
880: nsize: 3
881: 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
882: requires: !single
884: test:
885: suffix: 5
886: nsize: 3
887: 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
888: requires: !single
890: test:
891: suffix: 6
892: nsize: 3
893: 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
894: requires: !single
896: test:
897: suffix: 7
898: nsize: 3
899: 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
900: requires: !single
902: test:
903: suffix: 8
904: 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
905: requires: !single
907: test:
908: suffix: 9
909: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bncg -tao_gatol 1e-4
910: requires: !single
912: test:
913: suffix: 10
914: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bnls -tao_gatol 1e-5
915: requires: !single
917: test:
918: suffix: 11
919: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntr -tao_gatol 1e-5
920: requires: !single
922: test:
923: suffix: 12
924: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntl -tao_gatol 1e-5
925: requires: !single
927: test:
928: suffix: 13
929: 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
930: requires: !single
932: test:
933: suffix: 14
934: 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
935: requires: !single
937: test:
938: suffix: 15
939: 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
940: requires: !single
942: test:
943: suffix: 16
944: args: -tao_smonitor -mx 8 -my 8 -bmx 2 -bmy 5 -bheight 0.3 -tao_gatol 1e-4 -tao_type bqnls
945: requires: !single
947: test:
948: suffix: 17
949: 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
950: requires: !single
952: test:
953: suffix: 18
954: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bnls -tao_gatol 1e-5 -tao_mf_hessian
955: requires: !single
957: test:
958: suffix: 19
959: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntr -tao_gatol 1e-5 -tao_mf_hessian
960: requires: !single
962: test:
963: suffix: 20
964: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntl -tao_gatol 1e-5 -tao_mf_hessian
966: TEST*/