Actual source code: eptorsion3.c

  1: /*$Id: eptorsion.c, v 1.1 2002/08/08 10:30 lopezca@mauddib.mcs.anl.gov $*/

  3: /* Program usage: mpirun -np <proc> eptorsion [-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[] = "This example is based on the Elastic-Plastic Torsion (dept)\n\
 16: problem from the MINPACK-2 test suite.\n\
 17: The command line options are:\n\
 18:   -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
 19:   -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
 20:   -nlevels <nlevels>, where <nlevels> = number of levels in multigrid\n\
 21:   -byelement, if computation is made by functions on rectangular elements\n\
 22:   -adic, if AD is used (AD is not used by default)\n\
 23:   -u1 <u1>, where <u1> = upper limit in the 1st coordinate direction\n\
 24:   -u2 <u2>, where <u2> = upper limit in the 2nd coordinate direction\n\
 25:   -par <param>, where <param> = angle of twist per unit length\n\n";

 27: /*T
 28:    Concepts: TAO - Solving a bounded minimization problem
 29:    Routines: TaoInitialize(); TaoFinalize();
 30:    Routines: TaoCreate(); TaoDestroy();
 31:    Routines: DAApplicationCreate(); DAApplicationDestroy();
 32:    Routines: DAAppSetVariableBoundsRoutine();
 33:    Routines: DAAppSetElementObjectiveAndGradientRoutine();
 34:    Routines: DAAppSetElementHessianRoutine();
 35:    Routines: DAAppSetObjectiveAndGradientRoutine();
 36:    Routines: DAAppSetADElementFunctionGradient();
 37:    Routines: DAAppSetHessianRoutine();
 38:    Routines: TaoSetOptions();
 39:    Routines: TaoGetSolutionStatus(); TaoDAAppSolve();
 40:    Routines: DAAppSetMonitor(); TaoView();
 41:    Routines: DAAppGetSolution();
 42:    Routines: DAAppGetInterpolationMatrix();
 43:    Processors: n
 44: T*/

 46: /*
 47:    User-defined application context - contains data needed by the
 48:    application-provided call-back routines.
 49: */
 50: typedef struct {
 51:   InactiveDouble      param;
 52:   InactiveDouble      hx, hy;        /* increment size in both directions */
 53:   InactiveDouble      area;          /* area of the triangles */
 54: } ADFGCtx;

 56: typedef struct {
 57:   PetscReal      param;          /* 'c' parameter */
 58:   PetscReal      u1, u2;         /* domain upper limits (lower limits = 0) */
 59:   double      hx, hy;        /* increment size in both directions */
 60:   double      area;          /* area of the triangles */
 61:   ADFGCtx     fgctx;         /* Used only when an ADIC generated gradient is used */
 62: } AppCtx;
 63: int ad_EPTorsLocalFunction(PetscInt[2], DERIV_TYPE[4], DERIV_TYPE*, void*);

 65: /* User-defined routines found in this file */
 66: static int AppCtxInitialize(void *ptr);
 67: static int FormInitialGuess(DA, Vec);

 69: static int EPTorsLocalFunctionGradient(PetscInt[2], double x[4], double *f, double g[4], void *ptr);
 70: static int EPTorsLocalHessian(PetscInt[2], double x[4], double H[4][4], void *ptr);

 72: static int WholeEPTorsFunctionGradient(TAO_APPLICATION,DA,Vec,double *,Vec,void*);
 73: static int WholeEPTorsHessian(TAO_APPLICATION,DA,Vec,Mat,void*);

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

 77: static int MyGridMonitorBefore(TAO_APPLICATION, DA, PetscInt, void *);


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

 84:   int             info;                           /* used to check for functions returning nonzeros */
 85:   PetscInt             mx,my,Nx,Ny;
 86:   double          ff,gnorm;
 87:   PetscInt             iter, nlevels;                                                /* multigrid levels */
 88:   DA              DAarray[20];
 89:   Vec             X;
 90:   PetscTruth      flg, PreLoad = PETSC_TRUE;                                               /* flags */
 91:   TaoMethod       method = "tao_gpcg";                                       /* minimization method */
 92:   AppCtx          user;                                                /* user-defined work context */
 93:   TAO_SOLVER      tao;                                                 /* TAO_SOLVER solver context */
 94:   TAO_APPLICATION EPTorsApp;                                            /* The PETSc application */
 95:   TaoTerminateReason reason;
 96:   KSP ksp;
 97:   PC pc;

 99:   /* Initialize TAO */
100:   PetscInitialize(&argc, &argv, (char *)0, help);
101:   TaoInitialize(&argc, &argv, (char *)0, help);

103:   PreLoadBegin(PreLoad,"Solve");
104:   
105:   info = AppCtxInitialize((void*)&user); CHKERRQ(info);

107:   nlevels=5;
108:   info = PetscOptionsGetInt(PETSC_NULL,"-nlevels",&nlevels,&flg); CHKERRQ(info);
109:   mx = my = 11;                               /* these correspond to 10 segments on each dimension */
110:   info = PetscOptionsGetInt(TAO_NULL, "-mx", &mx, &flg); CHKERRQ(info);
111:   info = PetscOptionsGetInt(TAO_NULL, "-my", &my, &flg); CHKERRQ(info);
112:   if (PreLoadIt == 0) {
113:     nlevels = 1; mx = 11; my = 11; }

115:   PetscPrintf(MPI_COMM_WORLD,"\n---- Elastic-Plastic Torsion Problem -----\n\n");

117:   /* Let PETSc determine the vector distribution */
118:   Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;

120:   /* Create distributed array (DA) to manage parallel grid and vectors  */
121:   info = DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_BOX,mx,
122:                     my,Nx,Ny,1,1,PETSC_NULL,PETSC_NULL,&DAarray[0]); CHKERRQ(info);
123:   for (iter=1;iter<nlevels;iter++){
124:     info = DARefine(DAarray[iter-1],PETSC_COMM_WORLD,&DAarray[iter]); CHKERRQ(info);
125:   }

127:   /* Create TAO solver and set desired solution method */
128:   info = TaoCreate(MPI_COMM_WORLD,method,&tao); CHKERRQ(info);
129:   info = TaoApplicationCreate(PETSC_COMM_WORLD,&EPTorsApp); CHKERRQ(info);
130:   info = TaoAppSetDAApp(EPTorsApp, DAarray, nlevels ); CHKERRQ(info);
131:   /* Sets routines for function, gradient and bounds evaluation */
132:   info = DAAppSetVariableBoundsRoutine(EPTorsApp,DASetBounds,(void *)&user); CHKERRQ(info);

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

137:     /* Sets routines for function and gradient evaluation, element by element */
138:     info = PetscOptionsHasName(TAO_NULL, "-adic", &flg); CHKERRQ(info);
139:     if (flg) {
140:       info = DAAppSetADElementFunctionGradient(EPTorsApp,ad_EPTorsLocalFunction,192,(void *)&user.fgctx); CHKERRQ(info);
141:     } else {
142:       info = DAAppSetElementObjectiveAndGradientRoutine(EPTorsApp,EPTorsLocalFunctionGradient,42,(void *)&user); CHKERRQ(info);
143:     }
144:     /* Sets routines for Hessian evaluation, element by element */
145:     info = DAAppSetElementHessianRoutine(EPTorsApp,EPTorsLocalHessian,6,(void*)&user); CHKERRQ(info);

147:   } else {

149:     /* Sets routines for function and gradient evaluation, all in one routine */
150:     info = DAAppSetObjectiveAndGradientRoutine(EPTorsApp,WholeEPTorsFunctionGradient,(void *)&user); CHKERRQ(info);

152:     /* Sets routines for Hessian evaluation, all in one routine */
153:     info = DAAppSetHessianRoutine(EPTorsApp,WholeEPTorsHessian,(void*)&user); CHKERRQ(info);    
154:     
155:   }

157:   info = DAAppSetBeforeMonitor(EPTorsApp,MyGridMonitorBefore,(void*)&user); CHKERRQ(info);
158:   info = PetscOptionsHasName(TAO_NULL,"-tao_monitor", &flg); CHKERRQ(info);
159:   if (flg){
160:     info = DAAppPrintInterpolationError(EPTorsApp); CHKERRQ(info);
161:     info = DAAppPrintStageTimes(EPTorsApp); CHKERRQ(info);
162:   }
163:   info = TaoAppSetRelativeTolerance(EPTorsApp,1.0e-6); CHKERRQ(info);
164:   info = TaoSetTolerances(tao,0,0,0,0); CHKERRQ(info);
165:   info = TaoSetGradientTolerances(tao,0,0,0); CHKERRQ(info);

167:   info = TaoAppGetKSP(EPTorsApp,&ksp);CHKERRQ(info);
168:   info = KSPSetType(ksp,KSPCG); CHKERRQ(info);
169:   info = KSPSetTolerances(ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,100);CHKERRQ(info);
170:   info = KSPGetPC(ksp,&pc);CHKERRQ(info);
171:   info = PCSetType(pc,PCBJACOBI);CHKERRQ(info);

173:   /* Check for any tao command line options */
174:   info = TaoSetOptions(EPTorsApp, tao); CHKERRQ(info);

176:   info = DAAppGetSolution(EPTorsApp,0,&X); CHKERRQ(info);
177:   info = FormInitialGuess(DAarray[0],X); CHKERRQ(info);
178:   info = DAAppSetInitialSolution(EPTorsApp,X); CHKERRQ(info);
179:   
180:   /* SOLVE THE APPLICATION */
181:   info = TaoDAAppSolve(EPTorsApp, tao);  CHKERRQ(info);

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

190:   info = PetscOptionsHasName(PETSC_NULL,"-view_sol",&flg); CHKERRQ(info);
191:   if (flg){
192:     info = DAAppGetSolution(EPTorsApp,nlevels-1,&X); CHKERRQ(info);
193:     info=VecView(X,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(info);
194:   }

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

199:   /* Free TAO data structures */
200:   info = TaoDestroy(tao); CHKERRQ(info);
201:   info = TaoAppDestroy(EPTorsApp); CHKERRQ(info);

203:   /* Free PETSc data structures */
204:   for (iter=0;iter<nlevels;iter++){
205:     info = DADestroy(DAarray[iter]); CHKERRQ(info);
206:   }

208:   PreLoadEnd();

210:   /* Finalize TAO */
211:   TaoFinalize();
212:   PetscFinalize();

214:   return 0;
215: } /* main */



219: /*----- The following two routines
220:   MyGridMonitorBefore    MyGridMonitorAfter
221:   help diplay info of iterations at every grid level
222: */
225: static int MyGridMonitorBefore(TAO_APPLICATION myapp, DA da, PetscInt level, void *ctx) {

227:   AppCtx *user = (AppCtx*)ctx;
228:   PetscInt mx,my;
229:   int info;

231:   info = DAGetInfo(da,PETSC_NULL,&mx,&my,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
232:                    PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(info);
233:   user->hx = user->u1 / (mx - 1);
234:   user->hy = user->u2 / (my - 1);
235:   user->area = 0.5 * user->hx * user->hy;
236:   user->fgctx.hx   = user->hx;
237:   user->fgctx.hy   = user->hy;
238:   user->fgctx.area = user->area;
239:   user->fgctx.param = user->param;

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

243:   return 0;
244: }


247: /*------- USER-DEFINED: initialize the application context information -------*/
250: /*
251:   AppCtxInitialize - Sets initial values for the application context parameters

253:   Input:
254:     ptr - void user-defined application context

256:   Output:
257:     ptr - user-defined application context with the default or user-provided
258:              parameters
259: */
260: static int AppCtxInitialize(void *ptr) {

262:   AppCtx *user = (AppCtx*)ptr;
263:   PetscTruth    flg;            /* flag for PETSc calls */
264:   int info;

266:   /* Specify default parameters */
267:   user->param = 25.0;
268:   user->u1 = user->u2 = 1.0;

270:   /* Check for command line arguments that override defaults */
271:   info = PetscOptionsGetReal(TAO_NULL, "-par", &user->param, &flg); CHKERRQ(info);
272:   info = PetscOptionsGetReal(TAO_NULL, "-u1", &user->u1, &flg); CHKERRQ(info);
273:   info = PetscOptionsGetReal(TAO_NULL, "-u2", &user->u2, &flg); CHKERRQ(info);

275:   return 0;
276: } /* AppCtxInitialize */

280: static int FormInitialGuess(DA da, Vec X)
281: {
282:   int info;
283:   PetscInt    i, j, mx, my;
284:   PetscInt    xs, ys, xm, ym, xe, ye;
285:   PetscReal hx, hy, temp, val;
286:   double **x;

288:   /* Get local mesh boundaries */
289:   info = DAGetInfo(da,PETSC_NULL,&mx,&my,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
290:                    PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(info);
291:   hx = 1.0/(mx-1);  hy = 1.0/(my-1);

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

296:   info = DAVecGetArray(da, X, (void**)&x); CHKERRQ(info);
297:   /* Compute initial guess over locally owned part of mesh */
298:   for (j=ys; j<ye; j++) {  /*  for (j=0; j<my; j++) */
299:     temp = PetscMin(j+1,my-j)*hy;
300:     for (i=xs; i<xe; i++) {  /*  for (i=0; i<mx; i++) */
301:       val = PetscMin((PetscMin(i+1,mx-i))*hx,temp);
302:       x[j][i] = val;
303:     }
304:   }
305:   info = DAVecRestoreArray(da, X, (void**)&x); CHKERRQ(info);

307:   return 0;
308: }


311: /*------- USER-DEFINED: set the upper and lower bounds for the variables  -------*/

315: /*
316:   FormBounds - Forms bounds on the variables

318:   Input:
319:     user - user-defined application context

321:   Output:
322:     XL - vector of lower bounds
323:     XU - vector of upper bounds
324: */
325: static int DASetBounds(TAO_APPLICATION daapplication, DA da, Vec XL, Vec XU, void *ptr)
326: {
327:   AppCtx *user = (AppCtx*)ptr;
328:   int info;
329:   PetscInt i, j, xs, xm, ys, ym;
330:   double hx, hy, u1, u2, dist, d1, d2, hd, vd;
331:   double **xl, **xu;

333:   hx = user->hx;
334:   hy = user->hy;
335:   u1 = user->u1;
336:   u2 = user->u2;

338:   info = DAVecGetArray(da, XL, (void**)&xl); CHKERRQ(info);
339:   info = DAVecGetArray(da, XU, (void**)&xu); CHKERRQ(info);
340:   info = DAGetCorners(da, &xs, &ys, TAO_NULL, &xm, &ym, TAO_NULL); CHKERRQ(info);

342:   for (j = ys; j < ys+ym; j++){
343:     for (i = xs; i < xs+xm; i++){
344:       d1 = i * hx; d2 = u1 - d1; hd = PetscMin(d1,d2);
345:       d1 = j * hy; d2 = u2 - d1; vd = PetscMin(d1,d2);
346:       dist = PetscMin(hd,vd);
347:       xl[j][i] = -dist;
348:       xu[j][i] = dist;
349:     }
350:   }

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

355:   info = PetscLogFlops(xm * ym * 4); CHKERRQ(info);
356:   return 0;

358: } /* DASetBounds */


363: /*
364:   EPTorsLocalFunctionGradient - Evaluates function and gradient over the 
365:       local rectangular element

367:   Input:
368:     coor - vector with the indices of the position of current element
369:              in the first, second and third directions
370:     x - current point (values over the current rectangular element)
371:     df - degrees of freedom at each point
372:     ptr - user-defined application context

374:   Output:
375:     f - value of the objective funtion at the local rectangular element
376:     g - gradient of the local function
377: */
378: static int EPTorsLocalFunctionGradient(PetscInt coor[2], double x[4], double *f, double g[4], void *ptr) {

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

382:   double fquad, flin;
383:   double hx, hy, dvdx, dvdy, area;
384:   double cdiv3, cnt;

386:   cdiv3 = user->param / 3.0;
387:   hx = user->hx;
388:   hy = user->hy;
389:   area = user->area;
390:   cnt = area * cdiv3;

392:   /* lower triangle contribution */
393:   dvdx = (x[0] - x[1]) / hx;
394:   dvdy = (x[0] - x[2]) / hy;
395:   fquad = dvdx * dvdx + dvdy * dvdy;
396:   flin = x[0] + x[1] + x[2];

398:   dvdx = 0.5 * dvdx * hy;
399:   dvdy = 0.5 * dvdy * hx;
400:   g[0] = dvdx + dvdy - cnt;
401:   g[1] = -dvdx - 2.0 * cnt;
402:   g[2] = -dvdy - 2.0 * cnt;

404:   /* upper triangle contribution */
405:   dvdx = (x[3] - x[2]) / hx;
406:   dvdy = (x[3] - x[1]) / hy;
407:   fquad += dvdx * dvdx + dvdy * dvdy;
408:   flin += x[1] + x[2] + x[3];

410:   dvdx = 0.5 * dvdx * hy;
411:   dvdy = 0.5 * dvdy * hx;
412:   g[1] += -dvdy;
413:   g[2] += -dvdx;
414:   g[3] = dvdx + dvdy - cnt;

416:   *f = area * (0.5 * fquad - flin * cdiv3);

418:   return 0;
419: } /* EPTorsLocalFunctionGradient */



423: /*------- USER-DEFINED: routine to evaluate the Hessian
424:            at a local (rectangular element) level       -------*/
427: /*
428:   EPTorsLocalHessian - Computes the Hessian of the local (partial) function
429:          defined over the current rectangle

431:   Input:
432:     coor - vector with the indices of the position of current element
433:              in the first, second and third directions
434:     x - current local solution (over the rectangle only)
435:     df - degrees of freedom at each point
436:     ptr - user-defined application context

438:   Output:
439:     H - Hessian matrix of the local function (wrt the four
440:            points of the rectangle only)
441: */
442: static int EPTorsLocalHessian(PetscInt coor[2], double x[4], double H[4][4], void *ptr) {

444:   AppCtx *user = (AppCtx*)ptr;
445:   double hx, hy, dxdy, dydx;
446:   double diagxy, bandxy, bandyx;

448:   hx = user->hx;
449:   hy = user->hy;
450:   dxdy = hx/hy;
451:   dydx = hy/hx;
452:   diagxy = 0.5 * (dxdy + dydx);
453:   bandxy = -0.5 * dxdy;
454:   bandyx = -0.5 * dydx;

456:           /* Hessian contribution at 0,0 */
457:   H[0][0] = diagxy;
458:   H[0][1] =  H[1][0] = bandyx;
459:   H[0][2] =  H[2][0] = bandxy;
460:   H[0][3] =  H[3][0] = 0.0;

462:           /* Hessian contribution at 1,0 */
463:   H[1][1] = diagxy;
464:   H[1][2] =  H[2][1] = 0.0;
465:   H[1][3] =  H[3][1] = bandxy;

467:           /* Hessian contribution at 0,1 */
468:   H[2][2] = diagxy;
469:   H[2][3] =  H[3][2] = bandyx;

471:           /* Hessian contribution at 1,1 */
472:   H[3][3] = diagxy;

474:   return 0;

476: } /* EPTorsLocalHessian */


479: /*------- USER-DEFINED: routine to evaluate the function 
480:           and gradient at the whole grid             -------*/
483: /*
484:   WholeEPTorsFunctionGradient - Evaluates function and gradient over the 
485:       whole grid

487:   Input:
488:     daapplication - TAO application object
489:     da  - distributed array
490:     X   - the current point, at which the function and gradient are evaluated
491:     ptr - user-defined application context

493:   Output:
494:     f - value of the objective funtion at X
495:     G - gradient at X
496: */
497: static int WholeEPTorsFunctionGradient(TAO_APPLICATION daapplication, DA da, Vec X, double *f, Vec G, void *ptr) {

499:   AppCtx *user = (AppCtx*)ptr;
500:   Vec localX, localG;
501:   int info;
502:   PetscInt i, j;
503:   PetscInt xs, xm, gxs, gxm, xe, ys, ym, gys, gym, ye;
504:   double **x, **g;
505:   double floc = 0.0;
506:   PetscScalar zero = 0.0;

508:   double fquad, flin;
509:   double hx, hy, dvdx, dvdy, area;
510:   double cdiv3, cnt;

512:   cdiv3 = user->param / 3.0;
513:   hx = user->hx;
514:   hy = user->hy;
515:   area = user->area;
516:   cnt = area * cdiv3;

518:   info = DAGetLocalVector(da, &localX); CHKERRQ(info);
519:   info = DAGetLocalVector(da, &localG); CHKERRQ(info);
520:   info = VecSet(G, zero); CHKERRQ(info);
521:   info = VecSet(localG, zero); CHKERRQ(info);

523:   info = DAGlobalToLocalBegin(da, X, INSERT_VALUES, localX); CHKERRQ(info);
524:   info = DAGlobalToLocalEnd(da, X, INSERT_VALUES, localX); CHKERRQ(info);

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

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

532:   xe = gxs + gxm - 1;
533:   ye = gys + gym - 1;
534:   for (j = ys; j < ye; j++) {
535:     for (i = xs; i < xe; i++) {

537:       /* lower triangle contribution */
538:       dvdx = (x[j][i] - x[j][i+1]) / hx;  
539:       dvdy = (x[j][i] - x[j+1][i]) / hy;
540:       fquad = dvdx * dvdx + dvdy * dvdy;
541:       flin = x[j][i] + x[j][i+1] + x[j+1][i];

543:       dvdx = 0.5 * dvdx * hy;
544:       dvdy = 0.5 * dvdy * hx;
545:       g[j][i] += dvdx + dvdy - cnt;
546:       g[j][i+1] += -dvdx - 2.0 * cnt;
547:       g[j+1][i] += -dvdy - 2.0 * cnt;

549:       /* upper triangle contribution */
550:       dvdx = (x[j+1][i+1] - x[j+1][i]) / hx;
551:       dvdy = (x[j+1][i+1] - x[j][i+1]) / hy;
552:       fquad += dvdx * dvdx + dvdy * dvdy;
553:       flin += x[j][i+1] + x[j+1][i] + x[j+1][i+1];

555:       dvdx = 0.5 * dvdx * hy;
556:       dvdy = 0.5 * dvdy * hx;
557:       g[j][i+1] += -dvdy;
558:       g[j+1][i] += -dvdx;
559:       g[j+1][i+1] += dvdx + dvdy - cnt;

561:       floc += area * (0.5 * fquad - flin * cdiv3);

563:     }
564:   }

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

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

571:   info = DALocalToGlobalBegin(da, localG, G); CHKERRQ(info);
572:   info = DALocalToGlobalEnd(da, localG, G); CHKERRQ(info);

574:   info = DARestoreLocalVector(da, &localX); CHKERRQ(info);
575:   info = DARestoreLocalVector(da, &localG); CHKERRQ(info);

577:   info = PetscLogFlops((xe-xs) * (ye-ys) * 47 + 2); CHKERRQ(info);
578:   return 0;
579: } /* WholeEPTorsFunctionGradient */


582: /*------- USER-DEFINED: routine to evaluate the Hessian 
583:           at the whole grid             -------*/
586: /*
587:   WholeEPTorsHessian - Evaluates Hessian over the whole grid

589:   Input:
590:     daapplication - TAO application object
591:     da  - distributed array
592:     X   - the current point, at which the function and gradient are evaluated
593:     ptr - user-defined application context

595:   Output:
596:     H - Hessian at X
597: */
598: static int WholeEPTorsHessian(TAO_APPLICATION daapplication, DA da, Vec X, Mat H, void *ptr) {

600:   AppCtx *user = (AppCtx*)ptr;
601:   int info;
602:   PetscInt i, j, ind[4];
603:   PetscInt xs, xm, gxs, gxm, xe, ys, ym, gys, gym, ye;
604:   double smallH[4][4];

606:   double hx, hy, dxdy, dydx;
607:   double diagxy, bandxy, bandyx;
608:   PetscTruth assembled;

610:   hx = user->hx;
611:   hy = user->hy;
612:   dxdy = hx/hy;
613:   dydx = hy/hx;
614:   diagxy = 0.5 * (dxdy + dydx);
615:   bandxy = -0.5 * dxdy;
616:   bandyx = -0.5 * dydx;

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


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

625:   xe = gxs + gxm - 1;
626:   ye = gys + gym - 1;
627:   for (j = ys; j < ye; j++) {
628:     for (i = xs; i < xe; i++) {

630:           /* Hessian contribution at 0,0 */
631:       smallH[0][0] = diagxy;
632:       smallH[0][1] = smallH[1][0] = bandyx;
633:       smallH[0][2] = smallH[2][0] = bandxy;
634:       smallH[0][3] = smallH[3][0] = 0.0;

636:           /* Hessian contribution at 1,0 */
637:       smallH[1][1] = diagxy;
638:       smallH[1][2] = smallH[2][1] = 0.0;
639:       smallH[1][3] = smallH[3][1] = bandxy;

641:           /* Hessian contribution at 0,1 */
642:       smallH[2][2] = diagxy;
643:       smallH[2][3] = smallH[3][2] = bandyx;

645:           /* Hessian contribution at 1,1 */
646:       smallH[3][3] = diagxy;

648:       ind[0] = (j-gys) * gxm + (i-gxs);
649:       ind[1] = ind[0] + 1;
650:       ind[2] = ind[0] + gxm;
651:       ind[3] = ind[2] + 1;
652:       info = MatSetValuesLocal(H,4,ind,4,ind,(PetscScalar*)smallH,ADD_VALUES); CHKERRQ(info);

654:     }
655:   }

657:   info = MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY); CHKERRQ(info);
658:   info = MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY); CHKERRQ(info);
659:   info = MatSetOption(H, MAT_SYMMETRIC, PETSC_TRUE); CHKERRQ(info);


662:   info = PetscLogFlops(6); CHKERRQ(info);
663:   return 0;

665: } /* WholeEPTorsHessian */