#include "taoapp.h"         /*I  "tao.h"  I*/
#include "petscsnes.h"

extern int TaoAppSetColoring(TAO_APPLICATION, ISColoring);
extern int TaoAppGetColoring(TAO_APPLICATION, ISColoring*);


#undef __FUNCT__  
#define __FUNCT__ "TAOFiniteDifferenceGradient"
static int Ftemp(SNES snes ,Vec X,Vec G,void*ctx){
  int info;
  TAO_APPLICATION taoapp = (TAO_APPLICATION)ctx;
  PetscFunctionBegin;
  info=TaoAppComputeGradient(taoapp,X,G); CHKERRQ(info);
  PetscFunctionReturn(0);
}

EXTERN_C_BEGIN


#undef __FUNCT__  
#define __FUNCT__ "TaoAppDefaultComputeHessian"
/*@C
   TaoAppDefaultComputeHessian - Computes the Hessian using finite differences. 

   Collective on TAO_APPLICATION

   Input Parameters:
+  taoapp - the TAO_APPLICATION context 
.  V - compute Hessian at this point
-  ctx - the TAO_APPLICATION structure, cast to (void*)

   Output Parameters:
+  H - Hessian matrix (not altered in this routine)
.  B - newly computed Hessian matrix to use with preconditioner (generally the same as H)
-  flag - flag indicating whether the matrix sparsity structure has changed

   Options Database Key:
+  -tao_fd - Activates TaoAppDefaultComputeHessian()
-  -tao_view_hessian - view the hessian after each evaluation using PETSC_VIEWER_STDOUT_WORLD

   Level: intermediate

   Notes:
   This routine is slow and expensive, and is not currently optimized
   to take advantage of sparsity in the problem.  Although
   TaoAppDefaultComputeHessian() is not recommended for general use
   in large-scale applications, It can be useful in checking the
   correctness of a user-provided Hessian.

   Note:
   The gradient evaluation must be set using the routine TaoAppSetGradientRoutine().

.keywords: TAO_APPLICATION, finite differences, Hessian

.seealso: TaoAppSetHessianRoutine(), TaoAppDefaultComputeHessianColor(), SNESDefaultComputeJacobian(),
          TaoAppSetGradientRoutine()

@*/
int TaoAppDefaultComputeHessian(TAO_APPLICATION taoapp,Vec V,Mat *H,Mat *B,
			     MatStructure *flag,void *ctx){
  int                  info;
  MPI_Comm             comm;
  Vec                  G;
  SNES                 snes;

  PetscFunctionBegin;

  PetscValidHeaderSpecific(V,VEC_COOKIE,2);
  info = VecDuplicate(V,&G);CHKERRQ(info);

  info = PetscLogInfo((G,"TAO Using finite differences w/o coloring to compute matrix\n")); CHKERRQ(info);

  info = TaoAppComputeGradient(taoapp,V,G); CHKERRQ(info);

  info = PetscObjectGetComm((PetscObject)(*H),&comm);CHKERRQ(info);
  info = SNESCreate(comm,&snes);CHKERRQ(info);

  info = SNESSetFunction(snes,G,Ftemp,ctx);CHKERRQ(info);
  info = SNESSetJacobian(snes,*H,*B,SNESDefaultComputeJacobian,ctx);CHKERRQ(info);
  //  info = SNESComputeJacobian(snes,V,H,B,flag);CHKERRQ(info);
  info = SNESDefaultComputeJacobian(snes,V,H,B,flag,ctx);CHKERRQ(info);

  info = SNESDestroy(snes);CHKERRQ(info);
  
  info = VecDestroy(G);CHKERRQ(info);
  
  PetscFunctionReturn(0);
}

EXTERN_C_END



EXTERN_C_BEGIN
#undef __FUNCT__  
#define __FUNCT__ "TaoAppDefaultComputeHessianColor"
/*@C
   TaoAppDefaultComputeHessianColor - Computes the Hessian using colored finite differences. 

   Collective on TAO_APPLICATION

   Input Parameters:
+  tao - the TAO_APPLICATION context
.  V - compute Hessian at this point
-  ctx - the TAO_APPLICATION structure, cast to (void*)

   Output Parameters:
+  H - Hessian matrix (not altered in this routine)
.  B - newly computed Hessian matrix to use with preconditioner (generally the same as H)
-  flag - flag indicating whether the matrix sparsity structure has changed

   Options Database Keys:
+  -mat_fd_coloring_freq <freq>
-  -tao_view_hessian - view the hessian after each evaluation using PETSC_VIEWER_STDOUT_WORLD

   Level: intermediate

   Note:
   The gradient evaluation must be set using the routine TaoSetPetscGradient().

 .keywords: TAO_APPLICATION, finite differences, Hessian, coloring, sparse

.seealso: TaoAppSetHessianRoutine(), TaoAppDefaultComputeHessian(),SNESDefaultComputeJacobianColor(), 
          TaoAppSetGradientRoutine(), TaoAppSetColoring()

@*/
int TaoAppDefaultComputeHessianColor(TAO_APPLICATION taoapp, Vec V, Mat *HH,Mat *BB,
				  MatStructure *flag,void *ctx){
  int                 info;
  MPI_Comm            comm;
  Vec                 G=0;
  Mat                 H=*HH,B=*BB;
  SNES                snes;
  ISColoring          iscoloring;
  MatFDColoring       matfdcoloring;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(H,MAT_COOKIE,3);  
  PetscValidHeaderSpecific(B,MAT_COOKIE,4);  
  PetscCheckSameComm(V,2,H,3);
  PetscCheckSameComm(H,3,B,4);

  info = TaoAppGetColoring(taoapp,&iscoloring); CHKERRQ(info);
  if (!iscoloring){
    SETERRQ(1,"Must set coloring before using this routine.  Try Finite Differences without coloring\n");
  }
  info = VecDuplicate(V,&G);CHKERRQ(info);

  info=PetscLogInfo((G,"TAO computing matrix using finite differences and coloring\n")); CHKERRQ(info);

  info=TaoAppComputeGradient(taoapp,V,G); CHKERRQ(info);

  info = PetscObjectGetComm((PetscObject)(H),&comm);CHKERRQ(info);
  info = SNESCreate(comm,&snes);CHKERRQ(info);

  info = MatFDColoringCreate(H,iscoloring,&matfdcoloring);CHKERRQ(info);
  info = MatFDColoringSetFunction(matfdcoloring,(int (*)(void)) Ftemp,ctx);CHKERRQ(info);
  info = MatFDColoringSetFromOptions(matfdcoloring);CHKERRQ(info);

  info = SNESSetFunction(snes,G,Ftemp,ctx);CHKERRQ(info);
  info = SNESSetJacobian(snes,H,B,SNESDefaultComputeJacobianColor,(void*)matfdcoloring);CHKERRQ(info);
  info = SNESDefaultComputeJacobianColor(snes,V,HH,BB,flag,matfdcoloring);CHKERRQ(info);

  info = MatFDColoringDestroy(matfdcoloring);CHKERRQ(info);
  info = SNESDestroy(snes);CHKERRQ(info);
  
  info = VecDestroy(G);CHKERRQ(info);
  PetscFunctionReturn(0);
}
EXTERN_C_END

#undef __FUNCT__  
#define __FUNCT__ "TaoAppSetFiniteDifferencesOptions"
/* @
  TaoAppSetFiniteDifferencesOptions - Sets various TAO parameters from user options

   Collective on TAO_APPLICATION

   Input Parameters:
+  taoapp - the TAO Application (optional)

   Level: beginner

.keywords:  options, finite differences

.seealso: TaoSolveApplication();

@ */
int TaoAppSetFiniteDifferencesOptions(TAO_APPLICATION taoapp){
  int info;
  PetscTruth flg;

  PetscFunctionBegin;

  PetscValidHeaderSpecific(taoapp,TAO_APP_COOKIE,1);

  flg=PETSC_FALSE;
  info = PetscOptionsName("-tao_fd","use finite differences for Hessian","TaoAppDefaultComputeHessian",&flg);CHKERRQ(info);
  if (flg) {
    info = TaoAppSetHessianRoutine(taoapp,TaoAppDefaultComputeHessian,(void*)taoapp);CHKERRQ(info);
    info = PetscLogInfo((taoapp,"TaoSetOptions: Setting default finite difference Hessian matrix\n")); CHKERRQ(info);
  }

  flg=PETSC_FALSE;
  info = PetscOptionsName("-tao_fd_coloring","use finite differences with coloring to compute Hessian","TaoAppDefaultComputeHessianColor",&flg);CHKERRQ(info);
  if (flg) {
    info = TaoAppSetHessianRoutine(taoapp,TaoAppDefaultComputeHessianColor,(void*)taoapp);CHKERRQ(info);
    info = PetscLogInfo((taoapp,"TaoAppSetFromOptions: Use finite differencing with coloring to compute Hessian \n")); CHKERRQ(info);
  }
    
  PetscFunctionReturn(0);
}


static char TaoAppColoringXXX[] = "TaoColoring";

typedef struct {
  ISColoring coloring;
} TaoAppColorit;

#undef __FUNCT__  
#define __FUNCT__ "TaoAppDestroyColoringXXX"
static int TaoAppDestroyColoringXXX(void*ctx){
  //  int info;
  TaoAppColorit *mctx=(TaoAppColorit*)ctx;
  PetscFunctionBegin;
  if (mctx){
    /*
    if (mctx->coloring){  
      info = ISColoringDestroy(mctx->coloring);CHKERRQ(info);
    }
    */
    PetscFree(mctx);
  }
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "TaoAppSetColoring"
/*@
   TaoAppSetColoring - Set the matrix coloring to be used when computing the
   Hessian by finite differences.

   Collective on TAO_APPLICATION

   Input Parameters:
+  app - the TAO_APPLICATION context
-  coloring - the coloring to be used

   Level: intermediate

   Note:
   The gradient evaluation must be set using the routine TaoSetPetscGradient().

 .keywords: TAO_APPLICATION, finite differences, Hessian, coloring, sparse

.seealso: TaoAppSetHessianRoutine(), TaoAppDefaultComputeHessianColor(),
          TaoAppSetGradientRoutine()

@*/
int TaoAppSetColoring(TAO_APPLICATION taoapp, ISColoring coloring){
  int ii,info;
  TaoAppColorit *mctx;
  PetscFunctionBegin;
  PetscValidHeaderSpecific(taoapp,TAO_APP_COOKIE,1);
  
  info = TaoAppQueryForObject(taoapp,TaoAppColoringXXX,(void**)&mctx); CHKERRQ(info);
  if (mctx==0){
    info=PetscNew(TaoAppColorit,(void**)&mctx); CHKERRQ(info);
    info=TaoAppAddObject(taoapp,TaoAppColoringXXX,(void*)mctx,&ii); CHKERRQ(info);
    info=TaoAppSetDestroyRoutine(taoapp,TaoAppDestroyColoringXXX, (void*)mctx); CHKERRQ(info);
  }
  /*
  if (coloring){
  info=PetscObjectReference((PetscObject)coloring); CHKERRQ(info);
  }
  if (mctx->coloring){
       info=ISColoringDestroy(mctx->coloring); CHKERRQ(info);
  }
  */
  mctx->coloring=coloring;  
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "TaoAppGetColoring"
/*@C
   TaoAppGetColoring - Set the matrix coloring to be used when computing the
   Hessian by finite differences.

   Collective on TAO_APPLICATION

   Input Parameters:
+  app - the TAO_APPLICATION context
-  coloring - the coloring to be used

   Level: advanced

   Note:
   The gradient evaluation must be set using the routine TaoAppSetGradientRoutine().

 .keywords: TAO_APPLICATION, finite differences, Hessian, coloring, sparse

.seealso: TaoAppSetHessianRoutine(), TaoAppDefaultComputeHessianColor(),
          TaoAppSetGradientRoutine()

@*/
int TaoAppGetColoring(TAO_APPLICATION taoapp, ISColoring *coloring){
  int info;
  TaoAppColorit *mctx;
  PetscFunctionBegin;
  PetscValidHeaderSpecific(taoapp,TAO_APP_COOKIE,1);
  if (coloring){
    info = TaoAppQueryForObject(taoapp,TaoAppColoringXXX,(void**)&mctx); CHKERRQ(info);
    if (mctx){
      *coloring=mctx->coloring;
    } else {
      *coloring=0;
    }
  }
  PetscFunctionReturn(0);
}


#undef __FUNCT__  
#define __FUNCT__ "TaoAppAddFiniteDifferences"
int TaoAppAddFiniteDifferences(TAO_APPLICATION taoapp){
  int info;
  PetscFunctionBegin;
  info = TaoAppSetOptionsRoutine(taoapp,TaoAppSetFiniteDifferencesOptions); CHKERRQ(info);
  PetscFunctionReturn(0);
}

int MatTAOMFSetGradient(Mat mat,Vec v,int (*func)(TAO_APPLICATION,Vec,Vec,void *),void *funcctx){
  PetscFunctionBegin;
  PetscFunctionReturn(0);
}

