Actual source code: minsurf3.c
1: /*$Id: minsurf3.c,v 1.1 2002/05/13 19:32:57 benson Exp $*/
3: /* Program usage: mpirun -np <proc> minsurf3 [-help] [all TAO options] */
5: /* minsurf1: boundary values like in COPS example:
6: - parabola for top and bottom boundaries
7: - zero o.w. */
9: /*
10: Include "tao.h" so we can use TAO solvers.
11: petscda.h for distributed array
12: ad_deriv.h for AD gradient
13: */
15: #include "petscda.h"
16: #include tao.h
17: #include taodaapplication.h
19: static char help[] =
20: "Given a rectangular 2-D domain and boundary values along \n\
21: the edges of the domain, the objective is to find the surface with \n\
22: the minimal area that satisfies the boundary conditions.\n\
23: The command line options are:\n\
24: -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
25: -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
26: -nlevels <nlevels>, where <nlevels> = number of levels in multigrid\n\
27: -byelement, if computation is made by functions on rectangular elements\n\
28: -adic, if AD is used (AD is not used by default)\n\
29: -bottom <b>, where <b> = bottom bound on the rectangular domain\n\
30: -top <t>, where <t> = bottom bound on the rectangular domain\n\
31: -left <l>, where <l> = left bound on the rectangular domain\n\
32: -right <r>, where <r> = right bound on the rectangular domain\n\
33: -cops <r>, if COPS boundaries should be use (default MINPACK bounds)\n\
34: -obst <obst>, where <obst> = 1 is obstacle is present, zero otherwise \n\n";
36: /*T
37: Concepts: TAO - Solving a bounded minimization problem
38: Routines: TaoInitialize(); TaoFinalize();
39: Routines: TaoCreate(); TaoDestroy();
40: Routines: DAAppSetVariableBoundsRoutine();
41: Routines: DAAppSetElementObjectiveAndGradientRoutine();
42: Routines: DAAppSetElementHessianRoutine();
43: Routines: DAAppSetObjectiveAndGradientRoutine();
44: Routines: DAAppSetADElementFunctionGradient();
45: Routines: DAAppSetHessianRoutine();
46: Routines: TaoSetOptions();
47: Routines: TaoAppGetSolutionStatus(); TaoDAAppSolve();
48: Routines: DAAppSetBeforeMonitor(); TaoView();
49: Routines: DAAppGetSolution();
50: Routines: DAAppGetInterpolationMatrix();
51: Processors: n
52: T*/
54: /*
55: User-defined application context - contains data needed by the
56: application-provided call-back routines.
57: */
60: /*
61: This structure is used only when an ADIC generated gradient is used.
62: An InactiveDouble type is a double
63: */
64: typedef struct {
65:
66: InactiveDouble hx, hy; /* increment size in both directions */
67: InactiveDouble area; /* area of the triangles */
69: } ADFGCtx;
70: int ad_MSurfLocalFunction(int[2], DERIV_TYPE[4], DERIV_TYPE*, void*);
72: typedef struct {
73: PetscReal b, t, l, r; /* domain boundaries */
74: PetscReal bheight; /* height of obstacle */
75: PetscReal fx, fy; /* relative size of obstacle */
76: double hx, hy; /* increment size in both directions */
77: double area; /* area of the triangles */
78: int mx, my; /* discretization including boundaries */
79: int bmx, bmy; /* discretization of obstacle */
80: ADFGCtx fgctx; /* Used only when an ADIC generated gradient is used */
81: } AppCtx;
83: /* User-defined routines */
84: static int AppCtxInitialize(void *);
86: static int MSurfLocalFunctionGradient(int[2], double[4], double *, double[4], void *);
87: static int WholeMSurfFunctionGradient(TAO_APPLICATION,DA,Vec,double *,Vec,void*);
89: static int MSurfLocalHessian(int[2], double[4], double[4][4], void *);
90: static int WholeMSurfHessian(TAO_APPLICATION,DA,Vec,Mat,void*);
92: static int COPS_Bounds(TAO_APPLICATION, DA, Vec, Vec, void*);
93: static int MINPACK_Bounds(TAO_APPLICATION, DA, Vec, Vec, void *);
95: static int MyGridMonitorBefore(TAO_APPLICATION, DA, int, void *);
100: int main( int argc, char **argv ) {
102: int info; /* used to check for functions returning nonzeros */
103: int iter,nlevels;
104: int Nx,Ny;
105: double ff,gnorm;
106: DA DAarray[20];
107: Vec X;
108: PetscTruth flg, PreLoad = PETSC_TRUE; /* flags */
109: TaoMethod method = "tao_bnls"; /* minimization method */
110: TaoTerminateReason reason;
111: TAO_SOLVER tao; /* TAO_SOLVER solver context */
112: TAO_APPLICATION MSurfApp; /* The PETSc application */
113: AppCtx user; /* user-defined work context */
115: /* Initialize TAO */
116: PetscInitialize(&argc, &argv, (char *)0, help);
117: TaoInitialize(&argc, &argv, (char *)0, help);
119: PreLoadBegin(PreLoad,"Solve");
121: info = AppCtxInitialize((void*)&user); CHKERRQ(info);
123: nlevels=4;
124: info = PetscOptionsGetInt(PETSC_NULL,"-nlevels",&nlevels,&flg); CHKERRQ(info);
125: if (PreLoadIt == 0) {
126: nlevels = 1; user.mx = 11; user.my = 11;
127: }
129: PetscPrintf(MPI_COMM_WORLD,"\n---- Minimal Surface Problem (simple boundary) -----\n\n");
131: /* Let PETSc determine the vector distribution */
132: Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;
134: /* Create distributed array (DA) to manage parallel grid and vectors */
135: info = DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_BOX,user.mx,
136: user.my,Nx,Ny,1,1,PETSC_NULL,PETSC_NULL,&DAarray[0]); CHKERRQ(info);
137: for (iter=1;iter<nlevels;iter++){
138: info = DARefine(DAarray[iter-1],PETSC_COMM_WORLD,&DAarray[iter]); CHKERRQ(info);
139: }
141: /* Create TAO solver and set desired solution method */
142: info = TaoCreate(MPI_COMM_WORLD,method,&tao); CHKERRQ(info);
143: info = TaoApplicationCreate(PETSC_COMM_WORLD, &MSurfApp); CHKERRQ(info);
144: info = TaoAppSetDAApp(MSurfApp,DAarray,nlevels); CHKERRQ(info);
146: /* Sets routines for function, gradient and bounds evaluation */
147: info = PetscOptionsHasName(TAO_NULL, "-cops", &flg); CHKERRQ(info);
148: if (flg){
149: info = DAAppSetVariableBoundsRoutine(MSurfApp,COPS_Bounds,(void *)&user); CHKERRQ(info);
150: } else {
151: info = DAAppSetVariableBoundsRoutine(MSurfApp,MINPACK_Bounds,(void *)&user); CHKERRQ(info);
152: }
154: info = PetscOptionsHasName(TAO_NULL, "-byelement", &flg); CHKERRQ(info);
155: if (flg) {
157: /* Sets routines for function and gradient evaluation, element by element */
158: info = PetscOptionsHasName(TAO_NULL, "-adic", &flg); CHKERRQ(info);
159: if (flg) {
160: info = DAAppSetADElementFunctionGradient(MSurfApp,ad_MSurfLocalFunction,150,(void *)&user.fgctx); CHKERRQ(info);
161: } else {
162: info = DAAppSetElementObjectiveAndGradientRoutine(MSurfApp,MSurfLocalFunctionGradient,36,(void *)&user); CHKERRQ(info);
163: }
165: /* Sets routines for Hessian evaluation, element by element */
166: info = DAAppSetElementHessianRoutine(MSurfApp,MSurfLocalHessian,87,(void*)&user); CHKERRQ(info);
168: } else {
170: /* Sets routines for function and gradient evaluation, all in one routine */
171: info = DAAppSetObjectiveAndGradientRoutine(MSurfApp,WholeMSurfFunctionGradient,(void *)&user); CHKERRQ(info);
173: /* Sets routines for Hessian evaluation, all in one routine */
174: info = DAAppSetHessianRoutine(MSurfApp,WholeMSurfHessian,(void*)&user); CHKERRQ(info);
175:
176: }
178: info = DAAppSetBeforeMonitor(MSurfApp,MyGridMonitorBefore,(void*)&user); CHKERRQ(info);
179: info = PetscOptionsHasName(TAO_NULL,"-tao_monitor", &flg); CHKERRQ(info);
180: if (flg){
181: info = DAAppPrintStageTimes(MSurfApp); CHKERRQ(info);
182: info = DAAppPrintInterpolationError(MSurfApp); CHKERRQ(info);
183: }
185: info = TaoAppSetRelativeTolerance(MSurfApp,1.0e-8); CHKERRQ(info);
186: info = TaoSetTolerances(tao,0,0,0,0); CHKERRQ(info);
187: info = TaoSetGradientTolerances(tao,0,0,0); CHKERRQ(info);
189: /* Check for any tao command line options */
190: info = TaoSetOptions(MSurfApp,tao); CHKERRQ(info);
192: /* SOLVE THE APPLICATION */
193: info = TaoDAAppSolve(MSurfApp,tao); CHKERRQ(info);
195: /* Get information on termination */
196: info = TaoGetSolutionStatus(tao,&iter,&ff,&gnorm,0,0,&reason); CHKERRQ(info);
197: if (reason <= 0 ){
198: PetscPrintf(MPI_COMM_WORLD,"Try a different TAO method, adjust some parameters, or check the function evaluation routines\n");
199: PetscPrintf(MPI_COMM_WORLD," Iterations: %d, Function Value: %4.2e, Residual: %4.2e \n",iter,ff,gnorm);
200: }
202: info = PetscOptionsHasName(PETSC_NULL,"-view_sol",&flg); CHKERRQ(info);
203: if (flg){
204: info = DAAppGetSolution(MSurfApp,nlevels-1,&X); CHKERRQ(info);
205: info=VecView(X,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(info);
206: }
208: /* To View TAO solver information */
209: // info = TaoView(tao); CHKERRQ(info);
211: /* Free TAO data structures */
212: info = TaoDestroy(tao); CHKERRQ(info);
213: info = TaoAppDestroy(MSurfApp); CHKERRQ(info);
215: /* Free PETSc data structures */
216: for (iter=0;iter<nlevels;iter++){
217: info = DADestroy(DAarray[iter]); CHKERRQ(info);
218: }
220: PreLoadEnd();
222: /* Finalize TAO */
223: TaoFinalize();
224: PetscFinalize();
226: return 0;
227: } /* main */
231: /*----- The following two routines
232: MyGridMonitorBefore MyGridMonitorAfter
233: help diplay info of iterations at every grid level
234: */
238: static int MyGridMonitorBefore(TAO_APPLICATION myapp, DA da, int level, void *ctx) {
240: AppCtx *user = (AppCtx*)ctx;
241: int info,mx,my;
243: info = DAGetInfo(da,PETSC_NULL,&mx,&my,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
244: PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(info);
245: user->mx = mx;
246: user->my = my;
247: user->hx = (user->r - user->l) / (user->mx - 1);
248: user->hy = (user->t - user->b) / (user->my - 1);
249: user->area = 0.5 * user->hx * user->hy;
250: user->fgctx.hx = user->hx;
251: user->fgctx.hy = user->hy;
252: user->fgctx.area = user->area;
254: user->bmx=(int)((mx+1)*user->fx); user->bmy=(int)((my+1)*user->fy);
256: PetscPrintf(MPI_COMM_WORLD,"Grid: %d, mx: %d my: %d \n",level,mx,my);
258: return 0;
259: }
263: /*
264: AppCtxInitialize - Sets initial values for the application context parameters
266: Input:
267: ptr - void user-defined application context
269: Output:
270: ptr - user-defined application context with the default or user-provided
271: parameters
272: */
273: static int AppCtxInitialize(void *ptr) {
275: AppCtx *user = (AppCtx*)ptr;
276: PetscTruth flg;
277: int info;
279: /* Specify default parameters */
280: user->mx = user->my = 11;
281: user->b = -0.5;
282: user->t = 0.5;
283: user->l = -0.5;
284: user->r = 0.5;
285: user->fx=0.5;
286: user->fy=0.5;
287: user->bheight=0.0;
289: /* Check for command line arguments that override defaults */
290: info = PetscOptionsGetInt(TAO_NULL, "-mx", &user->mx, &flg); CHKERRQ(info);
291: info = PetscOptionsGetInt(TAO_NULL, "-my", &user->my, &flg); CHKERRQ(info);
292: info = PetscOptionsGetReal(TAO_NULL, "-bottom", &user->b, &flg); CHKERRQ(info);
293: info = PetscOptionsGetReal(TAO_NULL, "-top", &user->t, &flg); CHKERRQ(info);
294: info = PetscOptionsGetReal(TAO_NULL, "-left", &user->l, &flg); CHKERRQ(info);
295: info = PetscOptionsGetReal(TAO_NULL, "-right", &user->r, &flg); CHKERRQ(info);
296: info = PetscOptionsGetReal(PETSC_NULL,"-bmx",&user->fx,&flg); CHKERRQ(info);
297: info = PetscOptionsGetReal(PETSC_NULL,"-bmy",&user->fy,&flg); CHKERRQ(info);
298: info = PetscOptionsGetReal(PETSC_NULL,"-bheight",&user->bheight,&flg); CHKERRQ(info);
300: user->hx = (user->r - user->l) / (user->mx - 1);
301: user->hy = (user->t - user->b) / (user->my - 1);
302: user->area = 0.5 * user->hx * user->hy;
303: info = PetscLogFlops(8); CHKERRQ(info);
305: return 0;
307: } /* AppCtxInitialize */
310: /*------- USER-DEFINED: set the upper and lower bounds for the variables -------*/
314: /*
315: FormBounds - Forms bounds on the variables
317: Input:
318: user - user-defined application context
320: Output:
321: XL - vector of lower bounds
322: XU - vector of upper bounds
324: */
325: static int COPS_Bounds(TAO_APPLICATION tao, DA da, Vec XL, Vec XU, void *ptr) {
327: AppCtx *user = (AppCtx*)ptr;
328: PetscTruth flg;
329: int i, j, mx, my, info, xs, xm, ys, ym;
330: double lb = -TAO_INFINITY;
331: double ub = TAO_INFINITY;
332: double **xl, **xu;
333: double xi, xi1, xi2;
334: double cx, cy, radx, rady, height = 1.0;
336: info = DAGetInfo(da,PETSC_NULL,&mx,&my,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
337: PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(info);
338: user->hx = (user->r - user->l) / (mx - 1);
339: user->hy = (user->t - user->b) / (my - 1);
340: user->area = 0.5 * user->hx * user->hy;
342: info = DAVecGetArray(da, XL, (void**)&xl); CHKERRQ(info);
343: info = DAVecGetArray(da, XU, (void**)&xu); CHKERRQ(info);
344: info = DAGetCorners(da, &xs, &ys, TAO_NULL, &xm, &ym, TAO_NULL); CHKERRQ(info);
346: /* Compute default bounds */
347: for (j = ys; j < ys+ym; j++){
348: for (i = xs; i < xs+xm; i++){
350: if (j == 0 || j == my - 1) {
351: xi = user->l + i * user->hx;
352: xl[j][i] = xu[j][i] = -4 * (xi - user->l) * (xi - user->r);
353: } else if (i == 0 || i == mx - 1) {
354: xl[j][i] = xu[j][i] = 0.0;
355: } else {
356: xl[j][i] = lb;
357: xu[j][i] = ub;
358: }
360: }
361: }
363: /* Adjust lower bound if obstacle is present */
364: info = PetscOptionsHasName(PETSC_NULL, "-obst", &flg); CHKERRQ(info);
365: if (flg) {
366: radx = (user->r - user->l) * 0.25;
367: cx = user->l + 2.0 * radx;
368: rady = (user->t - user->b) * 0.25;
369: cy = user->b + 2.0 * rady;
370: for (j = ys; j < ys+ym; j++){
371: for (i = xs; i < xs+xm; i++){
372: xi1 = user->l + i * user->hx;
373: xi2 = user->b + j * user->hy;
374: if ( fabs(xi1 - cx) <= radx && fabs(xi2 - cy) <= rady ) {
375: xl[j][i] = height;
376: }
377: }
378: }
379: info = PetscLogFlops(8 + xm * ym * 6); CHKERRQ(info);
380: }
382: info = DAVecRestoreArray(da, XL, (void**)&xl); CHKERRQ(info);
383: info = DAVecRestoreArray(da, XU, (void**)&xu); CHKERRQ(info);
385: info = PetscLogFlops(12 * ym + 6); CHKERRQ(info);
387: return 0;
389: } /* DAGetBounds2d */
394: static int MINPACK_Bounds(TAO_APPLICATION tao, DA da, Vec XL,Vec XU, void *user1){
396: AppCtx *user=(AppCtx*)user1;
397: int i,j,k,limit=0,info,maxits=5;
398: int xs,ys,xm,ym;
399: int mx, my, bmy, bmx;
400: int row=0, bsize=0, lsize=0, tsize=0, rsize=0;
401: double bheight=user->bheight;
402: double one=1.0, two=2.0, three=3.0, tol=1e-10;
403: double fnorm,det,hx,hy,xt=0,yt=0;
404: double u1,u2,nf1,nf2,njac11,njac12,njac21,njac22;
405: double b=-0.5, t=0.5, l=-0.5, r=0.5;
406: PetscScalar scl,*xl,*xu, lb=TAO_NINFINITY, ub=TAO_INFINITY;
407: PetscTruth flg;
408: double **xll;
410: info = DAGetInfo(da,PETSC_NULL,&mx,&my,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
411: PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(info);
413: bheight=user->bheight,lb=-1000; ub=1000;
414: bmx=user->bmx; bmy=user->bmy;
416: info = DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ(info);
418: info = VecSet(XL, lb); CHKERRQ(info);
419: info = VecSet(XU, ub); CHKERRQ(info);
421: /*
422: Get pointers to vector data
423: */
424: info = DAVecGetArray(da,XL,(void**)&xll); CHKERRQ(info);
426: if (ys==0) bsize=xm;
427: if (xs==0) lsize=ym;
428: if (xs+xm==mx) rsize=ym;
429: if (ys+ym==my) tsize=xm;
430:
431: hx= (r-l)/(mx-1); hy=(t-b)/(my-1);
433: /* Compute the optional lower box */
434: for (i=xs; i< xs+xm; i++){
435: for (j=ys; j<ys+ym; j++){
436:
437: if (i>=(mx-bmx)/2 && i<mx-(mx-bmx)/2 && j>=(my-bmy)/2 && j<my-(my-bmy)/2 ){
438: xll[j][i] = bheight;
439: }
441: }
442: }
444: info = DAVecRestoreArray(da,XL,(void**)&xll); CHKERRQ(info);
446: /* Boundary Values */
447: info = VecGetArray(XL,&xl); CHKERRQ(info);
448: info = VecGetArray(XU,&xu); CHKERRQ(info);
450: for (j=0; j<4; j++){
451: if (j==0){
452: yt=b;
453: xt=l+hx*xs;
454: limit=bsize;
455: scl=1.0;
456: info = PetscOptionsGetReal(PETSC_NULL,"-bottom",&scl,&flg);
457: } else if (j==1){
458: yt=t;
459: xt=l+hx*xs;
460: limit=tsize;
461: scl=1.0;
462: info = PetscOptionsGetReal(PETSC_NULL,"-top",&scl,&flg);
463: } else if (j==2){
464: yt=b+hy*ys;
465: xt=l;
466: limit=lsize;
467: scl=1.0;
468: info = PetscOptionsGetReal(PETSC_NULL,"-left",&scl,&flg);
469: } else if (j==3){
470: yt=b+hy*ys;
471: xt=r;
472: limit=rsize;
473: scl=1.0;
474: info = PetscOptionsGetReal(PETSC_NULL,"-right",&scl,&flg);
475: }
477: for (i=0; i<limit; i++){
478: u1=xt;
479: u2=-yt;
480: for (k=0; k<maxits; k++){
481: nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
482: nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
483: fnorm=sqrt(nf1*nf1+nf2*nf2);
484: if (fnorm <= tol) break;
485: njac11=one+u2*u2-u1*u1;
486: njac12=two*u1*u2;
487: njac21=-two*u1*u2;
488: njac22=-one - u1*u1 + u2*u2;
489: det = njac11*njac22-njac21*njac12;
490: u1 = u1-(njac22*nf1-njac12*nf2)/det;
491: u2 = u2-(njac11*nf2-njac21*nf1)/det;
492: }
494: if (j==0){
495: row = i;
496: } else if (j==1){
497: row = (ym-1)*xm+i;
498: } else if (j==2){
499: row = (xm*i);
500: } else if (j==3){
501: row = xm*(i+1)-1;
502: }
503:
504: xl[row]=(u1*u1-u2*u2)*scl;
505: xu[row]=(u1*u1-u2*u2)*scl;
507: if (j==0 || j==1) {
508: xt=xt+hx;
509: } else if (j==2 || j==3){
510: yt=yt+hy;
511: }
512:
513: }
514:
515: }
517: info = VecRestoreArray(XL,&xl); CHKERRQ(info);
518: info = VecRestoreArray(XU,&xu); CHKERRQ(info);
520: return 0;
521: }
523: /*------- USER-DEFINED: routine to evaluate the function and gradient
524: at a local (rectangular element) level -------*/
528: /*
529: MSurfLocalFunctionGradient - Evaluates function and gradient over the
530: local rectangular element
532: Input:
533: coor - vector with the indices of the position of current element
534: in the first, second and third directions
535: x - current point (values over the current rectangular element)
536: ptr - user-defined application context
538: Output:
539: f - value of the objective funtion at the local rectangular element
540: g - gradient of the local function
542: */
543: static int MSurfLocalFunctionGradient(int coor[2], double x[4], double *f, double g[4], void *ptr) {
545: AppCtx *user = (AppCtx*)ptr;
547: double hx, hy, area;
548: double dvdx, dvdy, flow, fup;
549: double d1,d2;
551: hx = user->hx;
552: hy = user->hy;
553: area = user->area;
555: /* lower triangle contribution */
556: dvdx = (x[0] - x[1]);
557: dvdy = (x[0] - x[2]);
558: g[1] = (dvdx * (hy/hx))/(-2);
559: g[2] = (dvdy * (hx/hy))/(-2);
560: dvdx /= hx;
561: dvdy /= hy;
562: flow = sqrt( 1 + dvdx * dvdx + dvdy * dvdy );
563: *f=flow;
564: g[1] /= flow;
565: g[2] /= flow;
566: g[0] = -(g[1]+g[2]);
568: /* upper triangle contribution */
569: dvdx = (x[3] - x[2]);
570: dvdy = (x[3] - x[1]);
571: d1 = (dvdy*(hx/hy))/(-2);
572: d2 = (dvdx*(hy/hx))/(-2);
573: dvdx /= hx;
574: dvdy /= hy;
575: fup = sqrt( 1 + dvdx * dvdx + dvdy * dvdy );
576: *f += fup;
577: g[1] += d1/fup;
578: g[2] += d2/fup;
579: g[3] = -(d1+d2)/fup;
581: *f *= area;
583: return 0;
584: } /* MSurfLocalFunctionGradient */
590: /*
591: MSurfLocalHessian - Computes the Hessian of the local (partial) function
592: defined over the current rectangle
594: Input:
595: coor - vector with the indices of the position of current element
596: in the first, second and third directions
597: x - current local solution (over the rectangle only)
598: ptr - user-defined application context
600: Output:
601: H - Hessian matrix of the local function (wrt the four
602: points of the rectangle only)
604: */
605: static int MSurfLocalHessian(int coor[2], double x[4], double H[4][4],void *ptr) {
607: AppCtx *user = (AppCtx*)ptr;
608: double hx, hy, area, byhxhx, byhyhy;
609: double dvdx, dvdy, flow, fup;
610: double areadivf, areadivf3;
612: hx = user->hx;
613: hy = user->hy;
614: area = user->area;
615:
616: byhxhx = 1.0 / (hx * hx);
617: byhyhy = 1.0 / (hy * hy);
619: /* 0 is 0,0; 1 is 1,0; 2 is 0,1; 3 is 1,1 */
620: dvdx = (x[0] - x[1]) / hx; /* lower triangle contribution */
621: dvdy = (x[0] - x[2]) / hy;
622: flow = sqrt( 1 + dvdx * dvdx + dvdy * dvdy );
623: dvdx = dvdx / hx;
624: dvdy = dvdy / hy;
625: areadivf = area / flow;
626: areadivf3 = areadivf / (flow * flow);
627: H[0][0] = areadivf * (byhxhx + byhyhy) - areadivf3 * (dvdx + dvdy) * (dvdx + dvdy);
628: H[0][1] = areadivf * (-byhxhx) + areadivf3 * (dvdx + dvdy) * (dvdx);
629: H[0][2] = areadivf * (-byhyhy) + areadivf3 * (dvdx + dvdy) * (dvdy);
630: H[0][3] = 0.0;
631: H[1][1] = areadivf * byhxhx - areadivf3 * dvdx * dvdx;
632: H[1][2] = areadivf3 * (-dvdx) * dvdy;
633: H[2][2] = areadivf * byhyhy - areadivf3 * dvdy * dvdy;
635: /* upper triangle contribution */
636: dvdx = (x[3] - x[2]) / hx;
637: dvdy = (x[3] - x[1]) / hy;
638: fup = sqrt( 1 + dvdx * dvdx + dvdy * dvdy );
639: dvdx = dvdx / hx;
640: dvdy = dvdy / hy;
641: areadivf = area / fup;
642: areadivf3 = areadivf / (fup * fup);
643: H[1][1] += areadivf * byhyhy - areadivf3 * dvdy * dvdy;
644: H[1][2] += areadivf3 * (-dvdy) * dvdx;
645: H[2][2] += areadivf * byhxhx - areadivf3 * (dvdx * dvdx);
646: H[1][3] = areadivf * (-byhyhy) + areadivf3 * (dvdx + dvdy) * dvdy;
647: H[2][3] = areadivf * (-byhxhx) + areadivf3 * (dvdx + dvdy) * dvdx;
648: H[3][3] = areadivf * (byhxhx + byhyhy) - areadivf3 * (dvdx + dvdy) * (dvdx + dvdy);
650: H[1][0] = H[0][1];
651: H[2][0] = H[0][2];
652: H[3][0] = H[0][3];
653: H[2][1] = H[1][2];
654: H[3][1] = H[1][3];
655: H[3][2] = H[2][3];
657: return 0;
659: } /* MSurfLocalHessian */
662: /*------- USER-DEFINED: routine to evaluate the function
663: and gradient at the whole grid -------*/
667: /*
668: WholeMSurfFunctionGradient - Evaluates function and gradient over the
669: whole grid
671: Input:
672: daapplication - TAO application object
673: da - distributed array
674: X - the current point, at which the function and gradient are evaluated
675: ptr - user-defined application context
677: Output:
678: f - value of the objective funtion at X
679: G - gradient at X
680: */
681: static int WholeMSurfFunctionGradient(TAO_APPLICATION daapplication, DA da, Vec X, double *f, Vec G, void *ptr) {
683: AppCtx *user = (AppCtx*)ptr;
684: Vec localX, localG;
685: int info, i, j;
686: int xs, xm, gxs, gxm, xe, ys, ym, gys, gym, ye;
687: double **x, **g;
688: double floc = 0.0;
689: PetscScalar zero = 0.0;
691: double hx, hy, area;
692: double dvdx, dvdy, flow, fup;
693: double areadivf;
695: hx = user->hx;
696: hy = user->hy;
697: area = user->area;
699: info = DAGetLocalVector(da, &localX); CHKERRQ(info);
700: info = DAGetLocalVector(da, &localG); CHKERRQ(info);
701: info = VecSet(G, zero); CHKERRQ(info);
702: info = VecSet(localG, zero); CHKERRQ(info);
704: info = DAGlobalToLocalBegin(da, X, INSERT_VALUES, localX); CHKERRQ(info);
705: info = DAGlobalToLocalEnd(da, X, INSERT_VALUES, localX); CHKERRQ(info);
707: info = DAVecGetArray(da, localX, (void**)&x); CHKERRQ(info);
708: info = DAVecGetArray(da, localG, (void**)&g); CHKERRQ(info);
710: info = DAGetCorners(da, &xs, &ys, TAO_NULL, &xm, &ym, TAO_NULL); CHKERRQ(info);
711: info = DAGetGhostCorners(da, &gxs, &gys, TAO_NULL, &gxm, &gym, TAO_NULL); CHKERRQ(info);
713: xe = gxs + gxm - 1;
714: ye = gys + gym - 1;
715: for (j = ys; j < ye; j++) {
716: for (i = xs; i < xe; i++) {
718: /* lower triangle contribution */
719: dvdx = (x[j][i] - x[j][i+1]) / hx;
720: dvdy = (x[j][i] - x[j+1][i]) / hy;
721: flow = sqrt( 1 + dvdx * dvdx + dvdy * dvdy );
722: areadivf = area / flow;
723: g[j][i] += (dvdx / hx + dvdy / hy) * areadivf;
724: g[j][i+1] += (-dvdx / hx) * areadivf;
725: g[j+1][i] += (-dvdy / hy) * areadivf;
727: /* upper triangle contribution */
728: dvdx = (x[j+1][i+1] - x[j+1][i]) / hx;
729: dvdy = (x[j+1][i+1] - x[j][i+1]) / hy;
730: fup = sqrt( 1 + dvdx * dvdx + dvdy * dvdy );
731: areadivf = area / fup;
732: g[j][i+1] += (-dvdy / hy) * areadivf;
733: g[j+1][i] += (-dvdx / hx) * areadivf;
734: g[j+1][i+1] += (dvdx / hx + dvdy / hy) * areadivf;
736: floc += area * (flow + fup);
738: }
739: }
741: info = MPI_Allreduce(&floc, f, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); CHKERRQ(info);
743: info = DAVecRestoreArray(da, localX, (void**)&x); CHKERRQ(info);
744: info = DAVecRestoreArray(da, localG, (void**)&g); CHKERRQ(info);
746: info = DALocalToGlobalBegin(da, localG, G); CHKERRQ(info);
747: info = DALocalToGlobalEnd(da, localG, G); CHKERRQ(info);
749: info = DARestoreLocalVector(da, &localX); CHKERRQ(info);
750: info = DARestoreLocalVector(da, &localG); CHKERRQ(info);
752: info = PetscLogFlops((xe-xs) * (ye-ys) * 42); CHKERRQ(info);
754: return 0;
755: } /* WholeMSurfFunctionGradient */
758: /*------- USER-DEFINED: routine to evaluate the Hessian
759: at the whole grid -------*/
763: /*
764: WholeMSurfHessian - Evaluates Hessian over the whole grid
766: Input:
767: daapplication - TAO application object
768: da - distributed array
769: X - the current point, at which the function and gradient are evaluated
770: ptr - user-defined application context
772: Output:
773: H - Hessian at X
774: */
775: static int WholeMSurfHessian(TAO_APPLICATION daapplication, DA da, Vec X, Mat H, void *ptr) {
777: AppCtx *user = (AppCtx*)ptr;
778: Vec localX;
779: int info, i, j, ind[4];
780: int xs, xm, gxs, gxm, xe, ys, ym, gys, gym, ye;
781: double smallH[4][4];
782: double **x;
784: double hx, hy, area, byhxhx, byhyhy;
785: double dvdx, dvdy, flow, fup;
786: double areadivf, areadivf3;
787: PetscTruth assembled;
789: hx = user->hx;
790: hy = user->hy;
791: area = user->area;
792:
793: byhxhx = 1.0 / (hx * hx);
794: byhyhy = 1.0 / (hy * hy);
796: info = DAGetLocalVector(da, &localX); CHKERRQ(info);
797: info = MatAssembled(H,&assembled); CHKERRQ(info);
798: if (assembled){info = MatZeroEntries(H); CHKERRQ(info);}
800: info = DAGlobalToLocalBegin(da, X, INSERT_VALUES, localX); CHKERRQ(info);
801: info = DAGlobalToLocalEnd(da, X, INSERT_VALUES, localX); CHKERRQ(info);
803: info = DAVecGetArray(da, localX, (void**)&x); CHKERRQ(info);
805: info = DAGetCorners(da, &xs, &ys, TAO_NULL, &xm, &ym, TAO_NULL); CHKERRQ(info);
806: info = DAGetGhostCorners(da, &gxs, &gys, TAO_NULL, &gxm, &gym, TAO_NULL); CHKERRQ(info);
808: xe = gxs + gxm - 1;
809: ye = gys + gym - 1;
810: for (j = ys; j < ye; j++) {
811: for (i = xs; i < xe; i++) {
813: /* 0 is 0,0; 1 is 1,0; 2 is 0,1; 3 is 1,1 */
814: dvdx = (x[j][i] - x[j][i+1]) / hx; /* lower triangle contribution */
815: dvdy = (x[j][i] - x[j+1][i]) / hy;
816: flow = sqrt( 1 + dvdx * dvdx + dvdy * dvdy );
817: dvdx = dvdx / hx;
818: dvdy = dvdy / hy;
819: areadivf = area / flow;
820: areadivf3 = areadivf / (flow * flow);
821: smallH[0][0] = areadivf * (byhxhx + byhyhy) - areadivf3 * (dvdx + dvdy) * (dvdx + dvdy);
822: smallH[0][1] = areadivf * (-byhxhx) + areadivf3 * (dvdx + dvdy) * (dvdx);
823: smallH[0][2] = areadivf * (-byhyhy) + areadivf3 * (dvdx + dvdy) * (dvdy);
824: smallH[0][3] = 0.0;
825: smallH[1][1] = areadivf * byhxhx - areadivf3 * dvdx * dvdx;
826: smallH[1][2] = areadivf3 * (-dvdx) * dvdy;
827: smallH[2][2] = areadivf * byhyhy - areadivf3 * dvdy * dvdy;
829: /* upper triangle contribution */
830: dvdx = (x[j+1][i+1] - x[j+1][i]) / hx;
831: dvdy = (x[j+1][i+1] - x[j][i+1]) / hy;
832: fup = sqrt( 1 + dvdx * dvdx + dvdy * dvdy );
833: dvdx = dvdx / hx;
834: dvdy = dvdy / hy;
835: areadivf = area / fup;
836: areadivf3 = areadivf / (fup * fup);
837: smallH[1][1] += areadivf * byhyhy - areadivf3 * dvdy * dvdy;
838: smallH[1][2] += areadivf3 * (-dvdy) * dvdx;
839: smallH[2][2] += areadivf * byhxhx - areadivf3 * (dvdx * dvdx);
840: smallH[1][3] = areadivf * (-byhyhy) + areadivf3 * (dvdx + dvdy) * dvdy;
841: smallH[2][3] = areadivf * (-byhxhx) + areadivf3 * (dvdx + dvdy) * dvdx;
842: smallH[3][3] = areadivf * (byhxhx + byhyhy) - areadivf3 * (dvdx + dvdy) * (dvdx + dvdy);
844: smallH[1][0] = smallH[0][1];
845: smallH[2][0] = smallH[0][2];
846: smallH[3][0] = smallH[0][3];
847: smallH[2][1] = smallH[1][2];
848: smallH[3][1] = smallH[1][3];
849: smallH[3][2] = smallH[2][3];
851: ind[0] = (j-gys) * gxm + (i-gxs);
852: ind[1] = ind[0] + 1;
853: ind[2] = ind[0] + gxm;
854: ind[3] = ind[2] + 1;
855: info = MatSetValuesLocal(H,4,ind,4,ind,(PetscScalar*)smallH,ADD_VALUES); CHKERRQ(info);
857: }
858: }
860: info = DAVecRestoreArray(da, localX, (void**)&x); CHKERRQ(info);
862: info = MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY); CHKERRQ(info);
863: info = MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY); CHKERRQ(info);
864: info = MatSetOption(H, MAT_SYMMETRIC); CHKERRQ(info);
866: info = DARestoreLocalVector(da, &localX); CHKERRQ(info);
868: info = PetscLogFlops((xe-xs) * (ye-ys) * 83 + 4); CHKERRQ(info);
869: return 0;
871: } /* WholeMSurfHessian */