#include "lmvmmat.h"
#include "tao_general.h"
#include "taovec.h"
#include "stdio.h"

#undef __FUNCT__
#define __FUNCT__ "TaoLMVMMat::TaoLMVMMat"
TaoLMVMMat::TaoLMVMMat( TaoVec *tv, int nlm){

  TaoFunctionBegin;

  lm = nlm;
  lmnow = 0;
  eps = 2.2e-16;
  eps = 2.2e-11;
  iter = 0;
  rejects = 0;

  tv->CloneVecs(nlm+1,&S);
  tv->CloneVecs(nlm+1,&Y);

  rho = new double[nlm+1];
  beta = new double[nlm+1];

  Gprev = 0;
  Xprev = 0;
  H0 = 0;
  H0Vec = 0;

  H0default=TAO_TRUE;

  return;
}

#undef __FUNCT__
#define __FUNCT__ "TaoLMVMMat::TaoLMVMMat"
TaoLMVMMat::~TaoLMVMMat(){
  int info;

  info = S[0]->DestroyVecs(lm+1,S);
  info = Y[0]->DestroyVecs(lm+1,Y);
  /*  if (H0) delete H0; */
  if (H0Vec) delete H0Vec;

  delete[] rho;
  delete[] beta;
  if (info) rho=0;
}


#undef __FUNCT__
#define __FUNCT__ "TaoLMVMMat::InitialApproximation()"
int TaoLMVMMat::InitialApproximation(TaoVec *HH0){
  TaoFunctionBegin;
  TaoFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "TaoLMVMMat::SetH0()"
int TaoLMVMMat::SetH0(TaoMat* HH0){
  int info;
  TaoFunctionBegin;
  if (this->H0){ 
    /*    delete this->H0; Application deletes it */
    this->H0=0;
  };
  if (this->H0Vec){ delete this->H0Vec;this->H0Vec=0;};
  this->H0=HH0;
  this->H0default=TAO_TRUE;
  if (HH0){
    info=this->S[0]->Clone(&this->H0Vec);CHKERRQ(info);
    this->H0default=TAO_FALSE;
  }
  TaoFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "TaoLMVMMat::reset"
int TaoLMVMMat::reset(){
  int i;
  TaoFunctionBegin;
  Gprev = Y[lm];
  Xprev = S[lm];
  for (i=0;i<lm;i++){
    rho[i]=0.0;
  }
  rho[0]=1.0;
  y0normsquared = 1.0;
  iter=0;
  rejects=0;
  lmnow = 0;
  TaoFunctionReturn(0);
}


#undef __FUNCT__
#define __FUNCT__ "TaoLMVMMat::update"
int TaoLMVMMat::update(TaoVec* x, TaoVec* g){
  int i,info;
  double rhotemp,rhotol;
  double y0temp;
  TaoFunctionBegin;

  if (this->iter==0){

    info=this->reset(); CHKERRQ(info);

  } else {
    
    info = Gprev->Aypx(-1.0,g); CHKERRQ(info);
    info = Xprev->Aypx(-1.0,x); CHKERRQ(info);
    info = Gprev->Dot(Xprev,&rhotemp); CHKERRQ(info);
    info = Gprev->Norm2squared(&y0temp); CHKERRQ(info);
    rhotol=eps*y0temp;
    if (rhotemp > rhotol){
      lmnow = TaoMin(lmnow+1,lm);
      for (i=lm-1;i>=0;i--){
	S[i+1]=S[i];
	Y[i+1]=Y[i];
	rho[i+1]=rho[i];
      }
      S[0]=Xprev;
      Y[0]=Gprev;
      rho[0]=1.0/rhotemp;
      y0normsquared=y0temp;
      Xprev=S[lm]; Gprev=Y[lm];

    } else { 
      info = PetscLogInfo((0,"LMVM update rejected: iter= %d, %5.2e<=%5.2e \n",iter,rhotemp,rhotol)); CHKERRQ(info);
      rejects++;
    }
      
  }
  
  iter++;
  info = Xprev->CopyFrom(x); CHKERRQ(info);
  info = Gprev->CopyFrom(g); CHKERRQ(info);

  TaoFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "TaoLMVMMat::Solve"
int TaoLMVMMat::Solve(TaoVec* tb, TaoVec* dx, TaoTruth *tt){

  int      ll,info;
  double   sq, yq,dd;

  TaoFunctionBegin;
  if (lmnow<1){
    rho[0]=1.0; y0normsquared = 1.0;
  }

  info = dx->CopyFrom(tb); CHKERRQ(info);
  
  for (ll = 0; ll<lmnow; ll++){
    
    info = dx->Dot(S[ll],&sq); CHKERRQ(info);
    beta[ll] = sq*rho[ll];
    info = dx->Axpy(-beta[ll], Y[ll]); CHKERRQ(info);

  }
  
  if (this->H0default==TAO_TRUE){
    info = dx->Scale( 1.0 / (rho[0]*y0normsquared) ); CHKERRQ(info);
  } else {
    info=this->H0Vec->CopyFrom(dx);CHKERRQ(info);
    info=dx->Scale( 1.0 / (rho[0]*y0normsquared) ); CHKERRQ(info);
    info=this->H0->Solve(H0Vec,dx,tt);CHKERRQ(info);
    info=dx->Dot(H0Vec,&dd);CHKERRQ(info);
    if (dd<=0){
      info=this->H0Vec->CopyFrom(dx);CHKERRQ(info);
      info=dx->Scale( 1.0 / (rho[0]*y0normsquared) ); CHKERRQ(info);
      info = PetscLogInfo((0,"REJECTING Application Hessian Solve: Inner product= %4.2e \n",dd)); CHKERRQ(info);
    } else {
      info = PetscLogInfo((0,"Using Application Hessian Solve: Inner product= %4.2e \n",dd)); CHKERRQ(info);
    }
  }

  for (ll=lmnow-1; ll>=0; ll--){
    
    info = dx->Dot(Y[ll],&yq); CHKERRQ(info);
    info = dx->Axpy( beta[ll] - yq * rho[ll] , S[ll]); CHKERRQ(info);
    
  }

  *tt = TAO_TRUE;
  TaoFunctionReturn(0);
}


#undef __FUNCT__
#define __FUNCT__ "TaoLMVMMat::GetX0"
int TaoLMVMMat::GetX0(TaoVec* x){
  int i,info;

  TaoFunctionBegin;
  info=x->CopyFrom(Xprev);CHKERRQ(info);
  for (i=0;i<this->lmnow;i++){
    info = x->Axpy(-1.0, S[i]); CHKERRQ(info);
  }
  TaoFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "TaoLMVMMat::Clone"
int TaoLMVMMat::Clone(TaoMat** tm){
  TaoFunctionBegin;
  SETERRQ(56,"Operation not defined");
  /* TaoFunctionReturn(1); */
}

#undef __FUNCT__
#define __FUNCT__ "TaoLMVMMat::Presolve"
int TaoLMVMMat::Presolve(){
  TaoFunctionBegin;
  TaoFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "TaoLMVMMat::CopyFrom"
int TaoLMVMMat::CopyFrom(TaoMat* tm){
  TaoFunctionBegin;
  SETERRQ(56,"Operation not defined");
  /* TaoFunctionReturn(1); */
}

#undef __FUNCT__
#define __FUNCT__ "TaoLMVMMat::GetDimension"
int TaoLMVMMat::GetDimensions( int* m, int* n ){
  TaoFunctionBegin;
  SETERRQ(56,"Operation not defined");
  /* TaoFunctionReturn(1); */
}

#undef __FUNCT__
#define __FUNCT__ "TaoLMVMMat::Multiply"
int TaoLMVMMat::Multiply(TaoVec* tx,TaoVec* ty){
  TaoFunctionBegin;
  SETERRQ(56,"Operation not defined");
  /* TaoFunctionReturn(1); */
}

#undef __FUNCT__
#define __FUNCT__ "TaoLMVMMat::MultiplyAdd"
int TaoLMVMMat::MultiplyAdd(TaoVec* tx,TaoVec* ty,TaoVec* tw){
  TaoFunctionBegin;
  SETERRQ(56,"Operation not defined");
  /* TaoFunctionReturn(1); */
}

#undef __FUNCT__
#define __FUNCT__ "TaoLMVMMat::MultiplyTranspose"
int TaoLMVMMat::MultiplyTranspose(TaoVec* tx,TaoVec* ty){
  TaoFunctionBegin;
  SETERRQ(56,"Operation not defined");
  /* TaoFunctionReturn(1); */
}

#undef __FUNCT__
#define __FUNCT__ "TaoLMVMMat::MultiplyTransposeAdd"
int TaoLMVMMat::MultiplyTransposeAdd(TaoVec* tx,TaoVec* ty,TaoVec* yw){
  TaoFunctionBegin;
  SETERRQ(56,"Operation not defined");
  /* TaoFunctionReturn(1); */
}


#undef __FUNCT__
#define __FUNCT__ "TaoLMVMMat::SetDiagonal"
int TaoLMVMMat::SetDiagonal(TaoVec* tv){
  TaoFunctionBegin;
  SETERRQ(56,"Operation not defined");
  /* TaoFunctionReturn(1); */
}

#undef __FUNCT__
#define __FUNCT__ "TaoLMVMMat::AddDiagonal"
int TaoLMVMMat::AddDiagonal(TaoVec* tv){
  TaoFunctionBegin;
  SETERRQ(56,"Operation not defined");
  /* TaoFunctionReturn(1); */
}

#undef __FUNCT__
#define __FUNCT__ "TaoLMVMMat::View"
int TaoLMVMMat::View(){
  int i,info;
  TaoFunctionBegin;
  printf("Xprev \n");
  info = Xprev->View(); CHKERRQ(info);
  printf("Gprev \n");
  info = Gprev->View(); CHKERRQ(info);
  for (i=0;i<lmnow; i++){
    printf("rho: %4.2e \n",this->rho[i]);
    printf("S %d \n",i);
    info = this->S[i]->View(); CHKERRQ(info);
    printf("Y %d \n",i);
    info = this->Y[i]->View(); CHKERRQ(info);
  }
  TaoFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "TaoLMVMMat::GetDiagonal"
int TaoLMVMMat::GetDiagonal(TaoVec* tv){
  TaoFunctionBegin;
  SETERRQ(56,"Operation not defined");
  /* TaoFunctionReturn(1); */
}

#undef __FUNCT__
#define __FUNCT__ "TaoLMVMMat::ShiftDiagonal"
int TaoLMVMMat::ShiftDiagonal(double c){
  TaoFunctionBegin;
  SETERRQ(56,"Operation not defined");
  /* TaoFunctionReturn(1); */
}

#undef __FUNCT__
#define __FUNCT__ "TaoLMVMMat::extractSubMatrix"
int TaoLMVMMat::extractSubMatrix(TaoMat**M, TaoIndexSet* row,TaoIndexSet* col){
  TaoFunctionBegin;
  SETERRQ(56,"Operation not defined");
  /* TaoFunctionReturn(1); */
};

#undef __FUNCT__
#define __FUNCT__ "TaoLMVMMat::RowScale"
int TaoLMVMMat::RowScale(TaoVec* tv){
  TaoFunctionBegin;
  SETERRQ(56,"Operation not defined");
  /* TaoFunctionReturn(1); */
}

#undef __FUNCT__
#define __FUNCT__ "TaoLMVMMat::ColScale"
int TaoLMVMMat::ColScale(TaoVec* tv){
  TaoFunctionBegin;
  SETERRQ(56,"Operation not defined");
  /* TaoFunctionReturn(1); */
}

#undef __FUNCT__
#define __FUNCT__ "TaoLMVMMat::D_Fischer"
int TaoLMVMMat::D_Fischer(TaoVec *tv1, TaoVec *tv2, TaoVec *tv3, TaoVec *tv4, 
			  TaoVec *tv5, TaoVec *tv6, TaoVec *tv7, TaoVec *tv8)
{
  TaoFunctionBegin;
  SETERRQ(56,"Operation not defined");
  /* TaoFunctionReturn(1); */
}

#undef __FUNCT__
#define __FUNCT__ "TaoLMVMMat::D_SFischer"
int TaoLMVMMat::D_SFischer(TaoVec *tv1, TaoVec *tv2, TaoVec *tv3, TaoVec *tv4, 
			   double, TaoVec *tv5, TaoVec *tv6, 
                           TaoVec *tv7, TaoVec *tv8, TaoVec *tv9)
{
  TaoFunctionBegin;
  SETERRQ(56,"Operation not defined");
  /* TaoFunctionReturn(1); */
}

#undef __FUNCT__
#define __FUNCT__ "TaoLMVMMat::minQuadraticTrustRegion"
int TaoLMVMMat::minQuadraticTrustRegion(TaoVec*,TaoVec*,double, TaoTruth*){
  TaoFunctionBegin;
  SETERRQ(56,"Operation not defined");
  /* TaoFunctionReturn(1); */
}

