/*$Id: lmvm.c 1.99 05/05/10 15:21:38-05:00 sarich@zorak.(none) $*/

#include "lmvm.h"             /*I "tao_solver.h" I*/
#include "src/matrix/lmvmmat.h"             /*I "tao_solver.h" I*/

static int TaoDestroy_LMVM(TAO_SOLVER, void*);

/* ---  TESTING ------------------------------------------------------- */
#undef __FUNCT__  
#define __FUNCT__ "TaoInitializeLMVMmatrix"
int TaoInitializeLMVMmatrix(TAO_SOLVER tao, TaoVec *HV){
  int info;
  TAO_LMVM *lmvm;
  
  TaoFunctionBegin;
  info = TaoGetSolverContext(tao,"tao_lmvm",(void**)&lmvm); CHKERRQ(info);
  if (lmvm && lmvm->M){
    info = lmvm->M->InitialApproximation( HV );CHKERRQ(info);
  }
  TaoFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "TaoSetUp_LMVM"
int TaoSetUp_LMVM(TAO_SOLVER tao, void* solver)
{
  TAO_LMVM *lmvm = (TAO_LMVM *)solver;
  TaoVec  *xx;
  int    info, lm;

  TaoFunctionBegin;

  info = TaoGetSolution(tao,&xx);CHKERRQ(info);

  info = TaoLMVMGetSize(tao,&lm); CHKERRQ(info);
  lmvm->M = new TaoLMVMMat(xx,lm);

  info = xx->Clone(&lmvm->G); CHKERRQ(info);
  info = xx->Clone(&lmvm->Work); CHKERRQ(info);
  info = xx->Clone(&lmvm->DX); CHKERRQ(info);

  info = TaoSetLagrangianGradientVector(tao,lmvm->G);CHKERRQ(info);
  info = TaoSetStepDirectionVector(tao,lmvm->DX);CHKERRQ(info);

  TaoFunctionReturn(0);
}

/* ---------------------------------------------------------- */
#undef __FUNCT__  
#define __FUNCT__ "TaoDestroy_LMVM"
static int TaoDestroy_LMVM(TAO_SOLVER tao, void*solver)
{
  int info;
  TAO_LMVM *lmvm = (TAO_LMVM *)solver;

  TaoFunctionBegin;

  info=TaoVecDestroy(lmvm->G);CHKERRQ(info);
  info=TaoVecDestroy(lmvm->DX);CHKERRQ(info);
  info=TaoVecDestroy(lmvm->Work);CHKERRQ(info);
  info=TaoMatDestroy(lmvm->M);CHKERRQ(info);

  info = TaoSetLagrangianGradientVector(tao,0);CHKERRQ(info);
  info = TaoSetStepDirectionVector(tao,0);CHKERRQ(info);
  lmvm->DX=0;lmvm->Work=0;lmvm->M=0;

  TaoFunctionReturn(0);
}

/*------------------------------------------------------------*/
#undef __FUNCT__  
#define __FUNCT__ "TaoSolve_LMVM"
static int TaoSolve_LMVM(TAO_SOLVER tao, void *solver)
{
  TAO_LMVM *lmvm = (TAO_LMVM *)solver;
  TaoVec       *xx, *gg=lmvm->G;	 /* gradient vector */
  TaoVec       *dxdx=lmvm->DX, *ww=lmvm->Work;
  int          info,lsflag=0,iter=0;
  double       f, gdx, gnorm, dxnorm, r0, step=0;
  TaoTerminateReason reason;
  TaoTruth     success;

  TaoFunctionBegin;

  info = TaoGetSolution(tao,&xx);CHKERRQ(info);  

  info = TaoComputeMeritFunctionGradient(tao,xx,&f,gg);CHKERRQ(info);
  info = gg->Norm2(&gnorm);  CHKERRQ(info);

  /* Enter loop */
  while (1){

    /* Test for convergence */
    info=TaoMonitor(tao,iter++,f,gnorm,0,step,&reason);CHKERRQ(info);
    if (reason!=TAO_CONTINUE_ITERATING) break;

    /* Compute step direction */
    info = lmvm->M->update(xx,gg); CHKERRQ(info);
    info = lmvm->M->Solve(gg,dxdx,&success); CHKERRQ(info);
    info = dxdx->Negate(); CHKERRQ(info);
    info = dxdx->Dot(gg,&gdx); CHKERRQ(info);
    if (gdx>=0){
      info = dxdx->ScaleCopyFrom(-1.0,gg); CHKERRQ(info);
      gdx = -gnorm*gnorm;
    }

    /* Line Search */
    info = dxdx->Norm2(&dxnorm);  CHKERRQ(info);
    info = TaoGetInitialTrustRegionRadius(tao,&r0); CHKERRQ(info);
    step = TaoMin(1.0,r0/dxnorm);
    if (step < 1.0){
      info = PetscLogInfo((tao,"Restrict step size to %g so the change in solution <= %g\n",step,r0)); CHKERRQ(info);
    }
    info = TaoLineSearchApply(tao,xx,gg,dxdx,ww,&f,&step,&lsflag);
    CHKERRQ(info);
    info = gg->Norm2(&gnorm);CHKERRQ(info);

  }

  TaoFunctionReturn(0);

}


/*------------------------------------------------------------*/
#undef __FUNCT__  
#define __FUNCT__ "TaoSetOptions_LMVM"
static int TaoSetOptions_LMVM(TAO_SOLVER tao, void *solver)
{
  int        info;

  TaoFunctionBegin;
  info = TaoOptionsHead("Limited Memory Variable Metric method for unconstrained optimization");CHKERRQ(info);
  info = TaoLineSearchSetFromOptions(tao);CHKERRQ(info);
  info = TaoOptionsTail();CHKERRQ(info);

  TaoFunctionReturn(0);
}

/* ---------------------------------------------------------- */
#undef __FUNCT__  
#define __FUNCT__ "TaoLMVMGetX0"
int TaoLMVMGetX0(TAO_SOLVER tao,TaoVec *x0)
{
  TAO_LMVM *lmvm;
  int info;
  TaoFunctionBegin;
  TaoValidHeaderSpecific(tao,TAO_COOKIE,1);
  info=TaoGetSolverContext(tao,"tao_lmvm",(void**)&lmvm);CHKERRQ(info);
  if (lmvm && lmvm->M){
    info=lmvm->M->GetX0(x0);CHKERRQ(info);
  }
  TaoFunctionReturn(0);
}

/*------------------------------------------------------------*/
#undef __FUNCT__  
#define __FUNCT__ "TaoView_LMVM"
static int TaoView_LMVM(TAO_SOLVER tao,void *solver)
{
  int        info;

  TaoFunctionBegin;
  info = TaoLineSearchView(tao);CHKERRQ(info);
  TaoFunctionReturn(0);
}

/* ---------------------------------------------------------- */

EXTERN_C_BEGIN
#undef __FUNCT__  
#define __FUNCT__ "TaoCreate_LMVM"
int TaoCreate_LMVM(TAO_SOLVER tao)
{
  TAO_LMVM *lmvm;
  int info;

  TaoFunctionBegin;

  info = TaoNew(TAO_LMVM,&lmvm); CHKERRQ(info);
  info = PetscLogObjectMemory(tao,sizeof(TAO_LMVM)); CHKERRQ(info);

  info=TaoSetTaoSolveRoutine(tao,TaoSolve_LMVM,(void*)lmvm); CHKERRQ(info);
  info=TaoSetTaoSetUpDownRoutines(tao,TaoSetUp_LMVM,TaoDestroy_LMVM); CHKERRQ(info);
  info=TaoSetTaoOptionsRoutine(tao,TaoSetOptions_LMVM); CHKERRQ(info);
  info=TaoSetTaoViewRoutine(tao,TaoView_LMVM); CHKERRQ(info);

  info = TaoSetMaximumIterates(tao,2000); CHKERRQ(info);
  info = TaoSetMaximumFunctionEvaluations(tao,4000); CHKERRQ(info);
  info = TaoSetTolerances(tao,1e-4,1e-4,0,0); CHKERRQ(info);
  info = TaoLMVMSetSize(tao,5); CHKERRQ(info);
  info = TaoSetTrustRegionRadius(tao,1.0e16); CHKERRQ(info);
  info = TaoCreateMoreThuenteLineSearch(tao,0,0.99); CHKERRQ(info);
  lmvm->DX=0;lmvm->Work=0;lmvm->M=0;

  TaoFunctionReturn(0);
}
EXTERN_C_END



