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(int[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(int[2], double x[4], double *f, double g[4], void *ptr);
     70: static int EPTorsLocalHessian(int[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, int, void *);

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

     84:   int             info;                           /* used to check for functions returning nonzeros */
     85:   int             mx,my,Nx,Ny;
     86:   double          ff,gnorm;
     87:   int             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;

     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:   mx = my = 11;                               /* these correspond to 10 segments on each dimension */
    108:   info = PetscOptionsGetInt(TAO_NULL, "-mx", &mx, &flg); CHKERRQ(info);
    109:   info = PetscOptionsGetInt(TAO_NULL, "-my", &my, &flg); CHKERRQ(info);
    110:   if (PreLoadIt == 0) {
    111:     nlevels = 1; mx = 11; my = 11; }

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

    115:   /* Let PETSc determine the vector distribution */
    116:   Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;

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

    125:   /* Create TAO solver and set desired solution method */
    126:   info = TaoCreate(MPI_COMM_WORLD,method,&tao); CHKERRQ(info);
    127:   info = TaoApplicationCreate(PETSC_COMM_WORLD,&EPTorsApp); CHKERRQ(info);
    128:   info = TaoAppSetDAApp(EPTorsApp, DAarray, nlevels ); CHKERRQ(info);
    129:   /* Sets routines for function, gradient and bounds evaluation */
    130:   info = DAAppSetVariableBoundsRoutine(EPTorsApp,DASetBounds,(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(EPTorsApp,ad_EPTorsLocalFunction,192,(void *)&user.fgctx); CHKERRQ(info);
    139:     } else {
    140:       info = DAAppSetElementObjectiveAndGradientRoutine(EPTorsApp,EPTorsLocalFunctionGradient,42,(void *)&user); CHKERRQ(info);
    141:     }
    142:     /* Sets routines for Hessian evaluation, element by element */
    143:     info = DAAppSetElementHessianRoutine(EPTorsApp,EPTorsLocalHessian,6,(void*)&user); CHKERRQ(info);

    145:   } else {

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

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

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

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

    168:   info = DAAppGetSolution(EPTorsApp,0,&X); CHKERRQ(info);
    169:   info = FormInitialGuess(DAarray[0],X); CHKERRQ(info);
    170:   info = DAAppSetInitialSolution(EPTorsApp,X); CHKERRQ(info);
    171:
    172:   /* SOLVE THE APPLICATION */
    173:   info = TaoDAAppSolve(EPTorsApp, tao);  CHKERRQ(info);

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

    182:   info = PetscOptionsHasName(PETSC_NULL,"-view_sol",&flg); CHKERRQ(info);
    183:   if (flg){
    184:     info = DAAppGetSolution(EPTorsApp,nlevels-1,&X); CHKERRQ(info);
    185:     info=VecView(X,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(info);
    186:   }

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

    191:   /* Free TAO data structures */
    192:   info = TaoDestroy(tao); CHKERRQ(info);
    193:   info = TaoAppDestroy(EPTorsApp); CHKERRQ(info);

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

    200:   PreLoadEnd();

    202:   /* Finalize TAO */
    203:   TaoFinalize();
    204:   PetscFinalize();

    206:   return 0;
    207: } /* main */

    211: /*----- The following two routines
    212:   MyGridMonitorBefore    MyGridMonitorAfter
    213:   help diplay info of iterations at every grid level
    214: */
    217: static int MyGridMonitorBefore(TAO_APPLICATION myapp, DA da, int level, void *ctx) {

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

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

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

    234:   return 0;
    235: }

    238: /*------- USER-DEFINED: initialize the application context information -------*/
    241: /*
    242:   AppCtxInitialize - Sets initial values for the application context parameters

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

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

    253:   AppCtx *user = (AppCtx*)ptr;
    254:   PetscTruth    flg;            /* flag for PETSc calls */
    255:   int info;

    257:   /* Specify default parameters */
    258:   user->param = 25.0;
    259:   user->u1 = user->u2 = 1.0;

    261:   /* Check for command line arguments that override defaults */
    262:   info = PetscOptionsGetReal(TAO_NULL, "-par", &user->param, &flg); CHKERRQ(info);
    263:   info = PetscOptionsGetReal(TAO_NULL, "-u1", &user->u1, &flg); CHKERRQ(info);
    264:   info = PetscOptionsGetReal(TAO_NULL, "-u2", &user->u2, &flg); CHKERRQ(info);

    266:   return 0;
    267: } /* AppCtxInitialize */

    271: static int FormInitialGuess(DA da, Vec X)
    272: {
    273:   int    info, i, j, mx, my;
    274:   int    xs, ys, xm, ym, xe, ye;
    275:   PetscReal hx, hy, temp, val;
    276:   double **x;

    278:   /* Get local mesh boundaries */
    279:   info = DAGetInfo(da,PETSC_NULL,&mx,&my,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
    280:                    PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(info);
    281:   hx = 1.0/(mx-1);  hy = 1.0/(my-1);

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

    286:   info = DAVecGetArray(da, X, (void**)&x); CHKERRQ(info);
    287:   /* Compute initial guess over locally owned part of mesh */
    288:   for (j=ys; j<ye; j++) {  /*  for (j=0; j<my; j++) */
    289:     temp = PetscMin(j+1,my-j)*hy;
    290:     for (i=xs; i<xe; i++) {  /*  for (i=0; i<mx; i++) */
    291:       val = PetscMin((PetscMin(i+1,mx-i))*hx,temp);
    292:       x[j][i] = val;
    293:     }
    294:   }
    295:   info = DAVecRestoreArray(da, X, (void**)&x); CHKERRQ(info);

    297:   return 0;
    298: }

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

    305: /*
    306:   FormBounds - Forms bounds on the variables

    308:   Input:
    309:     user - user-defined application context

    311:   Output:
    312:     XL - vector of lower bounds
    313:     XU - vector of upper bounds
    314: */
    315: static int DASetBounds(TAO_APPLICATION daapplication, DA da, Vec XL, Vec XU, void *ptr)
    316: {
    317:   AppCtx *user = (AppCtx*)ptr;
    318:   int i, j, info, xs, xm, ys, ym;
    319:   double hx, hy, u1, u2, dist, d1, d2, hd, vd;
    320:   double **xl, **xu;

    322:   hx = user->hx;
    323:   hy = user->hy;
    324:   u1 = user->u1;
    325:   u2 = user->u2;

    327:   info = DAVecGetArray(da, XL, (void**)&xl); CHKERRQ(info);
    328:   info = DAVecGetArray(da, XU, (void**)&xu); CHKERRQ(info);
    329:   info = DAGetCorners(da, &xs, &ys, TAO_NULL, &xm, &ym, TAO_NULL); CHKERRQ(info);

    331:   for (j = ys; j < ys+ym; j++){
    332:     for (i = xs; i < xs+xm; i++){
    333:       d1 = i * hx; d2 = u1 - d1; hd = PetscMin(d1,d2);
    334:       d1 = j * hy; d2 = u2 - d1; vd = PetscMin(d1,d2);
    335:       dist = PetscMin(hd,vd);
    336:       xl[j][i] = -dist;
    337:       xu[j][i] = dist;
    338:     }
    339:   }

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

    344:   info = PetscLogFlops(xm * ym * 4); CHKERRQ(info);
    345:   return 0;

    347: } /* DASetBounds */

    352: /*
    353:   EPTorsLocalFunctionGradient - Evaluates function and gradient over the
    354:       local rectangular element

    356:   Input:
    357:     coor - vector with the indices of the position of current element
    358:              in the first, second and third directions
    359:     x - current point (values over the current rectangular element)
    360:     df - degrees of freedom at each point
    361:     ptr - user-defined application context

    363:   Output:
    364:     f - value of the objective funtion at the local rectangular element
    365:     g - gradient of the local function
    366: */
    367: static int EPTorsLocalFunctionGradient(int coor[2], double x[4], double *f, double g[4], void *ptr) {

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

    371:   double fquad, flin;
    372:   double hx, hy, dvdx, dvdy, area;
    373:   double cdiv3, cnt;

    375:   cdiv3 = user->param / 3.0;
    376:   hx = user->hx;
    377:   hy = user->hy;
    378:   area = user->area;
    379:   cnt = area * cdiv3;

    381:   /* lower triangle contribution */
    382:   dvdx = (x[0] - x[1]) / hx;
    383:   dvdy = (x[0] - x[2]) / hy;
    384:   fquad = dvdx * dvdx + dvdy * dvdy;
    385:   flin = x[0] + x[1] + x[2];

    387:   dvdx = 0.5 * dvdx * hy;
    388:   dvdy = 0.5 * dvdy * hx;
    389:   g[0] = dvdx + dvdy - cnt;
    390:   g[1] = -dvdx - 2.0 * cnt;
    391:   g[2] = -dvdy - 2.0 * cnt;

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

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

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

    407:   return 0;
    408: } /* EPTorsLocalFunctionGradient */

    412: /*------- USER-DEFINED: routine to evaluate the Hessian
    413:            at a local (rectangular element) level       -------*/
    416: /*
    417:   EPTorsLocalHessian - Computes the Hessian of the local (partial) function
    418:          defined over the current rectangle

    420:   Input:
    421:     coor - vector with the indices of the position of current element
    422:              in the first, second and third directions
    423:     x - current local solution (over the rectangle only)
    424:     df - degrees of freedom at each point
    425:     ptr - user-defined application context

    427:   Output:
    428:     H - Hessian matrix of the local function (wrt the four
    429:            points of the rectangle only)
    430: */
    431: static int EPTorsLocalHessian(int coor[2], double x[4], double H[4][4], void *ptr) {

    433:   AppCtx *user = (AppCtx*)ptr;
    434:   double hx, hy, dxdy, dydx;
    435:   double diagxy, bandxy, bandyx;

    437:   hx = user->hx;
    438:   hy = user->hy;
    439:   dxdy = hx/hy;
    440:   dydx = hy/hx;
    441:   diagxy = 0.5 * (dxdy + dydx);
    442:   bandxy = -0.5 * dxdy;
    443:   bandyx = -0.5 * dydx;

    445:           /* Hessian contribution at 0,0 */
    446:   H[0][0] = diagxy;
    447:   H[0][1] =  H[1][0] = bandyx;
    448:   H[0][2] =  H[2][0] = bandxy;
    449:   H[0][3] =  H[3][0] = 0.0;

    451:           /* Hessian contribution at 1,0 */
    452:   H[1][1] = diagxy;
    453:   H[1][2] =  H[2][1] = 0.0;
    454:   H[1][3] =  H[3][1] = bandxy;

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

    460:           /* Hessian contribution at 1,1 */
    461:   H[3][3] = diagxy;

    463:   return 0;

    465: } /* EPTorsLocalHessian */

    468: /*------- USER-DEFINED: routine to evaluate the function
    469:           and gradient at the whole grid             -------*/
    472: /*
    473:   WholeEPTorsFunctionGradient - Evaluates function and gradient over the
    474:       whole grid

    476:   Input:
    477:     daapplication - TAO application object
    478:     da  - distributed array
    479:     X   - the current point, at which the function and gradient are evaluated
    480:     ptr - user-defined application context

    482:   Output:
    483:     f - value of the objective funtion at X
    484:     G - gradient at X
    485: */
    486: static int WholeEPTorsFunctionGradient(TAO_APPLICATION daapplication, DA da, Vec X, double *f, Vec G, void *ptr) {

    488:   AppCtx *user = (AppCtx*)ptr;
    489:   Vec localX, localG;
    490:   int info, i, j;
    491:   int xs, xm, gxs, gxm, xe, ys, ym, gys, gym, ye;
    492:   double **x, **g;
    493:   double floc = 0.0;
    494:   PetscScalar zero = 0.0;

    496:   double fquad, flin;
    497:   double hx, hy, dvdx, dvdy, area;
    498:   double cdiv3, cnt;

    500:   cdiv3 = user->param / 3.0;
    501:   hx = user->hx;
    502:   hy = user->hy;
    503:   area = user->area;
    504:   cnt = area * cdiv3;

    506:   info = DAGetLocalVector(da, &localX); CHKERRQ(info);
    507:   info = DAGetLocalVector(da, &localG); CHKERRQ(info);
    508:   info = VecSet(G, zero); CHKERRQ(info);
    509:   info = VecSet(localG, zero); CHKERRQ(info);

    511:   info = DAGlobalToLocalBegin(da, X, INSERT_VALUES, localX); CHKERRQ(info);
    512:   info = DAGlobalToLocalEnd(da, X, INSERT_VALUES, localX); CHKERRQ(info);

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

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

    520:   xe = gxs + gxm - 1;
    521:   ye = gys + gym - 1;
    522:   for (j = ys; j < ye; j++) {
    523:     for (i = xs; i < xe; i++) {

    525:       /* lower triangle contribution */
    526:       dvdx = (x[j][i] - x[j][i+1]) / hx;
    527:       dvdy = (x[j][i] - x[j+1][i]) / hy;
    528:       fquad = dvdx * dvdx + dvdy * dvdy;
    529:       flin = x[j][i] + x[j][i+1] + x[j+1][i];

    531:       dvdx = 0.5 * dvdx * hy;
    532:       dvdy = 0.5 * dvdy * hx;
    533:       g[j][i] += dvdx + dvdy - cnt;
    534:       g[j][i+1] += -dvdx - 2.0 * cnt;
    535:       g[j+1][i] += -dvdy - 2.0 * cnt;

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

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

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

    551:     }
    552:   }

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

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

    559:   info = DALocalToGlobalBegin(da, localG, G); CHKERRQ(info);
    560:   info = DALocalToGlobalEnd(da, localG, G); CHKERRQ(info);

    562:   info = DARestoreLocalVector(da, &localX); CHKERRQ(info);
    563:   info = DARestoreLocalVector(da, &localG); CHKERRQ(info);

    565:   info = PetscLogFlops((xe-xs) * (ye-ys) * 47 + 2); CHKERRQ(info);
    566:   return 0;
    567: } /* WholeEPTorsFunctionGradient */

    570: /*------- USER-DEFINED: routine to evaluate the Hessian
    571:           at the whole grid             -------*/
    574: /*
    575:   WholeEPTorsHessian - Evaluates Hessian over the whole grid

    577:   Input:
    578:     daapplication - TAO application object
    579:     da  - distributed array
    580:     X   - the current point, at which the function and gradient are evaluated
    581:     ptr - user-defined application context

    583:   Output:
    584:     H - Hessian at X
    585: */
    586: static int WholeEPTorsHessian(TAO_APPLICATION daapplication, DA da, Vec X, Mat H, void *ptr) {

    588:   AppCtx *user = (AppCtx*)ptr;
    589:   int info, i, j, ind[4];
    590:   int xs, xm, gxs, gxm, xe, ys, ym, gys, gym, ye;
    591:   double smallH[4][4];

    593:   double hx, hy, dxdy, dydx;
    594:   double diagxy, bandxy, bandyx;
    595:   PetscTruth assembled;

    597:   hx = user->hx;
    598:   hy = user->hy;
    599:   dxdy = hx/hy;
    600:   dydx = hy/hx;
    601:   diagxy = 0.5 * (dxdy + dydx);
    602:   bandxy = -0.5 * dxdy;
    603:   bandyx = -0.5 * dydx;

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

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

    612:   xe = gxs + gxm - 1;
    613:   ye = gys + gym - 1;
    614:   for (j = ys; j < ye; j++) {
    615:     for (i = xs; i < xe; i++) {

    617:           /* Hessian contribution at 0,0 */
    618:       smallH[0][0] = diagxy;
    619:       smallH[0][1] = smallH[1][0] = bandyx;
    620:       smallH[0][2] = smallH[2][0] = bandxy;
    621:       smallH[0][3] = smallH[3][0] = 0.0;

    623:           /* Hessian contribution at 1,0 */
    624:       smallH[1][1] = diagxy;
    625:       smallH[1][2] = smallH[2][1] = 0.0;
    626:       smallH[1][3] = smallH[3][1] = bandxy;

    628:           /* Hessian contribution at 0,1 */
    629:       smallH[2][2] = diagxy;
    630:       smallH[2][3] = smallH[3][2] = bandyx;

    632:           /* Hessian contribution at 1,1 */
    633:       smallH[3][3] = diagxy;

    635:       ind[0] = (j-gys) * gxm + (i-gxs);
    636:       ind[1] = ind[0] + 1;
    637:       ind[2] = ind[0] + gxm;
    638:       ind[3] = ind[2] + 1;
    639:       info = MatSetValuesLocal(H,4,ind,4,ind,(PetscScalar*)smallH,ADD_VALUES); CHKERRQ(info);

    641:     }
    642:   }

    644:   info = MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY); CHKERRQ(info);
    645:   info = MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY); CHKERRQ(info);
    646:   info = MatSetOption(H, MAT_SYMMETRIC); CHKERRQ(info);

    649:   info = PetscLogFlops(6); CHKERRQ(info);
    650:   return 0;

    652: } /* WholeEPTorsHessian */