Actual source code: jbearing3.c

  1: /*$Id: jbearing.c, v 1.1 2002/08/08 11:35 lopezca@mauddib.mcs.anl.gov $*/

  3: /* Program usage: mpirun -np <proc> jbearing [-help] [all TAO options] */

  5: /*
  6:   Include "tao.h" so we can use TAO solvers.
  7:   petscda.h for distributed array
  8:   ad_deriv.h for AD gradient
  9: */

 11: #include "petscda.h"
 12: #include "tao.h"
 13: #include "taodaapplication.h"

 15: static char help[] ="Pressure distribution in a Journal Bearing. \n\
 16: This example is based on the problem DPJB from the MINPACK-2 test suite.\n\
 17: This pressure journal bearing problem is an example of elliptic variational\n\
 18: problem defined over a two dimensional rectangle. By discretizing the domain \n\
 19: into triangular elements, the pressure surrounding the journal bearing is\n\
 20: defined as the minimum of a quadratic function whose variables are bounded\n\
 21: below by zero. The command line options are:\n\
 22:   -ecc <ecc>, where <ecc> = epsilon parameter\n\
 23:   -b <b>, where <b> = half the upper limit in the 2nd coordinate direction\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\n";

 30: /*T
 31:    Concepts: TAO - Solving a bounded minimization problem
 32:    Routines: TaoInitialize(); TaoFinalize();
 33:    Routines: TaoCreate(); TaoDestroy();
 34:    Routines: DAApplicationCreate(); DAApplicationDestroy();
 35:    Routines: DAAppSetVariableBoundsRoutine;
 36:    Routines: DAAppSetElementObjectiveAndGradientRoutine();
 37:    Routines: DAAppSetElementHessianRoutine();
 38:    Routines: DAAppSetObjectiveAndGradientRoutine();
 39:    Routines: DAAppSetADElementFunctionGradient();
 40:    Routines: DAAppSetHessianRoutine();
 41:    Routines: TaoSetOptions();
 42:    Routines: TaoGetSolutionStatus(); TaoDAAppSolve();
 43:    Routines: DAAppSetBeforeMonitor(); DAAppSetAfterMonitor
 44:    Routines: DAAppGetSolution(); TaoView();
 45:    Routines: DAAppGetInterpolationMatrix();
 46:    Processors: n
 47: T*/
 48:  
 49: /*
 50:    User-defined application context - contains data needed by the
 51:    application-provided call-back routines.
 52: */

 54: int  ad_JBearLocalFunction(PetscInt[2] ,DERIV_TYPE[4], DERIV_TYPE *, void*);
 55: typedef struct {

 57:   InactiveDouble      *wq, *wl;      /* vectors with the parameters w_q(x) and w_l(x) */
 58:   InactiveDouble      hx, hy;        /* increment size in both directions */
 59:   InactiveDouble      area;          /* area of the triangles */

 61: } ADFGCtx;


 64: typedef struct {

 66:   PetscReal      ecc;           /* epsilon value */
 67:   PetscReal      b;             /* 0.5 * upper limit for 2nd variable */
 68:   double      *wq, *wl;      /* vectors with the parameters w_q(x) and w_l(x) */
 69:   double      hx, hy;        /* increment size in both directions */
 70:   double      area;          /* area of the triangles */

 72:   PetscInt    mx, my;        /* discretization including boundaries */

 74:   ADFGCtx     fgctx;         /* Used only when an ADIC generated gradient is used */

 76: } AppCtx;

 78: /* User-defined routines found in this file */
 79: static int AppCtxInitialize(void *ptr);
 80: static int FormInitialGuess(DA, Vec);

 82: static int JBearLocalFunctionGradient(PetscInt[2], double x[4], double *f, double g[4], void *ptr);
 83: static int JBearLocalHessian(PetscInt[2], double x[4], double H[4][4], void *ptr);

 85: static int WholeJBearFunctionGradient(TAO_APPLICATION,DA,Vec,double *,Vec,void*);
 86: static int WholeJBearHessian(TAO_APPLICATION,DA,Vec,Mat,void*);

 88: static int DASetBounds(TAO_APPLICATION, DA, Vec, Vec, void*);

 90: static int MyGridMonitorBefore(TAO_APPLICATION, DA, PetscInt, void *);
 91: static int MyGridMonitorAfter1(TAO_APPLICATION, DA, PetscInt, void *);

 95: int main( int argc, char **argv ) {

 97:   PetscInt             info,iter;             /* used to check for functions returning nonzeros */
 98:   PetscInt             nlevels;                                             /* multigrid levels */
 99:   PetscInt             Nx,Ny;
100:   double          ff,gnorm;
101:   DA              DAarray[20];
102:   Vec             X;
103:   PetscTruth      flg, PreLoad = PETSC_TRUE;                                      /* flags */
104:   AppCtx          user;                                       /* user-defined work context */
105:   TaoMethod       method = "tao_gpcg";                              /* minimization method */
106:   TAO_SOLVER      tao;                                        /* TAO_SOLVER solver context */
107:   TAO_APPLICATION JBearApp;                                    /* The PETSc application */
108:   TaoTerminateReason reason;
109:   KSP ksp;
110:   PC  pc;

112:   /* Initialize TAO */
113:   PetscInitialize(&argc, &argv, (char *)0, help);
114:   TaoInitialize(&argc, &argv, (char *)0, help);

116:   PreLoadBegin(PreLoad,"Solve");
117:   
118:   info = AppCtxInitialize((void*)&user); CHKERRQ(info);
119:   
120:   nlevels=5;
121:   info = PetscOptionsGetInt(PETSC_NULL,"-nlevels",&nlevels,&flg); CHKERRQ(info);
122:   if (PreLoadIt == 0) {
123:     nlevels = 1; user.mx = 11; user.my = 11; }

125:   PetscPrintf(MPI_COMM_WORLD,"\n---- Journal Bearing Problem -----\n\n");

127:   /* Let PETSc determine the vector distribution */
128:   Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;

130:   /* Create distributed array (DA) to manage parallel grid and vectors  */
131:   info = DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_BOX,user.mx,
132:                     user.my,Nx,Ny,1,1,PETSC_NULL,PETSC_NULL,&DAarray[0]); CHKERRQ(info);
133:   for (iter=1;iter<nlevels;iter++){
134:     info = DARefine(DAarray[iter-1],PETSC_COMM_WORLD,&DAarray[iter]); CHKERRQ(info);
135:   }

137:   /* Create TAO solver and set desired solution method */
138:   info = TaoCreate(MPI_COMM_WORLD,method,&tao); CHKERRQ(info);
139:   info = TaoApplicationCreate(PETSC_COMM_WORLD, &JBearApp); CHKERRQ(info);
140:   info = TaoAppSetDAApp(JBearApp,DAarray,nlevels); CHKERRQ(info);

142:   /* Sets routine bounds evaluation */
143:   info = DAAppSetVariableBoundsRoutine(JBearApp,DASetBounds,(void *)&user); CHKERRQ(info);

145:   info = PetscOptionsHasName(TAO_NULL, "-byelement", &flg); CHKERRQ(info);
146:   if (flg) {

148:     /* Sets routines for function and gradient evaluation, element by element */
149:     info = PetscOptionsHasName(TAO_NULL, "-adic", &flg); CHKERRQ(info);
150:     if (flg) {
151:       info = DAAppSetADElementFunctionGradient(JBearApp,ad_JBearLocalFunction,248,(void *)&user.fgctx); CHKERRQ(info);
152:     } else {
153:       info = DAAppSetElementObjectiveAndGradientRoutine(JBearApp,JBearLocalFunctionGradient,63,(void *)&user); CHKERRQ(info);
154:     }
155:     /* Sets routines for Hessian evaluation, element by element */
156:     info = DAAppSetElementHessianRoutine(JBearApp,JBearLocalHessian,16,(void*)&user); CHKERRQ(info);

158:   } else {

160:     /* Sets routines for function and gradient evaluation, all in one routine */
161:     info = DAAppSetObjectiveAndGradientRoutine(JBearApp,WholeJBearFunctionGradient,(void *)&user); CHKERRQ(info);

163:     /* Sets routines for Hessian evaluation, all in one routine */
164:     info = DAAppSetHessianRoutine(JBearApp,WholeJBearHessian,(void*)&user); CHKERRQ(info);    
165:     
166:   }

168:   info = DAAppSetBeforeMonitor(JBearApp,MyGridMonitorBefore,(void*)&user); CHKERRQ(info);
169:   info = DAAppSetAfterMonitor(JBearApp,MyGridMonitorAfter1,(void*)&user); CHKERRQ(info);
170:   info = PetscOptionsHasName(TAO_NULL,"-tao_monitor", &flg); CHKERRQ(info);
171:   if (flg){
172:     info = DAAppPrintInterpolationError(JBearApp); CHKERRQ(info);
173:     info = DAAppPrintStageTimes(JBearApp); CHKERRQ(info);
174:   }
175:   /* Check for any tao command line options */
176:   info = TaoAppSetRelativeTolerance(JBearApp,1.0e-6); CHKERRQ(info);
177:   info = TaoSetTolerances(tao,0,0,0,0); CHKERRQ(info);
178:   info = TaoSetGradientTolerances(tao,0,0,0); CHKERRQ(info);

180:   info = TaoAppGetKSP(JBearApp,&ksp);CHKERRQ(info);
181:   info = KSPSetType(ksp,KSPCG); CHKERRQ(info);
182:   info = KSPSetTolerances(ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,10);CHKERRQ(info);
183:   info = KSPGetPC(ksp,&pc);CHKERRQ(info);
184:   info = PCSetType(pc,PCBJACOBI);CHKERRQ(info);

186:   info = TaoSetOptions(JBearApp,tao); CHKERRQ(info);

188:   info = DAAppGetSolution(JBearApp,0,&X); CHKERRQ(info);
189:   info = FormInitialGuess(DAarray[0],X); CHKERRQ(info);
190:   info = DAAppSetInitialSolution(JBearApp,X); CHKERRQ(info);
191:   /* SOLVE THE APPLICATION */
192:   info = TaoDAAppSolve(JBearApp,tao);  CHKERRQ(info);

194:   /* Get information on termination */
195:   info = TaoGetSolutionStatus(tao,&iter,&ff,&gnorm,0,0,&reason); CHKERRQ(info);
196:   if (reason <= 0 ){
197:     PetscPrintf(MPI_COMM_WORLD,"Try a different TAO method, adjust some parameters, or check the function evaluation routines\n");
198:     PetscPrintf(MPI_COMM_WORLD," Iterations: %d,  Function Value: %4.2e, Residual: %4.2e \n",iter,ff,gnorm);
199:   }

201:   info = PetscOptionsHasName(PETSC_NULL,"-view_sol",&flg); CHKERRQ(info);
202:   if (flg){
203:     info = DAAppGetSolution(JBearApp,nlevels-1,&X); CHKERRQ(info);
204:     info=VecView(X,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(info);
205:   }

207:   /*  To View TAO solver information */
208:   // info = TaoView(tao); CHKERRQ(info);

210:   /* Free TAO data structures */
211:   info = TaoDestroy(tao); CHKERRQ(info);
212:   info = TaoAppDestroy(JBearApp); CHKERRQ(info);

214:   /* Free PETSc data structures */
215:   for (iter=0;iter<nlevels;iter++){
216:     info = DADestroy(DAarray[iter]); CHKERRQ(info);
217:   }

219:   PreLoadEnd();

221:   /* Finalize TAO */
222:   TaoFinalize();
223:   PetscFinalize();

225:   return 0;
226: } /* main */



230: /*----- The following two routines
231:    
232:   MyGridMonitorBefore    MyGridMonitorAfter

234:   help diplay info of iterations at every grid level -------*/

238: static int MyGridMonitorBefore(TAO_APPLICATION myapp, DA da, PetscInt level, void *ctx) {

240:   AppCtx *user = (AppCtx*)ctx;
241:   int info;
242:   PetscInt mx,my;
243:   double t;

245:   info = DAGetInfo(da,PETSC_NULL,&mx,&my,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
246:                    PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(info);
247:   user->mx = mx;
248:   user->my = my;

250:   user->hx = (2.0 * 3.14159265358979) / (user->mx - 1);
251:   user->hy = (2.0 * user->b) / (user->my - 1);
252:   user->area = 0.5 * user->hx * user->hy;

254:   user->wq = new double[user->mx];
255:   user->wl = new double[user->mx];
256:   for (PetscInt i=0; i<user->mx; i++) {
257:     t = 1.0 + user->ecc * cos(i*user->hx);
258:     user->wq[i] = t*t*t;
259:     user->wl[i] = user->ecc * sin(i*user->hx);
260:   }
261:   info = PetscLogFlops(8 + 7 * user->mx); CHKERRQ(info);

263:   user->fgctx.hx   = user->hx;
264:   user->fgctx.hy   = user->hy;
265:   user->fgctx.area = user->area;
266:   user->fgctx.wq = user->wq;
267:   user->fgctx.wl = user->wl;

269:   PetscPrintf(MPI_COMM_WORLD,"Grid: %d,    mx: %d     my: %d   \n",level,mx,my);

271:   return 0;
272: }

276: static int MyGridMonitorAfter1(TAO_APPLICATION myapp, DA da, PetscInt level, void *ctx){
277:   
278:   AppCtx *user = (AppCtx*)ctx;
279:   delete [] user->wq; delete [] user->wl;
280:   return 0;
281: }



285: /*------- USER-DEFINED: initialize the application context information -------*/

289: /*
290:   AppCtxInitialize - Sets initial values for the application context parameters

292:   Input:
293:     ptr - void user-defined application context

295:   Output:
296:     ptr - user-defined application context with the default or user-provided
297:              parameters
298: */
299: static int AppCtxInitialize(void *ptr) {

301:   AppCtx *user = (AppCtx*)ptr;
302:   PetscTruth    flg;            /* flag for PETSc calls */
303:   int info;

305:   /* Specify default parameters */
306:   user->ecc = 0.1;
307:   user->b = 10.0;
308:   user->mx = user->my = 11;

310:   /* Check for command line arguments that override defaults */
311:   info = PetscOptionsGetReal(TAO_NULL, "-ecc", &user->ecc, &flg); CHKERRQ(info);
312:   info = PetscOptionsGetReal(TAO_NULL, "-b", &user->b, &flg); CHKERRQ(info);
313:   info = PetscOptionsGetInt(TAO_NULL, "-mx", &user->mx, &flg); CHKERRQ(info);
314:   info = PetscOptionsGetInt(TAO_NULL, "-my", &user->my, &flg); CHKERRQ(info);

316:   return 0;
317: } /* AppCtxInitialize */


322: static int FormInitialGuess(DA da, Vec X)
323: {
324:   int info;
325:   PetscInt    i, j, mx;
326:   PetscInt    xs, ys, xm, ym, xe, ye;
327:   PetscReal hx, val;
328:   double **x;

330:   /* Get local mesh boundaries */
331:   info = DAGetInfo(da,PETSC_NULL,&mx,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
332:                    PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(info);
333:   hx = 2.0*4.0*atan(1.0)/(mx-1);

335:   info = DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ(info);
336:   xe = xs+xm; ye = ys+ym;

338:   info = DAVecGetArray(da, X, (void**)&x); CHKERRQ(info);
339:   /* Compute initial guess over locally owned part of mesh */
340:   for (j=ys; j<ye; j++) {  /*  for (j=0; j<my; j++) */
341:     for (i=xs; i<xe; i++) {  /*  for (i=0; i<mx; i++) */
342:       val = PetscMax(sin((i+1.0)*hx),0.0);
343:       x[j][i] = val;
344:       x[j][i] = 0;
345:     }
346:   }
347:   info = DAVecRestoreArray(da, X, (void**)&x); CHKERRQ(info);

349:   return 0;
350: }


353: /*------- USER-DEFINED: set the upper and lower bounds for the variables  -------*/
356: /*
357:   FormBounds - Forms bounds on the variables

359:   Input:
360:     user - user-defined application context

362:   Output:
363:     XL - vector of lower bounds
364:     XU - vector of upper bounds

366: */
367: static int DASetBounds(TAO_APPLICATION daapplication, DA da, Vec XL, Vec XU, void *ptr) {

369:   AppCtx *user = (AppCtx*)ptr;
370:   int info;
371:   PetscInt i, j, mx, my;
372:   PetscInt xs, xm, ys, ym;
373:   double **xl, **xu;

375:   mx = user->mx;
376:   my = user->my;

378:   info = DAVecGetArray(da, XL, (void**)&xl); CHKERRQ(info);
379:   info = DAVecGetArray(da, XU, (void**)&xu); CHKERRQ(info);
380:   info = DAGetCorners(da, &xs, &ys, TAO_NULL, &xm, &ym, TAO_NULL); CHKERRQ(info);

382:   for (j = ys; j < ys+ym; j++){
383:     for (i = xs; i < xs+xm; i++){
384:       xl[j][i] = 0.0;
385:       if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) {
386:         xu[j][i] = 0.0;
387:       } else {
388:         xu[j][i] = TAO_INFINITY;
389:       }
390:     }
391:   }

393:   info = DAVecRestoreArray(da, XL, (void**)&xl); CHKERRQ(info);
394:   info = DAVecRestoreArray(da, XU, (void**)&xu); CHKERRQ(info);

396:   return 0;

398: } /* DASetBounds */



404: /*
405:   JBearLocalFunctionGradient - Evaluates function and gradient over the 
406:       local rectangular element

408:   Input:
409:     coor - vector with the indices of the position of current element
410:              in the first, second and third directions
411:     x - current point (values over the current rectangular element)
412:     ptr - user-defined application context

414:   Output:
415:     f - value of the objective funtion at the local rectangular element
416:     g - gradient of the local function

418: */
419: static int JBearLocalFunctionGradient(PetscInt coor[2], double x[4], double *f, double g[4], void *ptr) {

421:   AppCtx *user = (AppCtx*)ptr;

423:   double avgWq, sqGrad, avgWV, fl, fu;
424:   double hx, hy, area, aread3, *wq, *wl;
425:   double dvdx, dvdy;
426:   PetscInt i;

428:   hx = user->hx;
429:   hy = user->hy;
430:   area = user->area;
431:   aread3 = area / 3.0;
432:   wq = user->wq;
433:   wl = user->wl;
434:   i = coor[0];

436:   /* lower triangle contribution */
437:   dvdx = (x[0] - x[1]) / hx;
438:   dvdy = (x[0] - x[2]) / hy;
439:   sqGrad = dvdx * dvdx + dvdy * dvdy;
440:   avgWq = (2.0 * wq[i] + wq[i+1]) / 6.0;
441:   avgWV = (wl[i]*x[0] + wl[i+1]*x[1] + wl[i]*x[2]) / 3.0;
442:   fl = avgWq * sqGrad - avgWV;

444:   dvdx = dvdx * hy * avgWq;
445:   dvdy = dvdy * hx * avgWq;
446:   g[0] = ( dvdx + dvdy ) - wl[i] * aread3;
447:   g[1] = ( -dvdx ) - wl[i+1] * aread3;
448:   g[2] = ( -dvdy ) - wl[i] * aread3;

450:   /* upper triangle contribution */
451:   dvdx = (x[3] - x[2]) / hx; 
452:   dvdy = (x[3] - x[1]) / hy;
453:   sqGrad = dvdx * dvdx + dvdy * dvdy;
454:   avgWq = (2.0 * wq[i+1] + wq[i]) / 6.0;
455:   avgWV = (wl[i+1]*x[1] + wl[i]*x[2] + wl[i+1]*x[3]) / 3.0;
456:   fu = avgWq * sqGrad - avgWV;

458:   dvdx = dvdx * hy * avgWq;
459:   dvdy = dvdy * hx * avgWq;
460:   g[1] += (-dvdy) - wl[i+1] * aread3;
461:   g[2] +=  (-dvdx) - wl[i] * aread3;
462:   g[3] = ( dvdx + dvdy ) - wl[i+1] * aread3;

464:   *f = area * (fl + fu);

466:   return 0;
467: } /* JBearLocalFunctionGradient */


470: /*------- USER-DEFINED: routine to evaluate the Hessian
471:            at a local (rectangular element) level       -------*/
474: /*
475:   JBearLocalHessian - Computes the Hessian of the local (partial) function
476:          defined over the current rectangle

478:   Input:
479:     coor - vector with the indices of the position of current element
480:              in the first, second and third directions
481:     x - current local solution (over the rectangle only)
482:     ptr - user-defined application context

484:   Output:
485:     H - Hessian matrix of the local function (wrt the four
486:            points of the rectangle only)

488: */
489: static int JBearLocalHessian(PetscInt coor[2], double x[4], double H[4][4],void *ptr) {

491:   AppCtx *user = (AppCtx*)ptr;
492:   double wql, wqu;
493:   double hx, hy, dydx, dxdy;
494:   double *wq;
495:   double wqldydx,wqldxdy,wqudydx,wqudxdy;
496:   PetscInt i;

498:   hx = user->hx;
499:   hy = user->hy;
500:   wq = user->wq;
501:   i = coor[0];

503:   dxdy = hx / hy;
504:   dydx = hy / hx;
505:   wql = (2.0 * wq[i] + wq[i+1]) / 6.0;
506:   wqu = (wq[i] + 2.0 * wq[i+1]) / 6.0;
507:   wqldydx = wql * dydx;
508:   wqldxdy = wql * dxdy;
509:   wqudydx = wqu * dydx;
510:   wqudxdy = wqu * dxdy;

512:           /* Hessian contribution at 0,0 */
513:   H[0][0] = wqldxdy + wqldydx;
514:   H[0][1] = H[1][0] = -wqldydx;
515:   H[0][2] = H[2][0] = -wqldxdy;
516:   H[0][3] = H[3][0] = 0.0;

518:           /* Hessian contribution at 1,0 */
519:   H[1][1] = wqldydx + wqudxdy;
520:   H[1][2] = H[2][1] = 0.0;
521:   H[1][3] = H[3][1] = -wqudxdy; 

523:           /* Hessian contribution at 0,1 */
524:   H[2][2] = wqldxdy + wqudydx;
525:   H[2][3] = H[3][2] = -wqudydx; 

527:           /* Hessian contribution at 1,1 */
528:   H[3][3] = wqudydx + wqudxdy;

530:   return 0;

532: } /* JBearLocalHessian */


535: /*------- USER-DEFINED: routine to evaluate the function 
536:           and gradient at the whole grid             -------*/

540: /*
541:   WholeJBearFunctionGradient - Evaluates function and gradient over the 
542:       whole grid

544:   Input:
545:     daapplication - TAO application object
546:     da  - distributed array
547:     X   - the current point, at which the function and gradient are evaluated
548:     ptr - user-defined application context

550:   Output:
551:     f - value of the objective funtion at X
552:     G - gradient at X
553: */
554: static int WholeJBearFunctionGradient(TAO_APPLICATION daapplication, DA da, Vec X, double *f, Vec G, void *ptr) {

556:   AppCtx *user = (AppCtx*)ptr;
557:   Vec localX, localG;
558:   int info;
559:   PetscInt i, j;
560:   PetscInt xs, xm, gxs, gxm, xe, ys, ym, gys, gym, ye;
561:   double **x, **g;
562:   double floc = 0.0;
563:   PetscScalar zero = 0.0;

565:   double avgWq, sqGrad, avgWV, fl, fu;
566:   double hx, hy, area, aread3, *wq, *wl;
567:   double dvdx, dvdy;

569:   hx = user->hx;
570:   hy = user->hy;
571:   area = user->area;
572:   aread3= area/3.0;
573:   wq = user->wq;
574:   wl = user->wl;

576:   info = DAGetLocalVector(da, &localX); CHKERRQ(info);
577:   info = DAGetLocalVector(da, &localG); CHKERRQ(info);
578:   info = VecSet(G, zero); CHKERRQ(info);
579:   info = VecSet(localG, zero); CHKERRQ(info);

581:   info = DAGlobalToLocalBegin(da, X, INSERT_VALUES, localX); CHKERRQ(info);
582:   info = DAGlobalToLocalEnd(da, X, INSERT_VALUES, localX); CHKERRQ(info);

584:   info = DAVecGetArray(da, localX, (void**)&x); CHKERRQ(info);
585:   info = DAVecGetArray(da, localG, (void**)&g); CHKERRQ(info);

587:   info = DAGetCorners(da, &xs, &ys, TAO_NULL, &xm, &ym, TAO_NULL); CHKERRQ(info);
588:   info = DAGetGhostCorners(da, &gxs, &gys, TAO_NULL, &gxm, &gym, TAO_NULL); CHKERRQ(info);

590:   xe = gxs + gxm - 1;
591:   ye = gys + gym - 1;
592:   for (j = ys; j < ye; j++) {
593:     for (i = xs; i < xe; i++) {

595:       /* lower triangle contribution */
596:       dvdx = (x[j][i] - x[j][i+1]) / hx;
597:       dvdy = (x[j][i] - x[j+1][i]) / hy;
598:       sqGrad = dvdx * dvdx + dvdy * dvdy;
599:       avgWq = (2.0 * wq[i] + wq[i+1]) / 6.0;
600:       avgWV = (wl[i]*x[j][i] + wl[i+1]*x[j][i+1] + wl[i]*x[j+1][i]) / 3.0;
601:       fl = avgWq * sqGrad - avgWV;

603:       dvdx = dvdx * hy * avgWq;
604:       dvdy = dvdy * hx * avgWq;
605:       g[j][i] +=  ( dvdx + dvdy ) - wl[i] * aread3;
606:       g[j][i+1] += ( -dvdx ) - wl[i+1] * aread3;
607:       g[j+1][i] += ( -dvdy ) - wl[i] * aread3;

609:       /* upper triangle contribution */
610:       dvdx = (x[j+1][i+1] - x[j+1][i]) / hx;
611:       dvdy = (x[j+1][i+1] - x[j][i+1]) / hy;
612:       sqGrad = dvdx * dvdx + dvdy * dvdy;
613:       avgWq = (2.0 * wq[i+1] + wq[i]) / 6.0;
614:       avgWV = (wl[i+1]*x[j][i+1] + wl[i]*x[j+1][i] + wl[i+1]*x[j+1][i+1]) / 3.0;
615:       fu = avgWq * sqGrad - avgWV;

617:       dvdx = dvdx * hy * avgWq;
618:       dvdy = dvdy * hx * avgWq;
619:       g[j][i+1] += (-dvdy) - wl[i+1] * aread3;
620:       g[j+1][i] +=  (-dvdx) - wl[i] * aread3;
621:       g[j+1][i+1] += ( dvdx + dvdy ) - wl[i+1] * aread3;

623:       floc += area * (fl + fu);

625:     }
626:   }

628:   info = MPI_Allreduce(&floc, f, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); CHKERRQ(info);

630:   info = DAVecRestoreArray(da, localX, (void**)&x); CHKERRQ(info);
631:   info = DAVecRestoreArray(da, localG, (void**)&g); CHKERRQ(info);

633:   info = DALocalToGlobalBegin(da, localG, G); CHKERRQ(info);
634:   info = DALocalToGlobalEnd(da, localG, G); CHKERRQ(info);

636:   info = DARestoreLocalVector(da, &localX); CHKERRQ(info);
637:   info = DARestoreLocalVector(da, &localG); CHKERRQ(info);

639:   info = PetscLogFlops((xe-xs) * (ye-ys) * 67 + 1); CHKERRQ(info);
640:   return 0;
641: } /* WholeJBearFunctionGradient */


646: /*
647:   WholeJBearHessian - Evaluates Hessian over the whole grid

649:   Input:
650:     daapplication - TAO application object
651:     da  - distributed array
652:     X   - the current point, at which the function and gradient are evaluated
653:     ptr - user-defined application context

655:   Output:
656:     H - Hessian at X
657: */
658: static int WholeJBearHessian(TAO_APPLICATION daapplication, DA da, Vec X, Mat H, void *ptr) {

660:   AppCtx *user = (AppCtx*)ptr;
661:   int info;
662:   PetscInt i, j, ind[4];
663:   PetscInt xs, xm, gxs, gxm, xe, ys, ym, gys, gym, ye;
664:   double smallH[4][4];
665:   double wql, wqu;
666:   double dydx, dxdy;
667:   double *wq;
668:   double wqldydx,wqldxdy,wqudydx,wqudxdy;
669:   PetscTruth assembled;

671:   wq = user->wq;
672:   dydx = user->hy / user->hx;
673:   dxdy = user->hx / user->hy;

675:   info = MatAssembled(H,&assembled); CHKERRQ(info);
676:   if (assembled){info = MatZeroEntries(H);  CHKERRQ(info);}


679:   info = DAGetCorners(da, &xs, &ys, TAO_NULL, &xm, &ym, TAO_NULL); CHKERRQ(info);
680:   info = DAGetGhostCorners(da, &gxs, &gys, TAO_NULL, &gxm, &gym, TAO_NULL); CHKERRQ(info);

682:   xe = gxs + gxm - 1;
683:   ye = gys + gym - 1;
684:   for (j = ys; j < ye; j++) {
685:     for (i = xs; i < xe; i++) {

687:       wql = (2.0 * wq[i] + wq[i+1]) / 6.0;
688:       wqu = (wq[i] + 2.0 * wq[i+1]) / 6.0;

690:       wqldydx = wql * dydx;
691:       wqldxdy = wql * dxdy;
692:       wqudydx = wqu * dydx;
693:       wqudxdy = wqu * dxdy;

695:           /* Hessian contribution at 0,0 */
696:       smallH[0][0] = wqldxdy + wqldydx;
697:       smallH[0][1] = smallH[1][0] = -wqldydx;
698:       smallH[0][2] = smallH[2][0] = -wqldxdy;
699:       smallH[0][3] = smallH[3][0] = 0.0;

701:           /* Hessian contribution at 1,0 */
702:       smallH[1][1] = (wqldydx + wqudxdy);
703:       smallH[1][2] = smallH[2][1] = 0.0;
704:       smallH[1][3] = smallH[3][1] = -wqudxdy;

706:           /* Hessian contribution at 0,1 */
707:       smallH[2][2] = (wqldxdy + wqudydx);
708:       smallH[2][3] = smallH[3][2] = -wqudydx;

710:           /* Hessian contribution at 1,1 */
711:       smallH[3][3] = wqudydx + wqudxdy;

713:       ind[0] = (j-gys) * gxm + (i-gxs);
714:       ind[1] = ind[0] + 1;
715:       ind[2] = ind[0] + gxm;
716:       ind[3] = ind[2] + 1;
717:       info = MatSetValuesLocal(H,4,ind,4,ind,(PetscScalar*)smallH,ADD_VALUES); CHKERRQ(info);

719:     }
720:   }


723:   info = MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY); CHKERRQ(info);
724:   info = MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY); CHKERRQ(info);
725:   info = MatSetOption(H, MAT_SYMMETRIC, PETSC_TRUE); CHKERRQ(info);


728:   info = PetscLogFlops((xe-xs) * (ye-ys) * 14 + 2); CHKERRQ(info);
729:   return 0;

731: } /* WholeJBearHessian */