Actual source code: jbearing3.c

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

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

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

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

     15: static char help[] ="Pressure distribution in a Journal Bearing. \n\
     16: This example is based on the problem DPJB from the MINPACK-2 test suite.\n\
     17: This pressure journal bearing problem is an example of elliptic variational\n\
     18: problem defined over a two dimensional rectangle. By discretizing the domain \n\
     19: into triangular elements, the pressure surrounding the journal bearing is\n\
     20: defined as the minimum of a quadratic function whose variables are bounded\n\
     21: below by zero. The command line options are:\n\
     22:   -ecc <ecc>, where <ecc> = epsilon parameter\n\
     23:   -b <b>, where <b> = half the upper limit in the 2nd coordinate direction\n\
     24:   -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
     25:   -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
     26:   -nlevels <nlevels>, where <nlevels> = number of levels in multigrid\n\
     27:   -byelement, if computation is made by functions on rectangular elements\n\
     28:   -adic, if AD is used (AD is not used by default)\n\n";

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

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

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

     61: } ADFGCtx;

     64: typedef struct {

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

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

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

     76: } AppCtx;

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

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

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

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

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

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

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

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

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

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

    125:   /* Let PETSc determine the vector distribution */
    126:   Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;

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

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

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

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

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

    156:   } else {

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

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

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

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

    180:   info = DAAppGetSolution(JBearApp,0,&X); CHKERRQ(info);
    181:   info = FormInitialGuess(DAarray[0],X); CHKERRQ(info);
    182:   info = DAAppSetInitialSolution(JBearApp,X); CHKERRQ(info);
    183:   /* SOLVE THE APPLICATION */
    184:   info = TaoDAAppSolve(JBearApp,tao);  CHKERRQ(info);

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

    193:   info = PetscOptionsHasName(PETSC_NULL,"-view_sol",&flg); CHKERRQ(info);
    194:   if (flg){
    195:     info = DAAppGetSolution(JBearApp,nlevels-1,&X); CHKERRQ(info);
    196:     info=VecView(X,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(info);
    197:   }

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

    202:   /* Free TAO data structures */
    203:   info = TaoDestroy(tao); CHKERRQ(info);
    204:   info = TaoAppDestroy(JBearApp); CHKERRQ(info);

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

    211:   PreLoadEnd();

    213:   /* Finalize TAO */
    214:   TaoFinalize();
    215:   PetscFinalize();

    217:   return 0;
    218: } /* main */

    222: /*----- The following two routines
    223:
    224:   MyGridMonitorBefore    MyGridMonitorAfter

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

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

    232:   AppCtx *user = (AppCtx*)ctx;
    233:   int info,mx,my;
    234:   double t;

    236:   info = DAGetInfo(da,PETSC_NULL,&mx,&my,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
    237:                    PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(info);
    238:   user->mx = mx;
    239:   user->my = my;

    241:   user->hx = (2.0 * 3.14159265358979) / (user->mx - 1);
    242:   user->hy = (2.0 * user->b) / (user->my - 1);
    243:   user->area = 0.5 * user->hx * user->hy;

    245:   user->wq = new double[user->mx];
    246:   user->wl = new double[user->mx];
    247:   for (int i=0; i<user->mx; i++) {
    248:     t = 1.0 + user->ecc * cos(i*user->hx);
    249:     user->wq[i] = t*t*t;
    250:     user->wl[i] = user->ecc * sin(i*user->hx);
    251:   }
    252:   info = PetscLogFlops(8 + 7 * user->mx); CHKERRQ(info);

    254:   user->fgctx.hx   = user->hx;
    255:   user->fgctx.hy   = user->hy;
    256:   user->fgctx.area = user->area;
    257:   user->fgctx.wq = user->wq;
    258:   user->fgctx.wl = user->wl;

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

    262:   return 0;
    263: }

    267: static int MyGridMonitorAfter1(TAO_APPLICATION myapp, DA da, int level, void *ctx){
    268:
    269:   AppCtx *user = (AppCtx*)ctx;
    270:   delete [] user->wq; delete [] user->wl;
    271:   return 0;
    272: }

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

    280: /*
    281:   AppCtxInitialize - Sets initial values for the application context parameters

    283:   Input:
    284:     ptr - void user-defined application context

    286:   Output:
    287:     ptr - user-defined application context with the default or user-provided
    288:              parameters
    289: */
    290: static int AppCtxInitialize(void *ptr) {

    292:   AppCtx *user = (AppCtx*)ptr;
    293:   PetscTruth    flg;            /* flag for PETSc calls */
    294:   int info;

    296:   /* Specify default parameters */
    297:   user->ecc = 0.1;
    298:   user->b = 10.0;
    299:   user->mx = user->my = 11;

    301:   /* Check for command line arguments that override defaults */
    302:   info = PetscOptionsGetReal(TAO_NULL, "-ecc", &user->ecc, &flg); CHKERRQ(info);
    303:   info = PetscOptionsGetReal(TAO_NULL, "-b", &user->b, &flg); CHKERRQ(info);
    304:   info = PetscOptionsGetInt(TAO_NULL, "-mx", &user->mx, &flg); CHKERRQ(info);
    305:   info = PetscOptionsGetInt(TAO_NULL, "-my", &user->my, &flg); CHKERRQ(info);

    307:   return 0;
    308: } /* AppCtxInitialize */

    313: static int FormInitialGuess(DA da, Vec X)
    314: {
    315:   int    info, i, j, mx;
    316:   int    xs, ys, xm, ym, xe, ye;
    317:   PetscReal hx, val;
    318:   double **x;

    320:   /* Get local mesh boundaries */
    321:   info = DAGetInfo(da,PETSC_NULL,&mx,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
    322:                    PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(info);
    323:   hx = 2.0*4.0*atan(1.0)/(mx-1);

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

    328:   info = DAVecGetArray(da, X, (void**)&x); CHKERRQ(info);
    329:   /* Compute initial guess over locally owned part of mesh */
    330:   for (j=ys; j<ye; j++) {  /*  for (j=0; j<my; j++) */
    331:     for (i=xs; i<xe; i++) {  /*  for (i=0; i<mx; i++) */
    332:       val = PetscMax(sin((i+1.0)*hx),0.0);
    333:       x[j][i] = val;
    334:       x[j][i] = 0;
    335:     }
    336:   }
    337:   info = DAVecRestoreArray(da, X, (void**)&x); CHKERRQ(info);

    339:   return 0;
    340: }

    343: /*------- USER-DEFINED: set the upper and lower bounds for the variables  -------*/
    346: /*
    347:   FormBounds - Forms bounds on the variables

    349:   Input:
    350:     user - user-defined application context

    352:   Output:
    353:     XL - vector of lower bounds
    354:     XU - vector of upper bounds

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

    359:   AppCtx *user = (AppCtx*)ptr;
    360:   int i, j, info, mx, my;
    361:   int xs, xm, ys, ym;
    362:   double **xl, **xu;

    364:   mx = user->mx;
    365:   my = user->my;

    367:   info = DAVecGetArray(da, XL, (void**)&xl); CHKERRQ(info);
    368:   info = DAVecGetArray(da, XU, (void**)&xu); CHKERRQ(info);
    369:   info = DAGetCorners(da, &xs, &ys, TAO_NULL, &xm, &ym, TAO_NULL); CHKERRQ(info);

    371:   for (j = ys; j < ys+ym; j++){
    372:     for (i = xs; i < xs+xm; i++){
    373:       xl[j][i] = 0.0;
    374:       if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) {
    375:         xu[j][i] = 0.0;
    376:       } else {
    377:         xu[j][i] = TAO_INFINITY;
    378:       }
    379:     }
    380:   }

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

    385:   return 0;

    387: } /* DASetBounds */

    393: /*
    394:   JBearLocalFunctionGradient - Evaluates function and gradient over the
    395:       local rectangular element

    397:   Input:
    398:     coor - vector with the indices of the position of current element
    399:              in the first, second and third directions
    400:     x - current point (values over the current rectangular element)
    401:     ptr - user-defined application context

    403:   Output:
    404:     f - value of the objective funtion at the local rectangular element
    405:     g - gradient of the local function

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

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

    412:   double avgWq, sqGrad, avgWV, fl, fu;
    413:   double hx, hy, area, aread3, *wq, *wl;
    414:   double dvdx, dvdy;
    415:   int i;

    417:   hx = user->hx;
    418:   hy = user->hy;
    419:   area = user->area;
    420:   aread3 = area / 3.0;
    421:   wq = user->wq;
    422:   wl = user->wl;
    423:   i = coor[0];

    425:   /* lower triangle contribution */
    426:   dvdx = (x[0] - x[1]) / hx;
    427:   dvdy = (x[0] - x[2]) / hy;
    428:   sqGrad = dvdx * dvdx + dvdy * dvdy;
    429:   avgWq = (2.0 * wq[i] + wq[i+1]) / 6.0;
    430:   avgWV = (wl[i]*x[0] + wl[i+1]*x[1] + wl[i]*x[2]) / 3.0;
    431:   fl = avgWq * sqGrad - avgWV;

    433:   dvdx = dvdx * hy * avgWq;
    434:   dvdy = dvdy * hx * avgWq;
    435:   g[0] = ( dvdx + dvdy ) - wl[i] * aread3;
    436:   g[1] = ( -dvdx ) - wl[i+1] * aread3;
    437:   g[2] = ( -dvdy ) - wl[i] * aread3;

    439:   /* upper triangle contribution */
    440:   dvdx = (x[3] - x[2]) / hx;
    441:   dvdy = (x[3] - x[1]) / hy;
    442:   sqGrad = dvdx * dvdx + dvdy * dvdy;
    443:   avgWq = (2.0 * wq[i+1] + wq[i]) / 6.0;
    444:   avgWV = (wl[i+1]*x[1] + wl[i]*x[2] + wl[i+1]*x[3]) / 3.0;
    445:   fu = avgWq * sqGrad - avgWV;

    447:   dvdx = dvdx * hy * avgWq;
    448:   dvdy = dvdy * hx * avgWq;
    449:   g[1] += (-dvdy) - wl[i+1] * aread3;
    450:   g[2] +=  (-dvdx) - wl[i] * aread3;
    451:   g[3] = ( dvdx + dvdy ) - wl[i+1] * aread3;

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

    455:   return 0;
    456: } /* JBearLocalFunctionGradient */

    459: /*------- USER-DEFINED: routine to evaluate the Hessian
    460:            at a local (rectangular element) level       -------*/
    463: /*
    464:   JBearLocalHessian - Computes the Hessian of the local (partial) function
    465:          defined over the current rectangle

    467:   Input:
    468:     coor - vector with the indices of the position of current element
    469:              in the first, second and third directions
    470:     x - current local solution (over the rectangle only)
    471:     ptr - user-defined application context

    473:   Output:
    474:     H - Hessian matrix of the local function (wrt the four
    475:            points of the rectangle only)

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

    480:   AppCtx *user = (AppCtx*)ptr;
    481:   double wql, wqu;
    482:   double hx, hy, dydx, dxdy;
    483:   double *wq;
    484:   double wqldydx,wqldxdy,wqudydx,wqudxdy;
    485:   int i;

    487:   hx = user->hx;
    488:   hy = user->hy;
    489:   wq = user->wq;
    490:   i = coor[0];

    492:   dxdy = hx / hy;
    493:   dydx = hy / hx;
    494:   wql = (2.0 * wq[i] + wq[i+1]) / 6.0;
    495:   wqu = (wq[i] + 2.0 * wq[i+1]) / 6.0;
    496:   wqldydx = wql * dydx;
    497:   wqldxdy = wql * dxdy;
    498:   wqudydx = wqu * dydx;
    499:   wqudxdy = wqu * dxdy;

    501:           /* Hessian contribution at 0,0 */
    502:   H[0][0] = wqldxdy + wqldydx;
    503:   H[0][1] = H[1][0] = -wqldydx;
    504:   H[0][2] = H[2][0] = -wqldxdy;
    505:   H[0][3] = H[3][0] = 0.0;

    507:           /* Hessian contribution at 1,0 */
    508:   H[1][1] = wqldydx + wqudxdy;
    509:   H[1][2] = H[2][1] = 0.0;
    510:   H[1][3] = H[3][1] = -wqudxdy;

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

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

    519:   return 0;

    521: } /* JBearLocalHessian */

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

    529: /*
    530:   WholeJBearFunctionGradient - Evaluates function and gradient over the
    531:       whole grid

    533:   Input:
    534:     daapplication - TAO application object
    535:     da  - distributed array
    536:     X   - the current point, at which the function and gradient are evaluated
    537:     ptr - user-defined application context

    539:   Output:
    540:     f - value of the objective funtion at X
    541:     G - gradient at X
    542: */
    543: static int WholeJBearFunctionGradient(TAO_APPLICATION daapplication, DA da, Vec X, double *f, Vec G, void *ptr) {

    545:   AppCtx *user = (AppCtx*)ptr;
    546:   Vec localX, localG;
    547:   int info, i, j;
    548:   int xs, xm, gxs, gxm, xe, ys, ym, gys, gym, ye;
    549:   double **x, **g;
    550:   double floc = 0.0;
    551:   PetscScalar zero = 0.0;

    553:   double avgWq, sqGrad, avgWV, fl, fu;
    554:   double hx, hy, area, aread3, *wq, *wl;
    555:   double dvdx, dvdy;

    557:   hx = user->hx;
    558:   hy = user->hy;
    559:   area = user->area;
    560:   aread3= area/3.0;
    561:   wq = user->wq;
    562:   wl = user->wl;

    564:   info = DAGetLocalVector(da, &localX); CHKERRQ(info);
    565:   info = DAGetLocalVector(da, &localG); CHKERRQ(info);
    566:   info = VecSet(G, zero); CHKERRQ(info);
    567:   info = VecSet(localG, zero); CHKERRQ(info);

    569:   info = DAGlobalToLocalBegin(da, X, INSERT_VALUES, localX); CHKERRQ(info);
    570:   info = DAGlobalToLocalEnd(da, X, INSERT_VALUES, localX); CHKERRQ(info);

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

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

    578:   xe = gxs + gxm - 1;
    579:   ye = gys + gym - 1;
    580:   for (j = ys; j < ye; j++) {
    581:     for (i = xs; i < xe; i++) {

    583:       /* lower triangle contribution */
    584:       dvdx = (x[j][i] - x[j][i+1]) / hx;
    585:       dvdy = (x[j][i] - x[j+1][i]) / hy;
    586:       sqGrad = dvdx * dvdx + dvdy * dvdy;
    587:       avgWq = (2.0 * wq[i] + wq[i+1]) / 6.0;
    588:       avgWV = (wl[i]*x[j][i] + wl[i+1]*x[j][i+1] + wl[i]*x[j+1][i]) / 3.0;
    589:       fl = avgWq * sqGrad - avgWV;

    591:       dvdx = dvdx * hy * avgWq;
    592:       dvdy = dvdy * hx * avgWq;
    593:       g[j][i] +=  ( dvdx + dvdy ) - wl[i] * aread3;
    594:       g[j][i+1] += ( -dvdx ) - wl[i+1] * aread3;
    595:       g[j+1][i] += ( -dvdy ) - wl[i] * aread3;

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

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

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

    613:     }
    614:   }

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

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

    621:   info = DALocalToGlobalBegin(da, localG, G); CHKERRQ(info);
    622:   info = DALocalToGlobalEnd(da, localG, G); CHKERRQ(info);

    624:   info = DARestoreLocalVector(da, &localX); CHKERRQ(info);
    625:   info = DARestoreLocalVector(da, &localG); CHKERRQ(info);

    627:   info = PetscLogFlops((xe-xs) * (ye-ys) * 67 + 1); CHKERRQ(info);
    628:   return 0;
    629: } /* WholeJBearFunctionGradient */

    634: /*
    635:   WholeJBearHessian - Evaluates Hessian over the whole grid

    637:   Input:
    638:     daapplication - TAO application object
    639:     da  - distributed array
    640:     X   - the current point, at which the function and gradient are evaluated
    641:     ptr - user-defined application context

    643:   Output:
    644:     H - Hessian at X
    645: */
    646: static int WholeJBearHessian(TAO_APPLICATION daapplication, DA da, Vec X, Mat H, void *ptr) {

    648:   AppCtx *user = (AppCtx*)ptr;
    649:   int info, i, j, ind[4];
    650:   int xs, xm, gxs, gxm, xe, ys, ym, gys, gym, ye;
    651:   double smallH[4][4];
    652:   double wql, wqu;
    653:   double dydx, dxdy;
    654:   double *wq;
    655:   double wqldydx,wqldxdy,wqudydx,wqudxdy;
    656:   PetscTruth assembled;

    658:   wq = user->wq;
    659:   dydx = user->hy / user->hx;
    660:   dxdy = user->hx / user->hy;

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

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

    669:   xe = gxs + gxm - 1;
    670:   ye = gys + gym - 1;
    671:   for (j = ys; j < ye; j++) {
    672:     for (i = xs; i < xe; i++) {

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

    677:       wqldydx = wql * dydx;
    678:       wqldxdy = wql * dxdy;
    679:       wqudydx = wqu * dydx;
    680:       wqudxdy = wqu * dxdy;

    682:           /* Hessian contribution at 0,0 */
    683:       smallH[0][0] = wqldxdy + wqldydx;
    684:       smallH[0][1] = smallH[1][0] = -wqldydx;
    685:       smallH[0][2] = smallH[2][0] = -wqldxdy;
    686:       smallH[0][3] = smallH[3][0] = 0.0;

    688:           /* Hessian contribution at 1,0 */
    689:       smallH[1][1] = (wqldydx + wqudxdy);
    690:       smallH[1][2] = smallH[2][1] = 0.0;
    691:       smallH[1][3] = smallH[3][1] = -wqudxdy;

    693:           /* Hessian contribution at 0,1 */
    694:       smallH[2][2] = (wqldxdy + wqudydx);
    695:       smallH[2][3] = smallH[3][2] = -wqudydx;

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

    700:       ind[0] = (j-gys) * gxm + (i-gxs);
    701:       ind[1] = ind[0] + 1;
    702:       ind[2] = ind[0] + gxm;
    703:       ind[3] = ind[2] + 1;
    704:       info = MatSetValuesLocal(H,4,ind,4,ind,(PetscScalar*)smallH,ADD_VALUES); CHKERRQ(info);

    706:     }
    707:   }

    710:   info = MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY); CHKERRQ(info);
    711:   info = MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY); CHKERRQ(info);
    712:   info = MatSetOption(H, MAT_SYMMETRIC); CHKERRQ(info);

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

    718: } /* WholeJBearHessian */