#include "taoapp_petsc.h"  /*I  "tao.h"  I*/
#include "taoapp.h"
#include "src/petsctao/include/taopetsc.h"

extern int TaoVecDestroy(TaoVec*);
extern int TaoMatDestroy(TaoMat*);

#undef __FUNCT__  
#define __FUNCT__ "TaoPetscApplication::TaoPetscApplication"
TaoPetscApplication::TaoPetscApplication(MPI_Comm mpicomm){
   this->comm=mpicomm;
   this->Setup();
   return;
}

#undef __FUNCT__  
#define __FUNCT__ "TaoPetscApplication::Setup"
int TaoPetscApplication::Setup(){

  PetscFunctionBegin;

  this->taox=0;this->taoh=0;
  this->taofv=0; this->taoj=0;
  this->taofvll=0; this->taofvuu=0; this->taoAA=0;

  this->ksptmp=0;
  this->tao=0;
  this->papp=0;

  PetscFunctionReturn(0);
}


#undef __FUNCT__  
#define __FUNCT__ "TaoPetscApplication::TakeDown"
int TaoPetscApplication::TakeDown(){
  int info;

  PetscFunctionBegin;
  info = TaoVecDestroy(this->taox); CHKERRQ(info);
  info = TaoMatDestroy(this->taoh); CHKERRQ(info);
  info = TaoVecDestroy(this->taofv); CHKERRQ(info);
  info = TaoMatDestroy(this->taoj); CHKERRQ(info);
  info = TaoVecDestroy(this->taofvll); CHKERRQ(info);
  info = TaoVecDestroy(this->taofvuu); CHKERRQ(info);
  info = TaoMatDestroy(this->taoAA); CHKERRQ(info);
  if (this->ksptmp){ info=KSPDestroy(this->ksptmp); CHKERRQ(info);}
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "TaoPetscApplication::~TaoPetscApplication"
TaoPetscApplication::~TaoPetscApplication(){
  this->TakeDown();
  return;
}

#undef __FUNCT__  
#define __FUNCT__ "TaoPetscApplication::GetVariableVector"
int TaoPetscApplication::GetVariableVector(TaoVec **xx){
  int info;
  Vec V;
  PetscFunctionBegin;
  info=TaoAppGetSolutionVec(this->papp,&V); CHKERRQ(info);
  if (this->taox==0){
    info = TaoWrapPetscVec(V,&this->taox); CHKERRQ(info);
  }
  info = this->taox->SetVec(V); CHKERRQ(info);
  *xx= this->taox;
  PetscFunctionReturn(0);
}



#undef __FUNCT__  
#define __FUNCT__ "TaoPetscApplication::EvaluateObjectiveFunction"
int TaoPetscApplication::EvaluateObjectiveFunction(TaoVec *xx, double *ff){
  int     info;
  Vec X = (Vec) xx->VecObject;  
  PetscFunctionBegin;
  info=TaoAppComputeObjective(this->papp,X,ff);CHKERRQ(info);
  PetscFunctionReturn(0);
}


#undef __FUNCT__  
#define __FUNCT__ "TaoPetscApplication::EvaluateGradient"
int TaoPetscApplication::EvaluateGradient(TaoVec *xx, TaoVec *gg){
  int     info;
  Vec X = (Vec) xx->VecObject, G = (Vec) gg->VecObject;
  PetscFunctionBegin;
  info=TaoAppComputeGradient(this->papp,X,G);CHKERRQ(info);
  PetscFunctionReturn(0);
}



#undef __FUNCT__  
#define __FUNCT__ "TaoPetscApplication::EvaluateObjectiveAndGradient"
int TaoPetscApplication::EvaluateObjectiveAndGradient(TaoVec *xx, double *ff, TaoVec *gg){
  int     info;
  Vec X = (Vec) xx->VecObject, G = (Vec) gg->VecObject;
  PetscFunctionBegin;
  info=TaoAppComputeObjectiveAndGradient(this->papp,X,ff,G); CHKERRQ(info);
  PetscFunctionReturn(0);
}



#undef __FUNCT__  
#define __FUNCT__ "TaoPetscApplication::GetHessianMatrix"
int TaoPetscApplication::GetHessianMatrix(TaoMat **HH){
  int info;
  Mat H,HP;
  MatStructure flag;

  PetscFunctionBegin;
  info=TaoAppGetHessianMat(this->papp,&H,&HP); CHKERRQ(info);
  if (this->taoh==0){
    info = TaoWrapPetscMat(H,&this->taoh); CHKERRQ(info);
  }
  info=this->taoh->GetMatrix(0,0,&flag);CHKERRQ(info);
  info=this->taoh->SetMatrix(H,HP,flag);CHKERRQ(info);
  *HH=this->taoh;
  PetscFunctionReturn(0);
}


#undef __FUNCT__  
#define __FUNCT__ "TaoPetscApplication::EvaluateHessian"
int TaoPetscApplication::EvaluateHessian(TaoVec *xx, TaoMat *HH){
  int     info;
  Vec X = (Vec) xx->VecObject;
  TaoMatPetsc *MM=(TaoMatPetsc*)HH;
  Mat H, HPre;
  MatStructure pflag;

  PetscFunctionBegin;
  //  info = TaoAppGetGradientVec(this,&this->FDGrad);CHKERRQ(info);
  info = MM->GetMatrix(&H,&HPre,&pflag);CHKERRQ(info);
  info = TaoAppComputeHessian(this->papp, X, &H, &HPre, &pflag);
  info = MM->SetMatrix(H,HPre,pflag);CHKERRQ(info);
  PetscFunctionReturn(0);
}


#undef __FUNCT__  
#define __FUNCT__ "TaoPetscApplication::HessianSolve"
int TaoPetscApplication::HessianSolve(TaoVec *vv, TaoVec *ww, TaoTruth *success){
  int     info;
  Vec V = (Vec) vv->VecObject;
  Vec W = (Vec) ww->VecObject;
  PetscTruth success2;
  PetscFunctionBegin;
  info=TaoAppHessianSolve(this->papp,V,W,&success2);CHKERRQ(info);
  if (success2==PETSC_TRUE){*success=TAO_TRUE;} else { *success=TAO_FALSE;}
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "TaoPetscApplication::EvaluateVariableBounds"
int TaoPetscApplication::EvaluateVariableBounds(TaoVec *xxll, TaoVec *xxuu){
  int info;
  Vec XL=(Vec)xxll->VecObject,XU=(Vec)xxuu->VecObject;
  PetscFunctionBegin;
  info=TaoAppComputeVariableBounds(this->papp,XL,XU); CHKERRQ(info);
  PetscFunctionReturn(0);
}



#undef __FUNCT__  
#define __FUNCT__ "TaoPetscApplication::GetConstraintVector"
int TaoPetscApplication::GetConstraintVector(TaoVec **rr){
  int info;
  Vec R;
  PetscFunctionBegin;
  info=TaoAppGetFunctionVec(this->papp,&R); CHKERRQ(info);
  if (this->taofv==0){
    info = TaoWrapPetscVec(R,&this->taofv); CHKERRQ(info);
  }
  info = this->taofv->SetVec(R); CHKERRQ(info);
  *rr= this->taofv;
  PetscFunctionReturn(0);
}


#undef __FUNCT__  
#define __FUNCT__ "TaoPetscApplication::EvaluateConstraints"
int TaoPetscApplication::EvaluateConstraints(TaoVec *xx, TaoVec *RR){
  int     info;
  Vec X = (Vec) xx->VecObject, R = (Vec) RR->VecObject;

  PetscFunctionBegin;
  info=TaoAppComputeFunction(this->papp,X,R);CHKERRQ(info);
  PetscFunctionReturn(0);
}


#undef __FUNCT__  
#define __FUNCT__ "TaoPetscApplication::GetJacobian"
int TaoPetscApplication::GetJacobianMatrix(TaoMat **JJ){
  int info;
  Mat J,JP;
  MatStructure flag;

  PetscFunctionBegin;
  info=TaoAppGetJacobianMat(this->papp,&J,&JP); CHKERRQ(info);
  if (this->taoj==0){
    info = TaoWrapPetscMat(J,&this->taoj); CHKERRQ(info);
  }
  info=this->taoj->GetMatrix(0,0,&flag);CHKERRQ(info);
  info=this->taoj->SetMatrix(J,JP,flag);CHKERRQ(info);
  *JJ=this->taoj;
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "TaoPetscApplication::EvaluateJacobian"
int TaoPetscApplication::EvaluateJacobian(TaoVec *xx, TaoMat *JJ){

  int     info;
  Vec X = (Vec) xx->VecObject;
  Mat J,JPre;
  MatStructure pflag;
  TaoMatPetsc *JJJ=(TaoMatPetsc*)JJ;

  PetscFunctionBegin;
  info = JJJ->GetMatrix(&J,&JPre,&pflag);CHKERRQ(info);
  info = TaoAppComputeJacobian(this->papp,X,&J,&JPre,&pflag);CHKERRQ(info);
  info = JJJ->SetMatrix(J,JPre,pflag);CHKERRQ(info);
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "TaoPetscApplication::GetInequalityConstraints"
int TaoPetscApplication::GetInequalityConstraints(TaoVec **LL, TaoMat **JJ, TaoVec**UU){
  PetscFunctionBegin;
  *LL=this->taofvll;
  *JJ=this->taoAA;
  *UU=this->taofvuu;
  PetscFunctionReturn(0);
}

int MatCreateADA(Mat,Vec,Vec,Mat*);

#undef __FUNCT__
#define __FUNCT__ "TaoPetscApplication::CreateATDAMatrix"
int TaoPetscApplication::CreateATDAMatrix(TaoMat *A, TaoVec* D,TaoVec*X, TaoMat**AADAA){
  int info;
  Vec D1 = (Vec) D->VecObject;
  Vec D2 = (Vec) X->VecObject;
  Mat ADA;
  Mat pm=(Mat)A->MatObject;
  TaoMatPetsc *TT;
  PetscFunctionBegin;
  info = MatCreateADA(pm,D1,D2, &ADA); CHKERRQ(info);
  info = TaoWrapPetscMat(ADA, &TT); CHKERRQ(info);
  info = PetscObjectDereference((PetscObject)ADA);
  *AADAA=TT;
  PetscFunctionReturn(0);
}



#undef __FUNCT__  
#define __FUNCT__ "TaoPetscApplication::GetLinearSolver"
int TaoPetscApplication::GetLinearSolver(TaoMat *H, int stype, TaoLinearSolver **SS){
  MPI_Comm comm2;
  int info;
  MatType mtype;
  PC pc;
  KSP newksp;
  TaoLinearSolver *S;
  PetscTruth flg1,flg2;
  Mat pm=(Mat)H->MatObject;

  PetscFunctionBegin;
  //  info = ((TaoMatPetsc*)H)->CreateLinearSolver(flag,tksp); CHKERRQ(info);
  //  info = TaoGetKSP(this->papp,&newksp); CHKERRQ(info);
  newksp=0;
  if (this->ksptmp){
    info = TaoWrapKSP( newksp, (TaoLinearSolverPetsc**)&S ); CHKERRQ(info);
    PetscObjectDereference((PetscObject)(ksptmp)); CHKERRQ(info);
    *SS=S;
    this->ksptmp=0;
    PetscFunctionReturn(0);
  } 

  if (pm){
    info = PetscObjectGetComm((PetscObject)pm,&comm2); CHKERRQ(info);
    info = MatGetType(pm,&mtype); CHKERRQ(info);
    info = PetscStrncmp(mtype,MATSEQDENSE,10,&flg1); CHKERRQ(info);
    info = PetscStrncmp(mtype,MATMPIDENSE,10,&flg2); CHKERRQ(info);
  } else {
    comm2 = this->comm;
  }
  info = KSPCreate(comm,&newksp); CHKERRQ(info);
  info = KSPGetPC(newksp,&pc); CHKERRQ(info);
  if (flg1==PETSC_TRUE || flg2==PETSC_TRUE){
    info=PCSetType(pc,PCJACOBI); CHKERRQ(info);
  } else {
    info=PCSetType(pc,PCBJACOBI); CHKERRQ(info);
  }

  if (stype==300){
    info=KSPSetType(newksp,KSPCG); CHKERRQ(info);
    info=KSPSetFromOptions(newksp); CHKERRQ(info);
  } else if (stype==310){
    info=KSPSetFromOptions(newksp); CHKERRQ(info);
    info=KSPSetType(newksp,KSPGMRES); CHKERRQ(info);
  } else if (stype==220){
    info=KSPGetPC(newksp,&pc); CHKERRQ(info);
    info=PCSetType(pc,PCNONE); CHKERRQ(info);
    info=PCSetFromOptions(pc); CHKERRQ(info);
    info=KSPSetType(newksp,KSPQCG); CHKERRQ(info);
  } else if (stype==110){
    info=KSPSetFromOptions(newksp); CHKERRQ(info);
    info=KSPSetType(newksp,KSPGMRES); CHKERRQ(info);
  } else {
    info=KSPSetType(newksp,KSPGMRES); CHKERRQ(info);
    info=KSPSetFromOptions(newksp); CHKERRQ(info);
  }

  info = TaoWrapKSP( newksp, (TaoLinearSolverPetsc**)&S ); CHKERRQ(info);
  //  info=KSPDestroy(newksp);CHKERRQ(info); newksp=0;
  PetscObjectDereference((PetscObject)(newksp)); CHKERRQ(info);
  *SS=S;
  //  info = TaoAppSetKSP(this->papp,newksp); CHKERRQ(info);

  PetscFunctionReturn(0);
}



#undef __FUNCT__  
#define __FUNCT__ "TaoPetscApplication::InitializeVariables"
int TaoPetscApplication::InitializeVariables(TaoVec *xx){
  PetscFunctionBegin; 
  PetscFunctionReturn(0);
}


#undef __FUNCT__  
#define __FUNCT__ "TaoPetscApplication::Monitor"
int TaoPetscApplication::Monitor(){
  int info;
  PetscFunctionBegin;
  info=TaoAppMonitor(this->papp); CHKERRQ(info);
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "TaoPetscApplication::Monitor2"
int TaoPetscApplication::Monitor2(TaoVec *xx, TaoVec *gl, TaoVec *dx, TaoTruth *flag){
  int     info;
  PetscTruth flag2;
  Vec GL = (Vec) gl->VecObject;
  PetscFunctionBegin;
  *flag=TAO_FALSE;
  info=TaoAppCheckConvergence(this->papp, GL,&flag2); CHKERRQ(info);
  if (flag2==PETSC_TRUE) *flag=TAO_TRUE;
  PetscFunctionReturn(0);
}


