Actual source code: combustion3.c

  1: /*$Id: combustion.c,v 1.1 2002/08/08 9:39 lopezca@mauddib.mcs.anl.gov $*/

  3: /* Program usage: mpirun -np <proc> combustion [-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[] = "Steady-State Combustion.\n\
 16: We solve the Steady-State Combustion problem (MINPACK-2 test suite) in a 2D\n\
 17: rectangular domain, using distributed arrays (DAs) to partition the parallel grid.\n\
 18: The command line options include:\n\
 19:   -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
 20:   -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
 21:   -nlevels <nlevels>, where <nlevels> = number of levels in multigrid\n\
 22:   -byelement, if computation is made by functions on rectangular elements\n\
 23:   -adic, if AD is used (AD is not used by default)\n\
 24:   -par <parameter>, where <parameter> indicates the problem's nonlinearity\n\
 25:      parameter lambda (0 <= par <= 6.81)\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: TaoAppSetOptions();
 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: */

 51: typedef struct {
 52:   InactiveDouble      param;
 53:   InactiveDouble      hx, hy;        /* increment size in both directions */
 54:   InactiveDouble      area;          /* area of the triangles */
 55: } ADFGCtx;

 57: typedef struct {
 58:   PetscReal  param;          /* nonlinearity parameter */
 59:   double  hx, hy, area;   /* increments and area of the triangle */
 60:   PetscInt     mx, my;         /* discretization including boundaries */
 61:   ADFGCtx fgctx;          /* Used only when an ADIC generated gradient is used */
 62: } AppCtx;
 63: int ad_CombLocalFunction(PetscInt[2], DERIV_TYPE[4], DERIV_TYPE*, void*);

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

 69: static int CombLocalFunctionGradient(PetscInt[3], double x[4], double *f, double g[4], void *ptr);
 70: static int WholeCombFunctionGradient(TAO_APPLICATION,DA,Vec,double *,Vec,void*);

 72: static int CombLocalHessian(PetscInt[3], double x[4], double H[4][4], void *ptr);
 73: static int WholeCombHessian(TAO_APPLICATION,DA,Vec,Mat,void*);

 75: static int DAFixBoundary(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        Nx,Ny,iter;
 86:   PetscInt        nlevels;                                           /* multigrid levels */
 87:   double          ff,gnorm;
 88:   DA              DAarray[20];
 89:   Vec             X;
 90:   KSP             ksp;
 91:   PetscTruth      flg, PreLoad = PETSC_TRUE;                                    /* flags */
 92:   AppCtx          user;                                     /* user-defined work context */
 93:   TaoMethod       method = "tao_tron";                            /* minimization method */
 94:   TAO_SOLVER      tao;                                      /* TAO_SOLVER solver context */
 95:   TAO_APPLICATION CombApp;                                   /* The PETSc application */
 96:   TaoTerminateReason reason;

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

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

106:   nlevels=5;
107:   info = PetscOptionsGetInt(PETSC_NULL,"-nlevels",&nlevels,&flg); CHKERRQ(info);
108:   if (PreLoadIt == 0) {
109:     nlevels = 1; user.mx = 11; user.my = 11;}

111:   PetscPrintf(MPI_COMM_WORLD,"\n---- Steady-State Combustion Problem -----\n\n");

113:   /* Let PETSc determine the vector distribution */
114:   Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;

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

123:   /* Create TAO solver and set desired solution method */
124:   info = TaoCreate(MPI_COMM_WORLD,method,&tao); CHKERRQ(info);

126:   info = TaoApplicationCreate(PETSC_COMM_WORLD, &CombApp); CHKERRQ(info);
127:   info = TaoAppSetDAApp(CombApp,DAarray,nlevels); CHKERRQ(info);

129:   /* Sets routines for function, gradient and bounds evaluation */
130:   info = DAAppSetVariableBoundsRoutine(CombApp,DAFixBoundary,(void *)&user); CHKERRQ(info);

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

135:     /* Sets routines for function and gradient evaluation, element by element */
136:     info = PetscOptionsHasName(TAO_NULL, "-adic", &flg); CHKERRQ(info);
137:     if (flg) {
138:       info = DAAppSetADElementFunctionGradient(CombApp,ad_CombLocalFunction,228,(void *)&user.fgctx); CHKERRQ(info);
139:     } else {
140:       info = DAAppSetElementObjectiveAndGradientRoutine(CombApp,CombLocalFunctionGradient,51,(void *)&user); CHKERRQ(info);
141:     }
142:     /* Sets routines for Hessian evaluation, element by element */
143:     info = DAAppSetElementHessianRoutine(CombApp,CombLocalHessian,21,(void*)&user); CHKERRQ(info);

145:   } else {

147:     /* Sets routines for function and gradient evaluation, all in one routine */
148:     info = DAAppSetObjectiveAndGradientRoutine(CombApp,WholeCombFunctionGradient,(void *)&user); CHKERRQ(info);

150:     /* Sets routines for Hessian evaluation, all in one routine */
151:     info = DAAppSetHessianRoutine(CombApp,WholeCombHessian,(void*)&user); CHKERRQ(info);    
152:     
153:   }

155:   info = DAAppSetBeforeMonitor(CombApp,MyGridMonitorBefore,(void*)&user); CHKERRQ(info);
156:   info = PetscOptionsHasName(TAO_NULL,"-tao_monitor", &flg); CHKERRQ(info);
157:   if (flg){
158:     info = DAAppPrintInterpolationError(CombApp); CHKERRQ(info);
159:     info = DAAppPrintStageTimes(CombApp); CHKERRQ(info);
160:   }


163:   info = TaoAppSetRelativeTolerance(CombApp,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(CombApp,&ksp);CHKERRQ(info);
168:   info = KSPSetType(ksp,KSPCG); CHKERRQ(info);
169:   info = KSPSetTolerances(ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,100);CHKERRQ(info);

171:   /* Check for any tao command line options */
172:   info = TaoSetOptions(CombApp, tao); CHKERRQ(info);

174:   info = DAAppGetSolution(CombApp,0,&X); CHKERRQ(info);
175:   info = FormInitialGuess(DAarray[0],X,&user); CHKERRQ(info);
176:   info = DAAppSetInitialSolution(CombApp,X); CHKERRQ(info);

178:   /* SOLVE THE APPLICATION */
179:   info = TaoDAAppSolve(CombApp, tao);  CHKERRQ(info);

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

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

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

197:   /* Free TAO data structures */
198:   info = TaoDestroy(tao); CHKERRQ(info);
199:   info = TaoAppDestroy(CombApp); CHKERRQ(info);

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

206:   PreLoadEnd();

208:   /* Finalize TAO */
209:   TaoFinalize();
210:   PetscFinalize();

212:   return 0;
213: } /* main */



217: /*----- The following two routines
218:    
219:   MyGridMonitorBefore    MyGridMonitorAfter

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

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

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

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->mx = mx;
234:   user->my = my;
235:   user->hx = 1.0 / (user->mx - 1);
236:   user->hy = 1.0 / (user->my - 1);
237:   user->area = 0.5 * user->hx * user->hy;
238:   user->fgctx.hx   = user->hx;
239:   user->fgctx.hy   = user->hy;
240:   user->fgctx.area = user->area;
241:   user->fgctx.param = user->param;

243:   PetscPrintf(MPI_COMM_WORLD,"Grid: %d,    mx: %d     my: %d   \n",level,mx,my);
244:   return 0;
245: }

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

251: /*
252:   AppCtxInitialize - Sets initial values for the application context parameters

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

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

263:   AppCtx *user = (AppCtx*)ptr;
264:   PetscReal     LambdaMax = 6.81, LambdaMin = 0.0;  /* bounds on parameter lambda */
265:   PetscTruth    flg;            /* flag for PETSc calls */
266:   int info;

268:   /* Specify dimension of the problem */
269:   user->param = 5.0;
270:   user->mx = 11;
271:   user->my = 11;

273:   /* Check for any command line arguments that override defaults */
274:   info = PetscOptionsGetReal(TAO_NULL, "-par", &user->param, &flg); CHKERRQ(info);
275:   if (user->param >= LambdaMax || user->param <= LambdaMin) {
276:     SETERRQ(1,"Lambda is out of range.");
277:   }
278:   info = PetscOptionsGetInt(PETSC_NULL,"-mx",&user->mx,&flg); CHKERRQ(info);
279:   info = PetscOptionsGetInt(PETSC_NULL,"-my",&user->my,&flg); CHKERRQ(info);

281:   user->hx = 1.0 / (user->mx - 1);
282:   user->hy = 1.0 / (user->my - 1);
283:   user->area = 0.5 * user->hx * user->hy;
284:   info = PetscLogFlops(6); CHKERRQ(info);

286:   return 0;
287: } /* AppCtxInitialize */



293: static int FormInitialGuess(DA da, Vec X, AppCtx *ctx)
294: {
295:   int    info;
296:   PetscInt i, j, mx, my;
297:   PetscInt xs, ys, xm, ym, xe, ye;
298:   PetscReal hx, hy, temp, val, lambda;
299:   double **x;

301:   lambda = ctx->param;
302:   lambda = lambda/(lambda+1.0);

304:   /* Get local mesh boundaries */
305:   info = DAGetInfo(da,PETSC_NULL,&mx,&my,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
306:                    PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(info);
307:   hx = 1.0/(mx-1);  hy = 1.0/(my-1);

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

312:   info = DAVecGetArray(da, X, (void**)&x); CHKERRQ(info);
313:   /* Compute initial guess over locally owned part of mesh */
314:   for (j=ys; j<ye; j++) {  /*  for (j=0; j<my; j++) */
315:     temp = PetscMin(j+1,my-j)*hy;
316:     for (i=xs; i<xe; i++) {  /*  for (i=0; i<mx; i++) */
317:       val = lambda*sqrt(PetscMin((PetscMin(i+1,mx-i))*hx,temp));
318:       x[j][i] = val;
319:     }
320:   }
321:   info = DAVecRestoreArray(da, X, (void**)&x); CHKERRQ(info);

323:   return 0;
324: }


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

331: /*
332:   FormBounds - Forms bounds on the variables

334:   Input:
335:     user - user-defined application context

337:   Output:
338:     XL - vector of lower bounds
339:     XU - vector of upper bounds

341: */
342: static int DAFixBoundary(TAO_APPLICATION daapplication, DA da, Vec XL, Vec XU, void *ptr)
343: {
344:   AppCtx *user = (AppCtx*)ptr;
345:   int info;
346:   PetscInt i, j, mx, my, xs, xm, ys, ym;
347:   double lb = -TAO_INFINITY;
348:   double ub = TAO_INFINITY;
349:   double **xl, **xu;

351:   mx = user->mx;  /* number of points including the boundary */
352:   my = user->my;

354:   info = DAVecGetArray(da, XL, (void**)&xl); CHKERRQ(info);
355:   info = DAVecGetArray(da, XU, (void**)&xu); CHKERRQ(info);
356:   info = DAGetCorners(da, &xs, &ys, TAO_NULL, &xm, &ym, TAO_NULL); CHKERRQ(info);

358:   /* Compute initial guess over locally owned part of the grid */
359:   for (j = ys; j < ys+ym; j++){
360:     for (i = xs; i < xs+xm; i++){
361:       if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) {
362:         xl[j][i] = xu[j][i] = 0.0;
363:       } else {
364:         xl[j][i] = lb;
365:         xu[j][i] = ub;
366:       }
367:     }
368:   }

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

373:   return 0;
374: } /* DAFixBoundary */


377: /*------- USER-DEFINED: routine to evaluate the function and gradient
378:            at a local (rectangular element) level              -------*/

382: /*
383:   CombLocalFunctionGradient - Evaluates function and gradient over the 
384:       local rectangular element

386:   Input:
387:     coor - vector with the indices of the position of current element
388:              in the first, second and third directions
389:     x - current point (values over the current rectangular element)
390:     ptr - user-defined application context

392:   Output:
393:     f - value of the objective funtion at the local rectangular element
394:     g - gradient of the local function

396: */
397: static int CombLocalFunctionGradient(PetscInt coor[2], double x[4], double *f, double g[4], void *ptr) {

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

401:   double lambdad3, hx, hy, area;
402:   double fquad, fexp, dvdx, dvdy;

404:   lambdad3 = user->param / 3.0;
405:   hx = user->hx;
406:   hy = user->hy;
407:   area = user->area;

409:   /* lower triangle contribution */
410:   dvdx = (x[0] - x[1]) / hx;
411:   dvdy = (x[0] - x[2]) / hy;
412:   fquad = dvdx * dvdx + dvdy * dvdy;
413:   fexp = exp(x[0]) + exp(x[1]) + exp(x[2]);

415:   dvdx = 0.5 * dvdx * hy;
416:   dvdy = 0.5 * dvdy * hx;
417:   g[0] = dvdx + dvdy - exp(x[0]) * area * lambdad3;
418:   g[1] = -dvdx - 2.0 * exp(x[1]) * area * lambdad3;
419:   g[2] = -dvdy - 2.0 * exp(x[2]) * area * lambdad3;

421:   /* upper triangle contribution */
422:   dvdx = (x[3] - x[2]) / hx;
423:   dvdy = (x[3] - x[1]) / hy;
424:   fquad += dvdx * dvdx + dvdy * dvdy;
425:   fexp += exp(x[1]) + exp(x[2]) + exp(x[3]);

427:   dvdx = 0.5 * dvdx * hy;
428:   dvdy = 0.5 * dvdy * hx;
429:   g[1] += -dvdy;
430:   g[2] += -dvdx;
431:   g[3] = dvdx + dvdy - exp(x[3]) * area * lambdad3;


434:   *f = area * (0.5 * fquad - lambdad3 * fexp);

436:   return 0;
437: } /* CombLocalFunctionGradient */


440: /*------- USER-DEFINED: routine to evaluate the Hessian
441:            at a local (rectangular element) level       -------*/

445: /*
446:   CombLocalHessian - Computes the Hessian of the local (partial) function
447:          defined over the current rectangle

449:   Input:
450:     coor - vector with the indices of the position of current element
451:              in the first, second and third directions
452:     x - current local solution (over the rectangle only)
453:     ptr - user-defined application context

455:   Output:
456:     H - Hessian matrix of the local function (wrt the four
457:            points of the rectangle only)

459: */
460: static int CombLocalHessian(PetscInt coor[2], double x[4], double H[4][4],void *ptr) {

462:   AppCtx *user = (AppCtx*)ptr;
463:   double hx, hy, lambdad3, area, dxdy, dydx;
464:   double diagxy, dexp, bandxy, bandyx;

466:   hx = user->hx;
467:   hy = user->hy;
468:   lambdad3 = user->param / 3.0;
469:   area = user->area;
470:   dxdy = hx/hy;
471:   dydx = hy/hx;
472:   diagxy = 0.5 * (dxdy + dydx);
473:   bandxy = -0.5 * dxdy;
474:   bandyx = -0.5 * dydx;

476:           /* Hessian contribution at 0,0 */
477:   dexp = exp(x[0]) * area * lambdad3;
478:   H[0][0] = diagxy - dexp;
479:   H[0][1] = H[1][0] = bandyx;
480:   H[0][2] = H[2][0] = bandxy;
481:   H[0][3] = H[3][0] = 0.0;

483:           /* Hessian contribution at 1,0 */
484:   dexp = exp(x[1]) * area * 2.0 * lambdad3;
485:   H[1][1] = diagxy - dexp;
486:   H[1][2] = H[2][1] = 0.0;
487:   H[1][3] = H[3][1] = bandxy;

489:           /* Hessian contribution at 0,1 */
490:   dexp = exp(x[2]) * area * 2.0 * lambdad3;
491:   H[2][2] = diagxy - dexp;
492:   H[2][3] = H[3][2] = bandyx;

494:           /* Hessian contribution at 1,1 */
495:   dexp = exp(x[3]) * area * lambdad3;
496:   H[3][3] = diagxy - dexp;

498:   return 0;
499: } /* CombLocalHessian */


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

507: /*
508:   WholeCombFunctionGradient - Evaluates function and gradient over the 
509:       whole grid

511:   Input:
512:     daapplication - TAO application object
513:     da  - distributed array
514:     X   - the current point, at which the function and gradient are evaluated
515:     ptr - user-defined application context

517:   Output:
518:     f - value of the objective funtion at X
519:     G - gradient at X
520: */
521: static int WholeCombFunctionGradient(TAO_APPLICATION daapplication, DA da, Vec X, double *f, Vec G, void *ptr) {

523:   AppCtx *user = (AppCtx*)ptr;
524:   Vec localX, localG;
525:   int info;
526:   PetscInt  i, j;
527:   PetscInt xs, xm, gxs, gxm, xe, ys, ym, gys, gym, ye;
528:   double **x, **g;
529:   double floc = 0.0;
530:   PetscScalar zero = 0.0;

532:   double lambdad3, hx, hy, area;
533:   double fquad, fexp, dvdx, dvdy;

535:   lambdad3 = user->param / 3.0;
536:   hx = user->hx;
537:   hy = user->hy;
538:   area = user->area;

540:   info = DAGetLocalVector(da, &localX); CHKERRQ(info);
541:   info = DAGetLocalVector(da, &localG); CHKERRQ(info);
542:   info = VecSet(G, zero); CHKERRQ(info);
543:   info = VecSet(localG, zero); CHKERRQ(info);

545:   info = DAGlobalToLocalBegin(da, X, INSERT_VALUES, localX); CHKERRQ(info);
546:   info = DAGlobalToLocalEnd(da, X, INSERT_VALUES, localX); CHKERRQ(info);

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

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

554:   xe = gxs + gxm - 1;
555:   ye = gys + gym - 1;
556:   for (j = ys; j < ye; j++) {
557:     for (i = xs; i < xe; i++) {

559:       /* lower triangle contribution */
560:       dvdx = (x[j][i] - x[j][i+1]) / hx;  
561:       dvdy = (x[j][i] - x[j+1][i]) / hy;
562:       fquad = dvdx * dvdx + dvdy * dvdy;
563:       fexp = exp(x[j][i]) + exp(x[j][i+1]) + exp(x[j+1][i]);

565:       dvdx = 0.5 * dvdx * hy;
566:       dvdy = 0.5 * dvdy * hx;
567:       g[j][i] += dvdx + dvdy - exp(x[j][i]) * area * lambdad3;
568:       g[j][i+1] += -dvdx - 2.0 * exp(x[j][i+1]) * area * lambdad3;
569:       g[j+1][i] += -dvdy - 2.0 * exp(x[j+1][i]) * area * lambdad3;

571:       /* upper triangle contribution */
572:       dvdx = (x[j+1][i+1] - x[j+1][i]) / hx;
573:       dvdy = (x[j+1][i+1] - x[j][i+1]) / hy;
574:       fquad += dvdx * dvdx + dvdy * dvdy;
575:       fexp += exp(x[j][i+1]) + exp(x[j+1][i]) + exp(x[j+1][i+1]);

577:       dvdx = 0.5 * dvdx * hy;
578:       dvdy = 0.5 * dvdy * hx;
579:       g[j][i+1] += -dvdy;
580:       g[j+1][i] += -dvdx;
581:       g[j+1][i+1] += dvdx + dvdy - exp(x[j+1][i+1]) * area * lambdad3;

583:       floc += area * (0.5 * fquad - fexp * lambdad3);

585:     }
586:   }

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

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

593:   info = DALocalToGlobalBegin(da, localG, G); CHKERRQ(info);
594:   info = DALocalToGlobalEnd(da, localG, G); CHKERRQ(info);

596:   info = DARestoreLocalVector(da, &localX); CHKERRQ(info);
597:   info = DARestoreLocalVector(da, &localG); CHKERRQ(info);

599:   info = PetscLogFlops((xe-xs) * (ye-ys) * 55 + 1); CHKERRQ(info);
600:   return 0;

602: } /* WholeCombFunctionGradient */


605: /*------- USER-DEFINED: routine to evaluate the Hessian 
606:           at the whole grid             -------*/
609: /*
610:   WholeCombHessian - Evaluates Hessian over the whole grid

612:   Input:
613:     daapplication - TAO application object
614:     da  - distributed array
615:     X   - the current point, at which the function and gradient are evaluated
616:     ptr - user-defined application context

618:   Output:
619:     H - Hessian at X
620: */
621: static int WholeCombHessian(TAO_APPLICATION daapplication, DA da, Vec X, Mat H, void *ptr) {

623:   AppCtx *user = (AppCtx*)ptr;
624:   Vec localX;
625:   int info;
626:   PetscInt i, j, ind[4];
627:   PetscInt xs, xm, gxs, gxm, xe, ys, ym, gys, gym, ye;
628:   double smallH[4][4];
629:   double **x;

631:   double hx, hy, lambdad3, area, dxdy, dydx;
632:   double diagxy, dexp, bandxy, bandyx;
633:   PetscTruth assembled;


636:   hx = user->hx;
637:   hy = user->hy;
638:   lambdad3 = user->param / 3.0;
639:   area = user->area;
640:   dxdy = hx/hy;
641:   dydx = hy/hx;
642:   diagxy = 0.5 * (dxdy + dydx);
643:   bandxy = -0.5 * dxdy;
644:   bandyx = -0.5 * dydx;

646:   info = DAGetLocalVector(da, &localX); CHKERRQ(info);
647:   info = MatAssembled(H,&assembled); CHKERRQ(info);
648:   if (assembled){info = MatZeroEntries(H);  CHKERRQ(info);}

650:   info = DAGlobalToLocalBegin(da, X, INSERT_VALUES, localX); CHKERRQ(info);
651:   info = DAGlobalToLocalEnd(da, X, INSERT_VALUES, localX); CHKERRQ(info);

653:   info = DAVecGetArray(da, localX, (void**)&x); CHKERRQ(info);

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

658:   xe = gxs + gxm - 1;
659:   ye = gys + gym - 1;
660:   for (j = ys; j < ye; j++) {
661:     for (i = xs; i < xe; i++) {

663:           /* Hessian contribution at 0,0 */
664:       dexp = exp(x[j][i]) * area * lambdad3;
665:       smallH[0][0] = diagxy - dexp;
666:       smallH[0][1] = smallH[1][0] = bandyx;
667:       smallH[0][2] = smallH[2][0] = bandxy;
668:       smallH[0][3] = smallH[3][0] = 0.0;

670:           /* Hessian contribution at 1,0 */
671:       dexp = exp(x[j][i+1]) * area * 2.0 * lambdad3;
672:       smallH[1][1] = diagxy - dexp;
673:       smallH[1][2] = smallH[2][1] = 0.0;
674:       smallH[1][3] = smallH[3][1] = bandxy;

676:           /* Hessian contribution at 0,1 */
677:       dexp = exp(x[j+1][i]) * area * 2.0 * lambdad3;
678:       smallH[2][2] = diagxy - dexp;
679:       smallH[2][3] = smallH[3][2] = bandyx;

681:           /* Hessian contribution at 1,1 */
682:       dexp = exp(x[j+1][i+1]) * area * lambdad3;
683:       smallH[3][3] = diagxy - dexp;

685:       ind[0] = (j-gys) * gxm + (i-gxs);
686:       ind[1] = ind[0] + 1;
687:       ind[2] = ind[0] + gxm;
688:       ind[3] = ind[2] + 1;
689:       info = MatSetValuesLocal(H,4,ind,4,ind,(PetscScalar*)smallH,ADD_VALUES); CHKERRQ(info);

691:     }
692:   }

694:   info = DAVecRestoreArray(da, localX, (void**)&x); CHKERRQ(info);

696:   info = MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY); CHKERRQ(info);
697:   info = MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY); CHKERRQ(info);
698:   info = MatSetOption(H, MAT_SYMMETRIC, PETSC_TRUE); CHKERRQ(info);

700:   info = DARestoreLocalVector(da, &localX); CHKERRQ(info);

702:   info = PetscLogFlops((xe-xs) * (ye-ys) * 14 + 7); CHKERRQ(info);
703:   return 0;

705: } /* WholeCombHessian */