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:   int     mx, my;         /* discretization including boundaries */
     61:   ADFGCtx fgctx;          /* Used only when an ADIC generated gradient is used */
     62: } AppCtx;
     63: int ad_CombLocalFunction(int[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(int[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(int[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, int, void *);

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

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

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

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

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

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

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

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

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

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

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

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

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

    144:   } else {

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

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

    154:   info = DAAppSetBeforeMonitor(CombApp,MyGridMonitorBefore,(void*)&user); CHKERRQ(info);
    155:   info = DAAppPrintStageTimes(CombApp); CHKERRQ(info);
    156:   info = DAAppPrintInterpolationError(CombApp); CHKERRQ(info);

    158:   info = TaoAppSetRelativeTolerance(CombApp,1.0e-8); CHKERRQ(info);
    159:   info = TaoSetTolerances(tao,0,0,0,0); CHKERRQ(info);
    160:   info = TaoSetGradientTolerances(tao,0,0,0); CHKERRQ(info);

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

    165:   info = DAAppGetSolution(CombApp,0,&X); CHKERRQ(info);
    166:   info = FormInitialGuess(DAarray[0],X,&user); CHKERRQ(info);
    167:   info = DAAppSetInitialSolution(CombApp,X); CHKERRQ(info);

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

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

    179:   info = PetscOptionsHasName(PETSC_NULL,"-view_sol",&flg); CHKERRQ(info);
    180:   if (flg){
    181:     info = DAAppGetSolution(CombApp,nlevels-1,&X); CHKERRQ(info);
    182:     info=VecView(X,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(info);
    183:   }

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

    188:   /* Free TAO data structures */
    189:   info = TaoDestroy(tao); CHKERRQ(info);
    190:   info = TaoAppDestroy(CombApp); CHKERRQ(info);

    192:   /* Free PETSc data structures */
    193:   for (iter=0;iter<nlevels;iter++){
    194:     info = DADestroy(DAarray[iter]); CHKERRQ(info);
    195:   }

    197:   PreLoadEnd();

    199:   /* Finalize TAO */
    200:   TaoFinalize();
    201:   PetscFinalize();

    203:   return 0;
    204: } /* main */

    208: /*----- The following two routines
    209:
    210:   MyGridMonitorBefore    MyGridMonitorAfter

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

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

    218:   AppCtx *user = (AppCtx*)ctx;
    219:   int info,mx,my;

    221:   info = DAGetInfo(da,PETSC_NULL,&mx,&my,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
    222:                    PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(info);
    223:   user->mx = mx;
    224:   user->my = my;
    225:   user->hx = 1.0 / (user->mx - 1);
    226:   user->hy = 1.0 / (user->my - 1);
    227:   user->area = 0.5 * user->hx * user->hy;
    228:   user->fgctx.hx   = user->hx;
    229:   user->fgctx.hy   = user->hy;
    230:   user->fgctx.area = user->area;
    231:   user->fgctx.param = user->param;

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

    235:   return 0;
    236: }

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

    242: /*
    243:   AppCtxInitialize - Sets initial values for the application context parameters

    245:   Input:
    246:     ptr - void user-defined application context

    248:   Output:
    249:     ptr - user-defined application context with the default or user-provided
    250:              parameters
    251: */
    252: static int AppCtxInitialize(void *ptr) {

    254:   AppCtx *user = (AppCtx*)ptr;
    255:   PetscReal     LambdaMax = 6.81, LambdaMin = 0.0;  /* bounds on parameter lambda */
    256:   PetscTruth    flg;            /* flag for PETSc calls */
    257:   int info;

    259:   /* Specify dimension of the problem */
    260:   user->param = 5.0;
    261:   user->mx = 11;
    262:   user->my = 11;

    264:   /* Check for any command line arguments that override defaults */
    265:   info = PetscOptionsGetReal(TAO_NULL, "-par", &user->param, &flg); CHKERRQ(info);
    266:   if (user->param >= LambdaMax || user->param <= LambdaMin) {
    267:     SETERRQ(1,"Lambda is out of range.");
    268:   }
    269:   info = PetscOptionsGetInt(PETSC_NULL,"-mx",&user->mx,&flg); CHKERRQ(info);
    270:   info = PetscOptionsGetInt(PETSC_NULL,"-my",&user->my,&flg); CHKERRQ(info);

    272:   user->hx = 1.0 / (user->mx - 1);
    273:   user->hy = 1.0 / (user->my - 1);
    274:   user->area = 0.5 * user->hx * user->hy;
    275:   info = PetscLogFlops(6); CHKERRQ(info);

    277:   return 0;
    278: } /* AppCtxInitialize */

    284: static int FormInitialGuess(DA da, Vec X, AppCtx *ctx)
    285: {
    286:   int    info, i, j, mx, my;
    287:   int    xs, ys, xm, ym, xe, ye;
    288:   PetscReal hx, hy, temp, val, lambda;
    289:   double **x;

    291:   lambda = ctx->param;
    292:   lambda = lambda/(lambda+1.0);

    294:   /* Get local mesh boundaries */
    295:   info = DAGetInfo(da,PETSC_NULL,&mx,&my,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
    296:                    PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(info);
    297:   hx = 1.0/(mx-1);  hy = 1.0/(my-1);

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

    302:   info = DAVecGetArray(da, X, (void**)&x); CHKERRQ(info);
    303:   /* Compute initial guess over locally owned part of mesh */
    304:   for (j=ys; j<ye; j++) {  /*  for (j=0; j<my; j++) */
    305:     temp = PetscMin(j+1,my-j)*hy;
    306:     for (i=xs; i<xe; i++) {  /*  for (i=0; i<mx; i++) */
    307:       val = lambda*sqrt(PetscMin((PetscMin(i+1,mx-i))*hx,temp));
    308:       x[j][i] = val;
    309:     }
    310:   }
    311:   info = DAVecRestoreArray(da, X, (void**)&x); CHKERRQ(info);

    313:   return 0;
    314: }

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

    321: /*
    322:   FormBounds - Forms bounds on the variables

    324:   Input:
    325:     user - user-defined application context

    327:   Output:
    328:     XL - vector of lower bounds
    329:     XU - vector of upper bounds

    331: */
    332: static int DAFixBoundary(TAO_APPLICATION daapplication, DA da, Vec XL, Vec XU, void *ptr)
    333: {
    334:   AppCtx *user = (AppCtx*)ptr;
    335:   int i, j, mx, my, info, xs, xm, ys, ym;
    336:   double lb = -TAO_INFINITY;
    337:   double ub = TAO_INFINITY;
    338:   double **xl, **xu;

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

    343:   info = DAVecGetArray(da, XL, (void**)&xl); CHKERRQ(info);
    344:   info = DAVecGetArray(da, XU, (void**)&xu); CHKERRQ(info);
    345:   info = DAGetCorners(da, &xs, &ys, TAO_NULL, &xm, &ym, TAO_NULL); CHKERRQ(info);

    347:   /* Compute initial guess over locally owned part of the grid */
    348:   for (j = ys; j < ys+ym; j++){
    349:     for (i = xs; i < xs+xm; i++){
    350:       if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) {
    351:         xl[j][i] = xu[j][i] = 0.0;
    352:       } else {
    353:         xl[j][i] = lb;
    354:         xu[j][i] = ub;
    355:       }
    356:     }
    357:   }

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

    362:   return 0;
    363: } /* DAFixBoundary */

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

    371: /*
    372:   CombLocalFunctionGradient - Evaluates function and gradient over the
    373:       local rectangular element

    375:   Input:
    376:     coor - vector with the indices of the position of current element
    377:              in the first, second and third directions
    378:     x - current point (values over the current rectangular element)
    379:     ptr - user-defined application context

    381:   Output:
    382:     f - value of the objective funtion at the local rectangular element
    383:     g - gradient of the local function

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

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

    390:   double lambdad3, hx, hy, area;
    391:   double fquad, fexp, dvdx, dvdy;

    393:   lambdad3 = user->param / 3.0;
    394:   hx = user->hx;
    395:   hy = user->hy;
    396:   area = user->area;

    398:   /* lower triangle contribution */
    399:   dvdx = (x[0] - x[1]) / hx;
    400:   dvdy = (x[0] - x[2]) / hy;
    401:   fquad = dvdx * dvdx + dvdy * dvdy;
    402:   fexp = exp(x[0]) + exp(x[1]) + exp(x[2]);

    404:   dvdx = 0.5 * dvdx * hy;
    405:   dvdy = 0.5 * dvdy * hx;
    406:   g[0] = dvdx + dvdy - exp(x[0]) * area * lambdad3;
    407:   g[1] = -dvdx - 2.0 * exp(x[1]) * area * lambdad3;
    408:   g[2] = -dvdy - 2.0 * exp(x[2]) * area * lambdad3;

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

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

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

    425:   return 0;
    426: } /* CombLocalFunctionGradient */

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

    434: /*
    435:   CombLocalHessian - Computes the Hessian of the local (partial) function
    436:          defined over the current rectangle

    438:   Input:
    439:     coor - vector with the indices of the position of current element
    440:              in the first, second and third directions
    441:     x - current local solution (over the rectangle only)
    442:     ptr - user-defined application context

    444:   Output:
    445:     H - Hessian matrix of the local function (wrt the four
    446:            points of the rectangle only)

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

    451:   AppCtx *user = (AppCtx*)ptr;
    452:   double hx, hy, lambdad3, area, dxdy, dydx;
    453:   double diagxy, dexp, bandxy, bandyx;

    455:   hx = user->hx;
    456:   hy = user->hy;
    457:   lambdad3 = user->param / 3.0;
    458:   area = user->area;
    459:   dxdy = hx/hy;
    460:   dydx = hy/hx;
    461:   diagxy = 0.5 * (dxdy + dydx);
    462:   bandxy = -0.5 * dxdy;
    463:   bandyx = -0.5 * dydx;

    465:           /* Hessian contribution at 0,0 */
    466:   dexp = exp(x[0]) * area * lambdad3;
    467:   H[0][0] = diagxy - dexp;
    468:   H[0][1] = H[1][0] = bandyx;
    469:   H[0][2] = H[2][0] = bandxy;
    470:   H[0][3] = H[3][0] = 0.0;

    472:           /* Hessian contribution at 1,0 */
    473:   dexp = exp(x[1]) * area * 2.0 * lambdad3;
    474:   H[1][1] = diagxy - dexp;
    475:   H[1][2] = H[2][1] = 0.0;
    476:   H[1][3] = H[3][1] = bandxy;

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

    483:           /* Hessian contribution at 1,1 */
    484:   dexp = exp(x[3]) * area * lambdad3;
    485:   H[3][3] = diagxy - dexp;

    487:   return 0;
    488: } /* CombLocalHessian */

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

    496: /*
    497:   WholeCombFunctionGradient - Evaluates function and gradient over the
    498:       whole grid

    500:   Input:
    501:     daapplication - TAO application object
    502:     da  - distributed array
    503:     X   - the current point, at which the function and gradient are evaluated
    504:     ptr - user-defined application context

    506:   Output:
    507:     f - value of the objective funtion at X
    508:     G - gradient at X
    509: */
    510: static int WholeCombFunctionGradient(TAO_APPLICATION daapplication, DA da, Vec X, double *f, Vec G, void *ptr) {

    512:   AppCtx *user = (AppCtx*)ptr;
    513:   Vec localX, localG;
    514:   int info, i, j;
    515:   int xs, xm, gxs, gxm, xe, ys, ym, gys, gym, ye;
    516:   double **x, **g;
    517:   double floc = 0.0;
    518:   PetscScalar zero = 0.0;

    520:   double lambdad3, hx, hy, area;
    521:   double fquad, fexp, dvdx, dvdy;

    523:   lambdad3 = user->param / 3.0;
    524:   hx = user->hx;
    525:   hy = user->hy;
    526:   area = user->area;

    528:   info = DAGetLocalVector(da, &localX); CHKERRQ(info);
    529:   info = DAGetLocalVector(da, &localG); CHKERRQ(info);
    530:   info = VecSet(G, zero); CHKERRQ(info);
    531:   info = VecSet(localG, zero); CHKERRQ(info);

    533:   info = DAGlobalToLocalBegin(da, X, INSERT_VALUES, localX); CHKERRQ(info);
    534:   info = DAGlobalToLocalEnd(da, X, INSERT_VALUES, localX); CHKERRQ(info);

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

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

    542:   xe = gxs + gxm - 1;
    543:   ye = gys + gym - 1;
    544:   for (j = ys; j < ye; j++) {
    545:     for (i = xs; i < xe; i++) {

    547:       /* lower triangle contribution */
    548:       dvdx = (x[j][i] - x[j][i+1]) / hx;
    549:       dvdy = (x[j][i] - x[j+1][i]) / hy;
    550:       fquad = dvdx * dvdx + dvdy * dvdy;
    551:       fexp = exp(x[j][i]) + exp(x[j][i+1]) + exp(x[j+1][i]);

    553:       dvdx = 0.5 * dvdx * hy;
    554:       dvdy = 0.5 * dvdy * hx;
    555:       g[j][i] += dvdx + dvdy - exp(x[j][i]) * area * lambdad3;
    556:       g[j][i+1] += -dvdx - 2.0 * exp(x[j][i+1]) * area * lambdad3;
    557:       g[j+1][i] += -dvdy - 2.0 * exp(x[j+1][i]) * area * lambdad3;

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

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

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

    573:     }
    574:   }

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

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

    581:   info = DALocalToGlobalBegin(da, localG, G); CHKERRQ(info);
    582:   info = DALocalToGlobalEnd(da, localG, G); CHKERRQ(info);

    584:   info = DARestoreLocalVector(da, &localX); CHKERRQ(info);
    585:   info = DARestoreLocalVector(da, &localG); CHKERRQ(info);

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

    590: } /* WholeCombFunctionGradient */

    593: /*------- USER-DEFINED: routine to evaluate the Hessian
    594:           at the whole grid             -------*/
    597: /*
    598:   WholeCombHessian - Evaluates Hessian over the whole grid

    600:   Input:
    601:     daapplication - TAO application object
    602:     da  - distributed array
    603:     X   - the current point, at which the function and gradient are evaluated
    604:     ptr - user-defined application context

    606:   Output:
    607:     H - Hessian at X
    608: */
    609: static int WholeCombHessian(TAO_APPLICATION daapplication, DA da, Vec X, Mat H, void *ptr) {

    611:   AppCtx *user = (AppCtx*)ptr;
    612:   Vec localX;
    613:   int info, i, j, ind[4];
    614:   int xs, xm, gxs, gxm, xe, ys, ym, gys, gym, ye;
    615:   double smallH[4][4];
    616:   double **x;

    618:   double hx, hy, lambdad3, area, dxdy, dydx;
    619:   double diagxy, dexp, bandxy, bandyx;
    620:   PetscTruth assembled;

    623:   hx = user->hx;
    624:   hy = user->hy;
    625:   lambdad3 = user->param / 3.0;
    626:   area = user->area;
    627:   dxdy = hx/hy;
    628:   dydx = hy/hx;
    629:   diagxy = 0.5 * (dxdy + dydx);
    630:   bandxy = -0.5 * dxdy;
    631:   bandyx = -0.5 * dydx;

    633:   info = DAGetLocalVector(da, &localX); CHKERRQ(info);
    634:   info = MatAssembled(H,&assembled); CHKERRQ(info);
    635:   if (assembled){info = MatZeroEntries(H);  CHKERRQ(info);}

    637:   info = DAGlobalToLocalBegin(da, X, INSERT_VALUES, localX); CHKERRQ(info);
    638:   info = DAGlobalToLocalEnd(da, X, INSERT_VALUES, localX); CHKERRQ(info);

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

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

    645:   xe = gxs + gxm - 1;
    646:   ye = gys + gym - 1;
    647:   for (j = ys; j < ye; j++) {
    648:     for (i = xs; i < xe; i++) {

    650:           /* Hessian contribution at 0,0 */
    651:       dexp = exp(x[j][i]) * area * lambdad3;
    652:       smallH[0][0] = diagxy - dexp;
    653:       smallH[0][1] = smallH[1][0] = bandyx;
    654:       smallH[0][2] = smallH[2][0] = bandxy;
    655:       smallH[0][3] = smallH[3][0] = 0.0;

    657:           /* Hessian contribution at 1,0 */
    658:       dexp = exp(x[j][i+1]) * area * 2.0 * lambdad3;
    659:       smallH[1][1] = diagxy - dexp;
    660:       smallH[1][2] = smallH[2][1] = 0.0;
    661:       smallH[1][3] = smallH[3][1] = bandxy;

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

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

    672:       ind[0] = (j-gys) * gxm + (i-gxs);
    673:       ind[1] = ind[0] + 1;
    674:       ind[2] = ind[0] + gxm;
    675:       ind[3] = ind[2] + 1;
    676:       info = MatSetValuesLocal(H,4,ind,4,ind,(PetscScalar*)smallH,ADD_VALUES); CHKERRQ(info);

    678:     }
    679:   }

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

    683:   info = MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY); CHKERRQ(info);
    684:   info = MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY); CHKERRQ(info);
    685:   info = MatSetOption(H, MAT_SYMMETRIC); CHKERRQ(info);

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

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

    692: } /* WholeCombHessian */