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: -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
17: -my <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, FormFunctionGradient()
24: 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: PetscErrorCode QuadraticH(AppCtx *, Vec, Mat);
38: PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);
39: PetscErrorCode FormGradient(Tao, Vec, Vec, void *);
40: PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void *);
41: PetscErrorCode My_Monitor(Tao, void *);
43: int main(int argc, char **argv)
44: {
45: PetscInt Nx, Ny; /* number of processors in x- and y- directions */
46: Vec x; /* solution, gradient vectors */
47: PetscBool flg, 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 */
56: PetscInitialize(&argc, &argv, (char *)0, help);
58: /* Specify dimension of the problem */
59: user.mx = 10;
60: user.my = 10;
62: /* Check for any command line arguments that override defaults */
63: PetscOptionsGetInt(NULL, NULL, "-mx", &user.mx, &flg);
64: PetscOptionsGetInt(NULL, NULL, "-my", &user.my, &flg);
66: PetscPrintf(MPI_COMM_WORLD, "\n---- Minimum Surface Area Problem -----\n");
67: PetscPrintf(MPI_COMM_WORLD, "mx: %" PetscInt_FMT " my: %" PetscInt_FMT " \n\n", user.mx, user.my);
69: /* Let PETSc determine the vector distribution */
70: Nx = PETSC_DECIDE;
71: Ny = PETSC_DECIDE;
73: /* Create distributed array (DM) to manage parallel grid and vectors */
74: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, user.mx, user.my, Nx, Ny, 1, 1, NULL, NULL, &user.dm);
75: DMSetFromOptions(user.dm);
76: DMSetUp(user.dm);
78: /* Create TAO solver and set desired solution method.*/
79: TaoCreate(PETSC_COMM_WORLD, &tao);
80: TaoSetType(tao, TAOCG);
82: /*
83: Extract global vector from DA for the vector of variables -- PETSC routine
84: Compute the initial solution -- application specific, see below
85: Set this vector for use by TAO -- TAO routine
86: */
87: DMCreateGlobalVector(user.dm, &x);
88: MSA_BoundaryConditions(&user);
89: MSA_InitialPoint(&user, x);
90: TaoSetSolution(tao, x);
92: /*
93: Initialize the Application context for use in function evaluations -- application specific, see below.
94: Set routines for function and gradient evaluation
95: */
96: TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user);
98: /*
99: Given the command line arguments, calculate the hessian with either the user-
100: provided function FormHessian, or the default finite-difference driven Hessian
101: functions
102: */
103: PetscOptionsHasName(NULL, NULL, "-tao_fddefault", &fddefault);
104: PetscOptionsHasName(NULL, NULL, "-tao_fdcoloring", &fdcoloring);
106: /*
107: Create a matrix data structure to store the Hessian and set
108: the Hessian evaluation routine.
109: Set the matrix structure to be used for Hessian evaluations
110: */
111: DMCreateMatrix(user.dm, &user.H);
112: MatSetOption(user.H, MAT_SYMMETRIC, PETSC_TRUE);
114: if (fdcoloring) {
115: DMCreateColoring(user.dm, IS_COLORING_GLOBAL, &iscoloring);
116: MatFDColoringCreate(user.H, iscoloring, &matfdcoloring);
117: ISColoringDestroy(&iscoloring);
118: MatFDColoringSetFunction(matfdcoloring, (PetscErrorCode(*)(void))FormGradient, (void *)&user);
119: MatFDColoringSetFromOptions(matfdcoloring);
120: TaoSetHessian(tao, user.H, user.H, TaoDefaultComputeHessianColor, (void *)matfdcoloring);
121: } else if (fddefault) {
122: TaoSetHessian(tao, user.H, user.H, TaoDefaultComputeHessian, (void *)NULL);
123: } else {
124: TaoSetHessian(tao, user.H, user.H, FormHessian, (void *)&user);
125: }
127: /*
128: If my_monitor option is in command line, then use the user-provided
129: monitoring function
130: */
131: PetscOptionsHasName(NULL, NULL, "-my_monitor", &viewmat);
132: if (viewmat) TaoSetMonitor(tao, My_Monitor, NULL, NULL);
134: /* Check for any tao command line options */
135: TaoSetFromOptions(tao);
137: /* SOLVE THE APPLICATION */
138: TaoSolve(tao);
140: TaoView(tao, PETSC_VIEWER_STDOUT_WORLD);
142: /* Free TAO data structures */
143: TaoDestroy(&tao);
145: /* Free PETSc data structures */
146: VecDestroy(&x);
147: MatDestroy(&user.H);
148: if (fdcoloring) MatFDColoringDestroy(&matfdcoloring);
149: PetscFree(user.bottom);
150: PetscFree(user.top);
151: PetscFree(user.left);
152: PetscFree(user.right);
153: DMDestroy(&user.dm);
154: PetscFinalize();
155: return 0;
156: }
158: PetscErrorCode FormGradient(Tao tao, Vec X, Vec G, void *userCtx)
159: {
160: PetscReal fcn;
162: FormFunctionGradient(tao, X, &fcn, G, userCtx);
163: return 0;
164: }
166: /* -------------------------------------------------------------------- */
167: /* FormFunctionGradient - Evaluates the function and corresponding gradient.
169: Input Parameters:
170: . tao - the Tao context
171: . XX - input vector
172: . userCtx - optional user-defined context, as set by TaoSetObjectiveAndGradient()
174: Output Parameters:
175: . fcn - the newly evaluated function
176: . GG - vector containing the newly evaluated gradient
177: */
178: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn, Vec G, void *userCtx)
179: {
180: AppCtx *user = (AppCtx *)userCtx;
181: PetscInt i, j;
182: PetscInt mx = user->mx, my = user->my;
183: PetscInt xs, xm, gxs, gxm, ys, ym, gys, gym;
184: PetscReal ft = 0.0;
185: PetscReal hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), hydhx = hy / hx, hxdhy = hx / hy, area = 0.5 * hx * hy;
186: PetscReal rhx = mx + 1, rhy = my + 1;
187: PetscReal f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb;
188: PetscReal df1dxc, df2dxc, df3dxc, df4dxc, df5dxc, df6dxc;
189: PetscReal **g, **x;
190: Vec localX;
192: /* Get local mesh boundaries */
193: DMGetLocalVector(user->dm, &localX);
194: DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL);
195: DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL);
197: /* Scatter ghost points to local vector */
198: DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, localX);
199: DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, localX);
201: /* Get pointers to vector data */
202: DMDAVecGetArray(user->dm, localX, (void **)&x);
203: DMDAVecGetArray(user->dm, G, (void **)&g);
205: /* Compute function and gradient over the locally owned part of the mesh */
206: for (j = ys; j < ys + ym; j++) {
207: for (i = xs; i < xs + xm; i++) {
208: xc = x[j][i];
209: xlt = xrb = xl = xr = xb = xt = xc;
211: if (i == 0) { /* left side */
212: xl = user->left[j - ys + 1];
213: xlt = user->left[j - ys + 2];
214: } else {
215: xl = x[j][i - 1];
216: }
218: if (j == 0) { /* bottom side */
219: xb = user->bottom[i - xs + 1];
220: xrb = user->bottom[i - xs + 2];
221: } else {
222: xb = x[j - 1][i];
223: }
225: if (i + 1 == gxs + gxm) { /* right side */
226: xr = user->right[j - ys + 1];
227: xrb = user->right[j - ys];
228: } else {
229: xr = x[j][i + 1];
230: }
232: if (j + 1 == gys + gym) { /* top side */
233: xt = user->top[i - xs + 1];
234: xlt = user->top[i - xs];
235: } else {
236: xt = x[j + 1][i];
237: }
239: if (i > gxs && j + 1 < gys + gym) xlt = x[j + 1][i - 1];
240: if (j > gys && i + 1 < gxs + gxm) xrb = x[j - 1][i + 1];
242: d1 = (xc - xl);
243: d2 = (xc - xr);
244: d3 = (xc - xt);
245: d4 = (xc - xb);
246: d5 = (xr - xrb);
247: d6 = (xrb - xb);
248: d7 = (xlt - xl);
249: d8 = (xt - xlt);
251: df1dxc = d1 * hydhx;
252: df2dxc = (d1 * hydhx + d4 * hxdhy);
253: df3dxc = d3 * hxdhy;
254: df4dxc = (d2 * hydhx + d3 * hxdhy);
255: df5dxc = d2 * hydhx;
256: df6dxc = d4 * hxdhy;
258: d1 *= rhx;
259: d2 *= rhx;
260: d3 *= rhy;
261: d4 *= rhy;
262: d5 *= rhy;
263: d6 *= rhx;
264: d7 *= rhy;
265: d8 *= rhx;
267: f1 = PetscSqrtReal(1.0 + d1 * d1 + d7 * d7);
268: f2 = PetscSqrtReal(1.0 + d1 * d1 + d4 * d4);
269: f3 = PetscSqrtReal(1.0 + d3 * d3 + d8 * d8);
270: f4 = PetscSqrtReal(1.0 + d3 * d3 + d2 * d2);
271: f5 = PetscSqrtReal(1.0 + d2 * d2 + d5 * d5);
272: f6 = PetscSqrtReal(1.0 + d4 * d4 + d6 * d6);
274: ft = ft + (f2 + f4);
276: df1dxc /= f1;
277: df2dxc /= f2;
278: df3dxc /= f3;
279: df4dxc /= f4;
280: df5dxc /= f5;
281: df6dxc /= f6;
283: g[j][i] = (df1dxc + df2dxc + df3dxc + df4dxc + df5dxc + df6dxc) * 0.5;
284: }
285: }
287: /* Compute triangular areas along the border of the domain. */
288: if (xs == 0) { /* left side */
289: for (j = ys; j < ys + ym; j++) {
290: d3 = (user->left[j - ys + 1] - user->left[j - ys + 2]) * rhy;
291: d2 = (user->left[j - ys + 1] - x[j][0]) * rhx;
292: ft = ft + PetscSqrtReal(1.0 + d3 * d3 + d2 * d2);
293: }
294: }
295: if (ys == 0) { /* bottom side */
296: for (i = xs; i < xs + xm; i++) {
297: d2 = (user->bottom[i + 1 - xs] - user->bottom[i - xs + 2]) * rhx;
298: d3 = (user->bottom[i - xs + 1] - x[0][i]) * rhy;
299: ft = ft + PetscSqrtReal(1.0 + d3 * d3 + d2 * d2);
300: }
301: }
303: if (xs + xm == mx) { /* right side */
304: for (j = ys; j < ys + ym; j++) {
305: d1 = (x[j][mx - 1] - user->right[j - ys + 1]) * rhx;
306: d4 = (user->right[j - ys] - user->right[j - ys + 1]) * rhy;
307: ft = ft + PetscSqrtReal(1.0 + d1 * d1 + d4 * d4);
308: }
309: }
310: if (ys + ym == my) { /* top side */
311: for (i = xs; i < xs + xm; i++) {
312: d1 = (x[my - 1][i] - user->top[i - xs + 1]) * rhy;
313: d4 = (user->top[i - xs + 1] - user->top[i - xs]) * rhx;
314: ft = ft + PetscSqrtReal(1.0 + d1 * d1 + d4 * d4);
315: }
316: }
318: if (ys == 0 && xs == 0) {
319: d1 = (user->left[0] - user->left[1]) / hy;
320: d2 = (user->bottom[0] - user->bottom[1]) * rhx;
321: ft += PetscSqrtReal(1.0 + d1 * d1 + d2 * d2);
322: }
323: if (ys + ym == my && xs + xm == mx) {
324: d1 = (user->right[ym + 1] - user->right[ym]) * rhy;
325: d2 = (user->top[xm + 1] - user->top[xm]) * rhx;
326: ft += PetscSqrtReal(1.0 + d1 * d1 + d2 * d2);
327: }
329: ft = ft * area;
330: MPI_Allreduce(&ft, fcn, 1, MPIU_REAL, MPIU_SUM, MPI_COMM_WORLD);
332: /* Restore vectors */
333: DMDAVecRestoreArray(user->dm, localX, (void **)&x);
334: DMDAVecRestoreArray(user->dm, G, (void **)&g);
336: /* Scatter values to global vector */
337: DMRestoreLocalVector(user->dm, &localX);
338: PetscLogFlops(67.0 * xm * ym);
339: return 0;
340: }
342: /* ------------------------------------------------------------------- */
343: /*
344: FormHessian - Evaluates Hessian matrix.
346: Input Parameters:
347: . tao - the Tao context
348: . x - input vector
349: . ptr - optional user-defined context, as set by TaoSetHessian()
351: Output Parameters:
352: . H - Hessian matrix
353: . Hpre - optionally different preconditioning matrix
354: . flg - flag indicating matrix structure
356: */
357: PetscErrorCode FormHessian(Tao tao, Vec X, Mat H, Mat Hpre, void *ptr)
358: {
359: AppCtx *user = (AppCtx *)ptr;
361: /* Evaluate the Hessian entries*/
362: QuadraticH(user, X, H);
363: return 0;
364: }
366: /* ------------------------------------------------------------------- */
367: /*
368: QuadraticH - Evaluates Hessian matrix.
370: Input Parameters:
371: . user - user-defined context, as set by TaoSetHessian()
372: . X - input vector
374: Output Parameter:
375: . H - Hessian matrix
376: */
377: PetscErrorCode QuadraticH(AppCtx *user, Vec X, Mat Hessian)
378: {
379: PetscInt i, j, k;
380: PetscInt mx = user->mx, my = user->my;
381: PetscInt xs, xm, gxs, gxm, ys, ym, gys, gym;
382: PetscReal hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), hydhx = hy / hx, hxdhy = hx / hy;
383: PetscReal f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb;
384: PetscReal hl, hr, ht, hb, hc, htl, hbr;
385: PetscReal **x, v[7];
386: MatStencil col[7], row;
387: Vec localX;
388: PetscBool assembled;
390: /* Get local mesh boundaries */
391: DMGetLocalVector(user->dm, &localX);
393: DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL);
394: DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL);
396: /* Scatter ghost points to local vector */
397: DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, localX);
398: DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, localX);
400: /* Get pointers to vector data */
401: DMDAVecGetArray(user->dm, localX, (void **)&x);
403: /* Initialize matrix entries to zero */
404: MatAssembled(Hessian, &assembled);
405: if (assembled) MatZeroEntries(Hessian);
407: /* Set various matrix options */
408: MatSetOption(Hessian, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE);
410: /* Compute Hessian over the locally owned part of the mesh */
412: for (j = ys; j < ys + ym; j++) {
413: for (i = xs; i < xs + xm; i++) {
414: xc = x[j][i];
415: xlt = xrb = xl = xr = xb = xt = xc;
417: /* Left side */
418: if (i == 0) {
419: xl = user->left[j - ys + 1];
420: xlt = user->left[j - ys + 2];
421: } else {
422: xl = x[j][i - 1];
423: }
425: if (j == 0) {
426: xb = user->bottom[i - xs + 1];
427: xrb = user->bottom[i - xs + 2];
428: } else {
429: xb = x[j - 1][i];
430: }
432: if (i + 1 == mx) {
433: xr = user->right[j - ys + 1];
434: xrb = user->right[j - ys];
435: } else {
436: xr = x[j][i + 1];
437: }
439: if (j + 1 == my) {
440: xt = user->top[i - xs + 1];
441: xlt = user->top[i - xs];
442: } else {
443: xt = x[j + 1][i];
444: }
446: if (i > 0 && j + 1 < my) xlt = x[j + 1][i - 1];
447: if (j > 0 && i + 1 < mx) xrb = x[j - 1][i + 1];
449: d1 = (xc - xl) / hx;
450: d2 = (xc - xr) / hx;
451: d3 = (xc - xt) / hy;
452: d4 = (xc - xb) / hy;
453: d5 = (xrb - xr) / hy;
454: d6 = (xrb - xb) / hx;
455: d7 = (xlt - xl) / hy;
456: d8 = (xlt - xt) / hx;
458: f1 = PetscSqrtReal(1.0 + d1 * d1 + d7 * d7);
459: f2 = PetscSqrtReal(1.0 + d1 * d1 + d4 * d4);
460: f3 = PetscSqrtReal(1.0 + d3 * d3 + d8 * d8);
461: f4 = PetscSqrtReal(1.0 + d3 * d3 + d2 * d2);
462: f5 = PetscSqrtReal(1.0 + d2 * d2 + d5 * d5);
463: f6 = PetscSqrtReal(1.0 + d4 * d4 + d6 * d6);
465: hl = (-hydhx * (1.0 + d7 * d7) + d1 * d7) / (f1 * f1 * f1) + (-hydhx * (1.0 + d4 * d4) + d1 * d4) / (f2 * f2 * f2);
466: hr = (-hydhx * (1.0 + d5 * d5) + d2 * d5) / (f5 * f5 * f5) + (-hydhx * (1.0 + d3 * d3) + d2 * d3) / (f4 * f4 * f4);
467: ht = (-hxdhy * (1.0 + d8 * d8) + d3 * d8) / (f3 * f3 * f3) + (-hxdhy * (1.0 + d2 * d2) + d2 * d3) / (f4 * f4 * f4);
468: hb = (-hxdhy * (1.0 + d6 * d6) + d4 * d6) / (f6 * f6 * f6) + (-hxdhy * (1.0 + d1 * d1) + d1 * d4) / (f2 * f2 * f2);
470: hbr = -d2 * d5 / (f5 * f5 * f5) - d4 * d6 / (f6 * f6 * f6);
471: htl = -d1 * d7 / (f1 * f1 * f1) - d3 * d8 / (f3 * f3 * f3);
473: 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);
475: hl /= 2.0;
476: hr /= 2.0;
477: ht /= 2.0;
478: hb /= 2.0;
479: hbr /= 2.0;
480: htl /= 2.0;
481: hc /= 2.0;
483: row.j = j;
484: row.i = i;
485: k = 0;
486: if (j > 0) {
487: v[k] = hb;
488: col[k].j = j - 1;
489: col[k].i = i;
490: k++;
491: }
493: if (j > 0 && i < mx - 1) {
494: v[k] = hbr;
495: col[k].j = j - 1;
496: col[k].i = i + 1;
497: k++;
498: }
500: if (i > 0) {
501: v[k] = hl;
502: col[k].j = j;
503: col[k].i = i - 1;
504: k++;
505: }
507: v[k] = hc;
508: col[k].j = j;
509: col[k].i = i;
510: k++;
512: if (i < mx - 1) {
513: v[k] = hr;
514: col[k].j = j;
515: col[k].i = i + 1;
516: k++;
517: }
519: if (i > 0 && j < my - 1) {
520: v[k] = htl;
521: col[k].j = j + 1;
522: col[k].i = i - 1;
523: k++;
524: }
526: if (j < my - 1) {
527: v[k] = ht;
528: col[k].j = j + 1;
529: col[k].i = i;
530: k++;
531: }
533: /*
534: Set matrix values using local numbering, which was defined
535: earlier, in the main routine.
536: */
537: MatSetValuesStencil(Hessian, 1, &row, k, col, v, INSERT_VALUES);
538: }
539: }
541: DMDAVecRestoreArray(user->dm, localX, (void **)&x);
542: DMRestoreLocalVector(user->dm, &localX);
544: MatAssemblyBegin(Hessian, MAT_FINAL_ASSEMBLY);
545: MatAssemblyEnd(Hessian, MAT_FINAL_ASSEMBLY);
547: PetscLogFlops(199.0 * xm * ym);
548: return 0;
549: }
551: /* ------------------------------------------------------------------- */
552: /*
553: MSA_BoundaryConditions - Calculates the boundary conditions for
554: the region.
556: Input Parameter:
557: . user - user-defined application context
559: Output Parameter:
560: . user - user-defined application context
561: */
562: static PetscErrorCode MSA_BoundaryConditions(AppCtx *user)
563: {
564: PetscInt i, j, k, limit = 0, maxits = 5;
565: PetscInt xs, ys, xm, ym, gxs, gys, gxm, gym;
566: PetscInt mx = user->mx, my = user->my;
567: PetscInt bsize = 0, lsize = 0, tsize = 0, rsize = 0;
568: PetscReal one = 1.0, two = 2.0, three = 3.0, tol = 1e-10;
569: PetscReal fnorm, det, hx, hy, xt = 0, yt = 0;
570: PetscReal u1, u2, nf1, nf2, njac11, njac12, njac21, njac22;
571: PetscReal b = -0.5, t = 0.5, l = -0.5, r = 0.5;
572: PetscReal *boundary;
573: PetscBool flg;
575: /* Get local mesh boundaries */
576: DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL);
577: DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL);
579: bsize = xm + 2;
580: lsize = ym + 2;
581: rsize = ym + 2;
582: tsize = xm + 2;
584: PetscMalloc1(bsize, &user->bottom);
585: PetscMalloc1(tsize, &user->top);
586: PetscMalloc1(lsize, &user->left);
587: PetscMalloc1(rsize, &user->right);
589: hx = (r - l) / (mx + 1);
590: hy = (t - b) / (my + 1);
592: for (j = 0; j < 4; j++) {
593: if (j == 0) {
594: yt = b;
595: xt = l + hx * xs;
596: limit = bsize;
597: boundary = user->bottom;
598: } else if (j == 1) {
599: yt = t;
600: xt = l + hx * xs;
601: limit = tsize;
602: boundary = user->top;
603: } else if (j == 2) {
604: yt = b + hy * ys;
605: xt = l;
606: limit = lsize;
607: boundary = user->left;
608: } else { /* if (j==3) */
609: yt = b + hy * ys;
610: xt = r;
611: limit = rsize;
612: boundary = user->right;
613: }
615: for (i = 0; i < limit; i++) {
616: u1 = xt;
617: u2 = -yt;
618: for (k = 0; k < maxits; k++) {
619: nf1 = u1 + u1 * u2 * u2 - u1 * u1 * u1 / three - xt;
620: nf2 = -u2 - u1 * u1 * u2 + u2 * u2 * u2 / three - yt;
621: fnorm = PetscSqrtReal(nf1 * nf1 + nf2 * nf2);
622: if (fnorm <= tol) break;
623: njac11 = one + u2 * u2 - u1 * u1;
624: njac12 = two * u1 * u2;
625: njac21 = -two * u1 * u2;
626: njac22 = -one - u1 * u1 + u2 * u2;
627: det = njac11 * njac22 - njac21 * njac12;
628: u1 = u1 - (njac22 * nf1 - njac12 * nf2) / det;
629: u2 = u2 - (njac11 * nf2 - njac21 * nf1) / det;
630: }
632: boundary[i] = u1 * u1 - u2 * u2;
633: if (j == 0 || j == 1) {
634: xt = xt + hx;
635: } else { /* if (j==2 || j==3) */
636: yt = yt + hy;
637: }
638: }
639: }
641: /* Scale the boundary if desired */
642: if (1 == 1) {
643: PetscReal scl = 1.0;
645: PetscOptionsGetReal(NULL, NULL, "-bottom", &scl, &flg);
646: if (flg) {
647: for (i = 0; i < bsize; i++) user->bottom[i] *= scl;
648: }
650: PetscOptionsGetReal(NULL, NULL, "-top", &scl, &flg);
651: if (flg) {
652: for (i = 0; i < tsize; i++) user->top[i] *= scl;
653: }
655: PetscOptionsGetReal(NULL, NULL, "-right", &scl, &flg);
656: if (flg) {
657: for (i = 0; i < rsize; i++) user->right[i] *= scl;
658: }
660: PetscOptionsGetReal(NULL, NULL, "-left", &scl, &flg);
661: if (flg) {
662: for (i = 0; i < lsize; i++) user->left[i] *= scl;
663: }
664: }
665: return 0;
666: }
668: /* ------------------------------------------------------------------- */
669: /*
670: MSA_InitialPoint - Calculates the initial guess in one of three ways.
672: Input Parameters:
673: . user - user-defined application context
674: . X - vector for initial guess
676: Output Parameters:
677: . X - newly computed initial guess
678: */
679: static PetscErrorCode MSA_InitialPoint(AppCtx *user, Vec X)
680: {
681: PetscInt start2 = -1, i, j;
682: PetscReal start1 = 0;
683: PetscBool flg1, flg2;
685: PetscOptionsGetReal(NULL, NULL, "-start", &start1, &flg1);
686: PetscOptionsGetInt(NULL, NULL, "-random", &start2, &flg2);
688: if (flg1) { /* The zero vector is reasonable */
690: VecSet(X, start1);
692: } else if (flg2 && start2 > 0) { /* Try a random start between -0.5 and 0.5 */
693: PetscRandom rctx;
694: PetscReal np5 = -0.5;
696: PetscRandomCreate(PETSC_COMM_WORLD, &rctx);
697: VecSetRandom(X, rctx);
698: PetscRandomDestroy(&rctx);
699: VecShift(X, np5);
701: } else { /* Take an average of the boundary conditions */
702: PetscInt xs, xm, ys, ym;
703: PetscInt mx = user->mx, my = user->my;
704: PetscReal **x;
706: /* Get local mesh boundaries */
707: DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL);
709: /* Get pointers to vector data */
710: DMDAVecGetArray(user->dm, X, (void **)&x);
712: /* Perform local computations */
713: for (j = ys; j < ys + ym; j++) {
714: 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;
715: }
716: DMDAVecRestoreArray(user->dm, X, (void **)&x);
717: PetscLogFlops(9.0 * xm * ym);
718: }
719: return 0;
720: }
722: /*-----------------------------------------------------------------------*/
723: PetscErrorCode My_Monitor(Tao tao, void *ctx)
724: {
725: Vec X;
727: TaoGetSolution(tao, &X);
728: VecView(X, PETSC_VIEWER_STDOUT_WORLD);
729: return 0;
730: }
732: /*TEST
734: build:
735: requires: !complex
737: test:
738: args: -tao_smonitor -tao_type lmvm -mx 10 -my 8 -tao_gatol 1.e-3
739: requires: !single
741: test:
742: suffix: 2
743: nsize: 2
744: args: -tao_smonitor -tao_type nls -tao_nls_ksp_max_it 15 -tao_gatol 1.e-4
745: filter: grep -v "nls ksp"
746: requires: !single
748: test:
749: suffix: 3
750: nsize: 3
751: args: -tao_smonitor -tao_type cg -tao_cg_type fr -mx 10 -my 10 -tao_gatol 1.e-3
752: requires: !single
754: test:
755: suffix: 5
756: nsize: 2
757: args: -tao_smonitor -tao_type bmrm -mx 10 -my 8 -tao_gatol 1.e-3
758: requires: !single
760: TEST*/