/*$Id: blmvm.c 1.87 05/05/09 11:09:17-05:00 sarich@zorak.(none) $*/

#include "blmvm.h"      /*I "tao_solver.h" I*/
static int TaoSetDown_BLMVM(TAO_SOLVER, void*);

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

#undef __FUNCT__  
#define __FUNCT__ "TaoSetUp_BLMVM"
int TaoSetUp_BLMVM(TAO_SOLVER tao, void*solver){

  TAO_BLMVM *blmvm = (TAO_BLMVM *)solver;
  TaoVec*  X;
  int info, lm;

  TaoFunctionBegin;

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

  info = X->Clone(&blmvm->DX); CHKERRQ(info);
  info = X->Clone(&blmvm->GP); CHKERRQ(info);
  info = X->Clone(&blmvm->G); CHKERRQ(info);
  info = X->Clone(&blmvm->XL); CHKERRQ(info);
  info = X->Clone(&blmvm->XU); CHKERRQ(info);

  info = TaoSetLagrangianGradientVector(tao,blmvm->GP);CHKERRQ(info);
  info = TaoSetStepDirectionVector(tao,blmvm->DX);CHKERRQ(info);
  info = TaoSetVariableBounds(tao,blmvm->XL,blmvm->XU);CHKERRQ(info);
  
  info = TaoLMVMGetSize(tao,&lm); CHKERRQ(info);
  blmvm->M = new TaoLMVMMat(X,lm);

  TaoFunctionReturn(0);
}

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

  TaoFunctionBegin;

  info=TaoVecDestroy(blmvm->DX);CHKERRQ(info);blmvm->DX=0;
  info=TaoVecDestroy(blmvm->GP);CHKERRQ(info);blmvm->GP=0;
  info=TaoVecDestroy(blmvm->G);CHKERRQ(info);
  info=TaoVecDestroy(blmvm->XL);CHKERRQ(info);
  info=TaoVecDestroy(blmvm->XU);CHKERRQ(info);
  info=TaoMatDestroy(blmvm->M);CHKERRQ(info);blmvm->M=0; 
  
  TaoFunctionReturn(0);
}

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

  TAO_BLMVM *blmvm = (TAO_BLMVM *) solver;
  TaoVec  *X,*G=blmvm->G,*GP=blmvm->GP; 	/* solution TaoVec */
  TaoVec  *DX=blmvm->DX;
  TaoVec  *XL=blmvm->XL,*XU=blmvm->XU;
  int       info,lsflag,iter=0;
  double    f,gnorm,step=0,gdx;
  TaoTerminateReason reason;
  TaoTruth success;

  TaoFunctionBegin;

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

  info = TaoEvaluateVariableBounds(tao,XL,XU); CHKERRQ(info);
  info = X->Median(XL,X,XU); CHKERRQ(info);

  info = TaoComputeMeritFunctionGradient(tao,X,&f,G); CHKERRQ(info);
  blmvm->pgits=0;

  /* Enter loop */
  while (1){

    info = GP->BoundGradientProjection(G,XL,X,XU); CHKERRQ(info);
    info = GP->Norm2(&gnorm); CHKERRQ(info);

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

    info = blmvm->M->update(X,GP); CHKERRQ(info);
    info = blmvm->M->Solve(G,DX,&success); CHKERRQ(info);
    info = GP->BoundGradientProjection(DX,XL,X,XU); CHKERRQ(info);

    info = GP->Dot(G,&gdx); CHKERRQ(info);
    if (gdx<=0){
      info = DX->ScaleCopyFrom(-1.0,G); CHKERRQ(info);
      gdx=-gnorm*gnorm;
      blmvm->pgits++;
    } else {
      info = DX->Negate(); CHKERRQ(info);
      gdx*= -1;
    }

    /* Line Search */
    if (iter==1){ 
      info = TaoGetInitialTrustRegionRadius(tao,&step); CHKERRQ(info);
      if (step <= 0 ) step=1.0; 
    } else step=1.0;

    info = TaoLineSearchApply(tao,X,G,DX,GP,&f,&step,&lsflag);
    CHKERRQ(info)

  }

  TaoFunctionReturn(0);
}


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

  TaoFunctionBegin;
  info = TaoOptionsHead("Variable metric method for bound constrained optimization");CHKERRQ(info);
  info = TaoLineSearchSetFromOptions(tao);CHKERRQ(info);
  info = TaoOptionsTail();CHKERRQ(info);
  TaoFunctionReturn(0);
}


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

  TaoFunctionBegin;
  
  info=TaoPrintInt(tao,"  Projected gradient iterates: %d\n",blmvm->pgits); CHKERRQ(info);
  info=TaoPrintInt(tao,"  Rejected matrix updates: %d\n",blmvm->M->countrejects());
  CHKERRQ(info);
  
  info = TaoLineSearchView(tao);CHKERRQ(info);

  TaoFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "TaoGetDualVariables_BLMVM" 
static int TaoGetDualVariables_BLMVM(TAO_SOLVER tao, TaoVec *DXL, TaoVec* DXU, void* solver)
{
  TAO_BLMVM *blmvm = (TAO_BLMVM *) solver;
  TaoVec  *G=blmvm->G,*GP=blmvm->GP;
  int       info;

  TaoFunctionBegin;
  
  info = DXL->Waxpby(-1.0,G,1.0,GP); CHKERRQ(info);
  info = DXU->SetToZero(); CHKERRQ(info);
  info = DXL->PointwiseMaximum(DXL,DXU); CHKERRQ(info);

  info = DXU->Waxpby(-1.0,GP,1.0,G); CHKERRQ(info);
  info = DXU->Axpy(1.0,DXL); CHKERRQ(info);

  TaoFunctionReturn(0);
}


/* ---------------------------------------------------------- */
EXTERN_C_BEGIN
#undef __FUNCT__  
#define __FUNCT__ "TaoCreate_BLMVM"
int TaoCreate_BLMVM(TAO_SOLVER tao)
{
  TAO_BLMVM *blmvm;
  int       info;

  TaoFunctionBegin;

  info = TaoNew(TAO_BLMVM,&blmvm); CHKERRQ(info);
  info = PetscLogObjectMemory(tao,sizeof(TAO_BLMVM)); CHKERRQ(info);

  info=TaoSetTaoSolveRoutine(tao,TaoSolve_BLMVM,(void*)blmvm); CHKERRQ(info);
  info=TaoSetTaoSetUpDownRoutines(tao,TaoSetUp_BLMVM,TaoSetDown_BLMVM); CHKERRQ(info);
  info=TaoSetTaoOptionsRoutine(tao,TaoSetOptions_BLMVM); CHKERRQ(info);
  info=TaoSetTaoViewRoutine(tao,TaoView_BLMVM); CHKERRQ(info);
  info=TaoSetTaoDualVariablesRoutine(tao,TaoGetDualVariables_BLMVM); CHKERRQ(info);

  info = TaoSetMaximumIterates(tao,10000); CHKERRQ(info);
  info = TaoSetMaximumFunctionEvaluations(tao,10000); CHKERRQ(info);
  info = TaoSetTolerances(tao,1e-4,1e-4,0,0); CHKERRQ(info);
  info = TaoLMVMSetSize(tao,5); CHKERRQ(info);
  info = TaoSetTrustRegionRadius(tao,0.01);CHKERRQ(info);

  //  info = TaoCreateProjectedLineSearch(tao); CHKERRQ(info);
  info = TaoCreateMoreThuenteBoundLineSearch(tao,0,0.99); CHKERRQ(info);

  blmvm->DX=0;blmvm->GP=0;blmvm->M=0;
  TaoFunctionReturn(0);
}
EXTERN_C_END

