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

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

/*
  Include "tao.h" so we can use TAO solvers.
  petscda.h for distributed array
  ad_deriv.h for AD gradient
*/

#include "petscda.h"
#include "tao.h"
#include "taodaapplication.h"

static char help[] = "Steady-State Combustion.\n\
We solve the Steady-State Combustion problem (MINPACK-2 test suite) in a 2D\n\
rectangular domain, using distributed arrays (DAs) to partition the parallel grid.\n\
The command line options include:\n\
  -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
  -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
  -nlevels <nlevels>, where <nlevels> = number of levels in multigrid\n\
  -byelement, if computation is made by functions on rectangular elements\n\
  -adic, if AD is used (AD is not used by default)\n\
  -par <parameter>, where <parameter> indicates the problem's nonlinearity\n\
     parameter lambda (0 <= par <= 6.81)\n\n";

/*T
   Concepts: TAO - Solving a bounded minimization problem
   Routines: TaoInitialize(); TaoFinalize();
   Routines: TaoCreate(); TaoDestroy();
   Routines: DAApplicationCreate(); DAApplicationDestroy();
   Routines: DAAppSetVariableBoundsRoutine();
   Routines: DAAppSetElementObjectiveAndGradientRoutine();
   Routines: DAAppSetElementHessianRoutine();
   Routines: DAAppSetObjectiveAndGradientRoutine();
   Routines: DAAppSetADElementFunctionGradient();
   Routines: DAAppSetHessianRoutine();
   Routines: TaoAppSetOptions();
   Routines: TaoGetSolutionStatus(); TaoDAAppSolve();
   Routines: DAAppSetMonitor(); TaoView();
   Routines: DAAppGetSolution();
   Routines: DAAppGetInterpolationMatrix();
   Processors: n
T*/

/*
   User-defined application context - contains data needed by the
   application-provided call-back routines.
*/

typedef struct {
  InactiveDouble      param;
  InactiveDouble      hx, hy;        /* increment size in both directions */
  InactiveDouble      area;          /* area of the triangles */
} ADFGCtx;

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

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

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

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

static int DAFixBoundary(TAO_APPLICATION, DA, Vec, Vec, void*);

static int MyGridMonitorBefore(TAO_APPLICATION, DA, int, void *);

#undef __FUNCT__
#define __FUNCT__ "main"

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

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

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

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

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

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

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

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

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

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

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

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

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

  } else {

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

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

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

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

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

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

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

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

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

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

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

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

  PreLoadEnd();

  /* Finalize TAO */
  TaoFinalize();
  PetscFinalize();

  return 0;
} /* main */



/*----- The following two routines
   
  MyGridMonitorBefore    MyGridMonitorAfter

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

#undef __FUNCT__
#define __FUNCT__ "MyGridMonitorBefore"
static int MyGridMonitorBefore(TAO_APPLICATION myapp, DA da, int level, void *ctx) {

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

  info = DAGetInfo(da,PETSC_NULL,&mx,&my,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
		   PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(info);
  user->mx = mx;
  user->my = my;
  user->hx = 1.0 / (user->mx - 1);
  user->hy = 1.0 / (user->my - 1);
  user->area = 0.5 * user->hx * user->hy;
  user->fgctx.hx   = user->hx;
  user->fgctx.hy   = user->hy;
  user->fgctx.area = user->area;
  user->fgctx.param = user->param;

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

  return 0;
}

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

#undef __FUNCT__
#define __FUNCT__ "AppCtxInitialize"
/*
  AppCtxInitialize - Sets initial values for the application context parameters

  Input:
    ptr - void user-defined application context

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

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

  /* Specify dimension of the problem */
  user->param = 5.0;
  user->mx = 11;
  user->my = 11;

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

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

  return 0;
} /* AppCtxInitialize */



#undef __FUNCT__
#define __FUNCT__ "FormInitialGuess"
static int FormInitialGuess(DA da, Vec X, AppCtx *ctx)
{
  int    info, i, j, mx, my;
  int    xs, ys, xm, ym, xe, ye;
  PetscReal hx, hy, temp, val, lambda;
  double **x;

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

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

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

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

  return 0;
}


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

#undef __FUNCT__
#define __FUNCT__ "DAFixBoundary"
/*
  FormBounds - Forms bounds on the variables

  Input:
    user - user-defined application context

  Output:
    XL - vector of lower bounds
    XU - vector of upper bounds

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

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

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

  /* Compute initial guess over locally owned part of the grid */
  for (j = ys; j < ys+ym; j++){
    for (i = xs; i < xs+xm; i++){
      if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) {
        xl[j][i] = xu[j][i] = 0.0;
      } else {
        xl[j][i] = lb;
        xu[j][i] = ub;
      }
    }
  }

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

  return 0;
} /* DAFixBoundary */


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

#undef __FUNCT__
#define __FUNCT__ "CombLocalFunctionGradient"
/*
  CombLocalFunctionGradient - Evaluates function and gradient over the 
      local rectangular element

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

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

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

  AppCtx *user = (AppCtx*)ptr;

  double lambdad3, hx, hy, area;
  double fquad, fexp, dvdx, dvdy;

  lambdad3 = user->param / 3.0;
  hx = user->hx;
  hy = user->hy;
  area = user->area;

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

  dvdx = 0.5 * dvdx * hy;
  dvdy = 0.5 * dvdy * hx;
  g[0] = dvdx + dvdy - exp(x[0]) * area * lambdad3;
  g[1] = -dvdx - 2.0 * exp(x[1]) * area * lambdad3;
  g[2] = -dvdy - 2.0 * exp(x[2]) * area * lambdad3;

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

  dvdx = 0.5 * dvdx * hy;
  dvdy = 0.5 * dvdy * hx;
  g[1] += -dvdy;
  g[2] += -dvdx;
  g[3] = dvdx + dvdy - exp(x[3]) * area * lambdad3;


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

  return 0;
} /* CombLocalFunctionGradient */


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

#undef __FUNCT__
#define __FUNCT__ "CombLocalHessian"
/*
  CombLocalHessian - Computes the Hessian of the local (partial) function
         defined over the current rectangle

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

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

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

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

  hx = user->hx;
  hy = user->hy;
  lambdad3 = user->param / 3.0;
  area = user->area;
  dxdy = hx/hy;
  dydx = hy/hx;
  diagxy = 0.5 * (dxdy + dydx);
  bandxy = -0.5 * dxdy;
  bandyx = -0.5 * dydx;

          /* Hessian contribution at 0,0 */
  dexp = exp(x[0]) * area * lambdad3;
  H[0][0] = diagxy - dexp;
  H[0][1] = H[1][0] = bandyx;
  H[0][2] = H[2][0] = bandxy;
  H[0][3] = H[3][0] = 0.0;

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

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

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

  return 0;
} /* CombLocalHessian */


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

#undef __FUNCT__
#define __FUNCT__ "WholeCombFunctionGradient"
/*
  WholeCombFunctionGradient - Evaluates function and gradient over the 
      whole grid

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

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

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

  double lambdad3, hx, hy, area;
  double fquad, fexp, dvdx, dvdy;

  lambdad3 = user->param / 3.0;
  hx = user->hx;
  hy = user->hy;
  area = user->area;

  info = DAGetLocalVector(da, &localX); CHKERRQ(info);
  info = DAGetLocalVector(da, &localG); CHKERRQ(info);
  info = VecSet(G, zero); CHKERRQ(info);
  info = VecSet(localG, zero); CHKERRQ(info);

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

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

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

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

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

      dvdx = 0.5 * dvdx * hy;
      dvdy = 0.5 * dvdy * hx;
      g[j][i] += dvdx + dvdy - exp(x[j][i]) * area * lambdad3;
      g[j][i+1] += -dvdx - 2.0 * exp(x[j][i+1]) * area * lambdad3;
      g[j+1][i] += -dvdy - 2.0 * exp(x[j+1][i]) * area * lambdad3;

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

      dvdx = 0.5 * dvdx * hy;
      dvdy = 0.5 * dvdy * hx;
      g[j][i+1] += -dvdy;
      g[j+1][i] += -dvdx;
      g[j+1][i+1] += dvdx + dvdy - exp(x[j+1][i+1]) * area * lambdad3;

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

    }
  }

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

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

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

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

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

} /* WholeCombFunctionGradient */


/*------- USER-DEFINED: routine to evaluate the Hessian 
          at the whole grid             -------*/
#undef __FUNCT__
#define __FUNCT__ "WholeCombHessian"
/*
  WholeCombHessian - Evaluates Hessian over the whole grid

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

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

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

  double hx, hy, lambdad3, area, dxdy, dydx;
  double diagxy, dexp, bandxy, bandyx;
  PetscTruth assembled;


  hx = user->hx;
  hy = user->hy;
  lambdad3 = user->param / 3.0;
  area = user->area;
  dxdy = hx/hy;
  dydx = hy/hx;
  diagxy = 0.5 * (dxdy + dydx);
  bandxy = -0.5 * dxdy;
  bandyx = -0.5 * dydx;

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

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

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

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

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

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

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

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

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

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

    }
  }

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

  info = MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY); CHKERRQ(info);
  info = MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY); CHKERRQ(info);
  info = MatSetOption(H, MAT_SYMMETRIC); CHKERRQ(info);

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

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

} /* WholeCombHessian */
