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