SOCIPBlock.hxx

Go to the documentation of this file.

#ifndef CONICBUNDLE_SOCIPBLOCK_HXX
#define CONICBUNDLE_SOCIPBLOCK_HXX

#include "InteriorPointBlock.hxx"

namespace ConicBundle {


class SOCIPBlock: public virtual InteriorPointBlock
{
protected:
  CH_Matrix_Classes::Integer vecdim; 
  CH_Matrix_Classes::Matrix x;     
  CH_Matrix_Classes::Matrix z;     
  CH_Matrix_Classes::Matrix dx;    
  CH_Matrix_Classes::Matrix dz;    

  CH_Matrix_Classes::Real gammaxsqr;  
  CH_Matrix_Classes::Real gammazsqr;  
  CH_Matrix_Classes::Real omega;      
  CH_Matrix_Classes::Matrix f;        
  CH_Matrix_Classes::Matrix scaled_point; 
  CH_Matrix_Classes::Matrix compl_rhs;  
  
  CH_Matrix_Classes::Real last_rhs_mu;  
  CH_Matrix_Classes::Real mu;           
  CH_Matrix_Classes::Real old_mu;       

  CH_Matrix_Classes::Real last_alpha;   
  
  CH_Matrix_Classes::Matrix oldx;       
  CH_Matrix_Classes::Matrix oldz;       
  
  mutable CH_Matrix_Classes::Matrix tmpvec; 
  mutable CH_Matrix_Classes::Matrix tmpmat; 
  mutable CH_Matrix_Classes::Matrix tmp_xdzpdxz; 

  void point_changed();

  CH_Matrix_Classes::Matrix& apply_Arw(const CH_Matrix_Classes::Matrix& x,CH_Matrix_Classes::Matrix& v) const;

  CH_Matrix_Classes::Matrix& apply_Arwinv(const CH_Matrix_Classes::Matrix& x,CH_Matrix_Classes::Matrix& v) const;

  CH_Matrix_Classes::Matrix& apply_F(CH_Matrix_Classes::Matrix& v) const;

  CH_Matrix_Classes::Matrix& apply_Finv(CH_Matrix_Classes::Matrix& v) const;

  CH_Matrix_Classes::Matrix& apply_Fsqr(CH_Matrix_Classes::Matrix& v,bool minus=false) const;

  int apply_Fsqr(CH_Matrix_Classes::Real* vp,bool minus=false) const;

  CH_Matrix_Classes::Matrix& apply_Finvsqr(CH_Matrix_Classes::Matrix& v,bool minus=false) const;

  int apply_Finvsqr(CH_Matrix_Classes::Real *vp,bool minus=false) const;

  int compute_NTscaling(void);
  
public:
  virtual void clear(CH_Matrix_Classes::Integer dim=0);

  SOCIPBlock(CH_Matrix_Classes::Integer dim=0, CBout* cb=0,int cbinc=-1);

  ~SOCIPBlock();

  virtual CH_Matrix_Classes::Integer get_vecdim() const;
  
  virtual int center_x(CH_Matrix_Classes::Real val,bool add=false);

  virtual int center_z(CH_Matrix_Classes::Real val,bool add=false);

  virtual int set_x(const CH_Matrix_Classes::Matrix& vec,CH_Matrix_Classes::Integer startindex,CH_Matrix_Classes::Real& add_center_value);

  virtual int set_z(const CH_Matrix_Classes::Matrix& vec,CH_Matrix_Classes::Integer startindex,CH_Matrix_Classes::Real& add_center_value);

  virtual int vecgetsax(CH_Matrix_Classes::Matrix& vec,
                        CH_Matrix_Classes::Integer startindex,
                        CH_Matrix_Classes::Real a=1.,
                        bool add=false);

  virtual int vecgetsaz(CH_Matrix_Classes::Matrix& vec,
                        CH_Matrix_Classes::Integer startindex,
                        CH_Matrix_Classes::Real a=1.,
                        bool add=false);

  virtual int get_mu_info(CH_Matrix_Classes::Integer& mudim,
                          CH_Matrix_Classes::Real& tr_xz,
                          CH_Matrix_Classes::Real& tr_xdzpdxz,
                          CH_Matrix_Classes::Real& tr_dxdz,
                          CH_Matrix_Classes::Real& min_xz,
                          CH_Matrix_Classes::Real& max_xz) const;

  virtual int get_nbh_info(CH_Matrix_Classes::Integer mudim,
                           CH_Matrix_Classes::Real tr_xz,
                           CH_Matrix_Classes::Real tr_xdzpdxz,
                           CH_Matrix_Classes::Real tr_dxdz,
                           CH_Matrix_Classes::Real nbh_ubnd,
                           CH_Matrix_Classes::Real& alpha,
                           CH_Matrix_Classes::Real& max_nbh,
                           CH_Matrix_Classes::Real& nrmsqr_xz,
                           CH_Matrix_Classes::Real& nrmsqr_xdzpdxz,
                           CH_Matrix_Classes::Real& nrmsqr_dxdz,
                           CH_Matrix_Classes::Real& ip_xz_xdzpdxz,
                           CH_Matrix_Classes::Real& ip_xz_dxdz,
                           CH_Matrix_Classes::Real& ip_dxdz_xdzpdxz) const;

  virtual int linesearch(CH_Matrix_Classes::Real& alpha) const;

  virtual int add_muxinv(CH_Matrix_Classes::Matrix& rhs,
                         CH_Matrix_Classes::Integer startindex,
                         CH_Matrix_Classes::Real rhsmu,
                         CH_Matrix_Classes::Real rhscorr,
                         bool minus=false);
  
  
  virtual int set_dx(const CH_Matrix_Classes::Matrix& rhs,CH_Matrix_Classes::Integer startindex);

  virtual int set_dx_xizsolverhs(const CH_Matrix_Classes::Matrix& rhs,CH_Matrix_Classes::Integer startindex);


  virtual int apply_xizinv(CH_Matrix_Classes::Matrix& rhs,
                           CH_Matrix_Classes::Integer startindex,
                           bool minus=false);

  virtual int apply_xiz(CH_Matrix_Classes::Matrix& rhs,
                        CH_Matrix_Classes::Integer startindex,
                        bool minus=false);

  virtual int do_step(CH_Matrix_Classes::Real alpha);

  virtual int add_AxizinvAt(const CH_Matrix_Classes::Matrix& A,
                            CH_Matrix_Classes::Symmatrix& globalsys,
                            bool minus=false,
                            bool Atrans=false);

  virtual int add_xiz(CH_Matrix_Classes::Symmatrix& globalsys,CH_Matrix_Classes::Integer startindex,bool minus=false);



  //---------------------- mainly for testing

  virtual int get_vecx(CH_Matrix_Classes::Matrix& vecx,CH_Matrix_Classes::Integer startindex)
  { CH_Matrix_Classes::mat_xey(vecdim,vecx.get_store()+startindex,x.get_store());return 0; }
  
  virtual int get_vecz(CH_Matrix_Classes::Matrix& vecz,CH_Matrix_Classes::Integer startindex)
  { CH_Matrix_Classes::mat_xey(vecdim,vecz.get_store()+startindex,z.get_store());return 0; }
  
  virtual int get_vecdx(CH_Matrix_Classes::Matrix& vecdx,CH_Matrix_Classes::Integer startindex)
  {
    if (dx.dim()!=vecdim)
      return 1;
    CH_Matrix_Classes::mat_xey(vecdim,vecdx.get_store()+startindex,dx.get_store());
    return 0;
  }
  
  
  virtual int get_vecdz(CH_Matrix_Classes::Matrix& vecdz,CH_Matrix_Classes::Integer startindex)
  {
    if (dz.dim()!=vecdim)
      return 1;
    CH_Matrix_Classes::mat_xey(vecdim,vecdz.get_store()+startindex,dz.get_store());
    return 0;
  }

};

  

}

#endif