UQPSolver.hxx

Go to the documentation of this file.

#ifndef CONICBUNDLE_UQPSOLVER_HXX
#define CONICBUNDLE_UQPSOLVER_HXX

#include <vector>
#include "QPSolverObject.hxx"
#include "UQPModelBlock.hxx"
#include "clock.hxx"
#include "Groundset.hxx"

namespace ConicBundle {


  class UQPSolver: public virtual UQPModelPointer, public virtual QPSolverObject
{
 private:
  BundleProxObject* Hp;
  const CH_Matrix_Classes::Matrix* center_yp;
  MinorantPointer gs_aggr;  
  
  //--- problem desription
  CH_Matrix_Classes::Symmatrix Q;  
  CH_Matrix_Classes::Matrix c;     
  CH_Matrix_Classes::Real offset;  

  CH_Matrix_Classes::Matrix A;     
  CH_Matrix_Classes::Matrix b;     

  //--- termination parameters
  CH_Matrix_Classes::Real lowerbound;  
  CH_Matrix_Classes::Real upperbound; 
  CH_Matrix_Classes::Real termeps;   
  CH_Matrix_Classes::Integer maxiter; 

  //--- global variables
  CH_Matrix_Classes::Matrix x;         
  CH_Matrix_Classes::Matrix y;         
  CH_Matrix_Classes::Matrix dx;        
  CH_Matrix_Classes::Matrix dy;        

  CH_Matrix_Classes::Real primalval;   
  CH_Matrix_Classes::Real dualval;     

  CH_Matrix_Classes::Real mu;          
  CH_Matrix_Classes::Matrix Qx;        
  CH_Matrix_Classes::Matrix rhs;       

  bool run_starting_point;         

  CH_Matrix_Classes::Integer iter;     
  CH_Matrix_Classes::Integer status;   

  //--- temporary variables, global only for memory managment purposes
  CH_Matrix_Classes::Symmatrix Qplus;   
  CH_Matrix_Classes::Matrix LinvAt;     
  CH_Matrix_Classes::Symmatrix sysdy;   
  CH_Matrix_Classes::Matrix rd;         
  CH_Matrix_Classes::Matrix xcorr;      
  mutable CH_Matrix_Classes::Matrix tmpvec;  
  mutable CH_Matrix_Classes::Matrix tmpxvec; 
  mutable CH_Matrix_Classes::Matrix tmpyvec; 


  //--- statistics
  CH_Matrix_Classes::Integer sum_iter;    
  CH_Matrix_Classes::Integer sum_choliter; 
  CH_Tools::Clock clock; 
  CH_Tools::Microseconds sum_choltime;  
  CH_Tools::Microseconds QPcoeff_time;  
  CH_Tools::Microseconds QPsolve_time; 

  
  //--- private routines
  void select_new_mu(const CH_Matrix_Classes::Matrix& dx,
                     const CH_Matrix_Classes::Matrix& dy,
                     const CH_Matrix_Classes::Matrix& rhs_residual);
  int predcorr_step(CH_Matrix_Classes::Real& alpha);
  int iterate();

 public:
  void clear();

  void set_defaults(); 

  UQPSolver(CBout* cb=0,int cbinc=-1):UQPModelPointer(cb,cbinc)
  {clear();set_defaults();}
  ~UQPSolver(){}

  
  //----  proper own routines
  void init_size(CH_Matrix_Classes::Integer maxdim){Q.newsize(maxdim); x.newsize(maxdim,1);}

  const CH_Matrix_Classes::Symmatrix& get_Q(void) const {return Q;}
  const CH_Matrix_Classes::Matrix& get_c(void) const {return c;}
  CH_Matrix_Classes::Real get_offset(void) const {return offset;}

  void set_termbounds(CH_Matrix_Classes::Real lb,CH_Matrix_Classes::Real ub){lowerbound=lb;upperbound=ub;}
  void set_termeps(CH_Matrix_Classes::Real te){termeps=te;}
  void set_maxiter(CH_Matrix_Classes::Integer mi){maxiter=mi;}

  CH_Matrix_Classes::Integer get_iter() const {return iter;}
  CH_Matrix_Classes::Integer get_status() const {return status;}
  CH_Matrix_Classes::Real get_termeps() const {return termeps;}
  CH_Matrix_Classes::Integer get_maxiter() const {return maxiter;}

  CH_Matrix_Classes::Real get_primalval()const {return primalval;}
  CH_Matrix_Classes::Real get_dualval()const {return dualval; }

  const CH_Matrix_Classes::Matrix& get_x() {return x;}
  const CH_Matrix_Classes::Matrix& get_y() {return y;}

  int solve(const CH_Matrix_Classes::Symmatrix& Q,const CH_Matrix_Classes::Matrix& c,CH_Matrix_Classes::Real offset);
  // returns status
  //   0 ... if solve terminated correctly for the standard stopping criterion
  //   1 ... if terminated by maxiter or precision tolerances
  //   2 ... if the factorization of the system failed
  //   3 ... if one of the blocks failed 
  
  int resolve();
  // returns status
  //   0 ... if solve terminated correctly for the standard stopping criterion
  //   1 ... if terminated by maxiter or precision tolerances
  //   2 ... if the factorization of the system failed
  //   3 ... if one of the blocks failed 

  int update(const CH_Matrix_Classes::Symmatrix& dQ,const CH_Matrix_Classes::Matrix& dc,CH_Matrix_Classes::Real doffset);
  // solves the system for c+=dc offset+=doffset and returns status
  //   0 ... if solve terminated correctly for the standard stopping criterion
  //   1 ... if terminated by maxiter or precision tolerances
  //   2 ... if the factorization of the system failed
  //   3 ... if one of the blocks failed 

  std::ostream& print_statistics(std::ostream& out) const;

  std::ostream& save(std::ostream& out) const;

  std::istream& restore(std::istream& in);

  //------------  QPSolverOject routines

  void QPclear(){clear();}

  int QPset_parameters(QPSolverParametersObject* /* params */)
  {
    //not used here
    return 0;
  }


  int apply_modification(const GroundsetModification& /* mdf */)
  {return 0;}

  bool QPsupports_yfixing()
  { return true; }

  virtual bool QPsupports_updates()
  { return true;}

  CH_Matrix_Classes::Real QPget_lower_bound()
  { return get_primalval(); }

  int QPsolve(const CH_Matrix_Classes::Matrix& center_y,
              CH_Matrix_Classes::Real lower_bound,
              CH_Matrix_Classes::Real upper_bound,
              CH_Matrix_Classes::Real relprec,
              QPSolverProxObject* Hp,
              const MinorantPointer& gs_aggr,
              CH_Matrix_Classes::Indexmatrix* yfixed);

  int QPupdate(const CH_Matrix_Classes::Matrix& center_y,
               CH_Matrix_Classes::Real lower_bound,
               CH_Matrix_Classes::Real upper_bound,
               CH_Matrix_Classes::Real relprec,
               QPSolverProxObject* Hp,
               const MinorantPointer& gs_aggr,
               CH_Matrix_Classes::Indexmatrix* yfixed,
               const MinorantPointer& delta_gs_subg,
               const CH_Matrix_Classes::Indexmatrix& delta_index);

  int QPresolve(CH_Matrix_Classes::Real lower_bound,
                CH_Matrix_Classes::Real upper_bound,
                CH_Matrix_Classes::Real relprec);

  int QPget_solution(CH_Matrix_Classes::Real&  augval_lb,
                     CH_Matrix_Classes::Real&  augval_ub,
                     CH_Matrix_Classes::Matrix& new_point,
                     CH_Matrix_Classes::Real&  gsaggr_offset,
                     CH_Matrix_Classes::Matrix& gsaggr_gradient);

  bool QPprefer_UQPSolver(QPSolverProxObject*) const
  { return true; }  
  
  bool QPconstrained() const
  {return false;}

  virtual GroundsetModification* QPstart_modification()
  {return 0;}

  int QPapply_modification(const GroundsetModification& /* mdf */)
  {return 0;}

  bool QPis_feasible(const CH_Matrix_Classes::Matrix& /* y */,
                     CH_Matrix_Classes::Real /* relprec */=1e-10)
  {return true;}

  virtual int QPensure_feasibility(CH_Matrix_Classes::Matrix& /* y */,
                                   bool& ychanged,
                                   QPSolverProxObject* /* Hp */,
                                   CH_Matrix_Classes::Real /* relprec */=1e-10)
  { ychanged=false; return 0;}

  bool QPboxconstrained(const CH_Matrix_Classes::Matrix*& lb,
                        const CH_Matrix_Classes::Matrix*& ub,
                        const CH_Matrix_Classes::Indexmatrix*& lbind,
                        const CH_Matrix_Classes::Indexmatrix*& ubind) const
  { lb=ub=0; lbind=ubind=0; return false;}

    std::ostream& QPprint_statistics(std::ostream& out,int /* printlevel */ =0)
  {return print_statistics(out);}
  

};


  
}

#endif