lanczpol.hxx

Go to the documentation of this file.

#ifndef CH_MATRIX_CLASSES__LANCZPOL_HXX
#define CH_MATRIX_CLASSES__LANCZPOL_HXX

#ifndef CH_MATRIX_CLASSES__LANCZOS_HXX
#include "lanczos.hxx"
#endif
#include "clock.hxx"

namespace CH_Matrix_Classes {




class Lanczpol:public Lanczos,protected Memarrayuser
{
private:
  int     ierr;        
  Integer maxop;       
  Integer maxiter;     
  Integer maxguessiter;
  Integer guessmult;   
  Integer choicencheb; 
  Integer nchebit;     
  Integer choicenbmult;
  Integer nblockmult;  
  Integer nlanczvecs;  
  Integer retlanvecs;  
  
  //--- global variables to allow reastarting
  Integer neigfound;  
  Integer blocksz;    
  Integer iter;       
  Integer nmult;      
  
  Real errc;   
  Real eps;    
  Real mcheps; 
  Real maxval; 
  Real minval; 
  Real polval; 
  
  Matrix X; 
  Matrix C; 
  Matrix d; 
  Matrix e; 
  Matrix Xqr; 
  Matrix u; 
  Matrix v; 
  Matrix w; 
  Matrix minvec;  
  
  const Lanczosmatrix* bigmatrix; 
    
  int stop_above;    
  Real upper_bound;  
  
  Integer ncalls;    
  Integer mymaxj;    

  CH_Tools::GB_rand randgen;  

  CH_Tools::Clock myclock;           
  CH_Tools::Microseconds time_mult;  
  CH_Tools::Microseconds time_mult_sum;  
  CH_Tools::Microseconds time_iter; 
  CH_Tools::Microseconds time_sum; 

  int print_level;   
  std::ostream* myout;   

  void error(const char *s){
    if (myout) 
      (*myout)<<"Lanczos Error: "<<s<<std::endl;
  }

  int guess_extremes(Integer nproposed);
    
  int dhua(Integer nreig,Integer nproposed, Integer maxmult);
    
  int bklanc(Integer neigfound,Integer blocksz,Integer s,
               const Matrix &d,Matrix &C,Matrix &X,
               Matrix &e,Matrix &u,Matrix &v);

  int bkqrlanc(Integer neigfound,Integer blocksz,Integer s,
               const Matrix &d,Matrix &C,Matrix &X,
               Matrix &e,Matrix &u,Matrix &v);

  int bklanccheb(Integer neigfound,Integer blocksz,Integer s,
                    Matrix &C,Matrix &X,Matrix &e,Matrix &v);

  int bkqrlanccheb(Integer neigfound,Integer blocksz,Integer s,
                    Matrix &C,Matrix &X,Matrix &e,Matrix &v);
  
  int err(Integer neigfound,Integer blocksz,const Matrix& X,Matrix &e);

  int cnvtst(Integer neigfound,Integer blocksz,Real& errc,Real eps,
               const Matrix& d,const Matrix &e,Integer &nconv);
    
  int eigen(Integer neigfound,Integer blocksz,Integer sbsz,Matrix& C,
            Matrix& d,Matrix& u,Matrix& v,Real& af);

  int sectn(Matrix& X,Integer neigfound,Integer blocksz,
            Matrix& C,Matrix& d,Matrix& u,Matrix& v,
            Real& af);

  int rotate_extremes(Integer neigfound,Integer sbs,Matrix& d,
                      const Matrix& C,Matrix &X,Matrix& v);
  
  int rotate(Integer neigfound,Integer sbs,Integer l,
             const Matrix& C,Matrix &X,Matrix& v);

  int random(Matrix& X,Integer j);

  Integer orthog(Integer offset,Integer blocksz,Matrix &X,Matrix& B);

  int blockcheby(Integer col_offset,const Matrix& X,Matrix& v);

  Real scalarcheby(Real xval);

  int tred2(Integer n,const Matrix& C,Matrix& u,Matrix& v,Matrix& Z);

  int tql2(Integer n,Matrix &u,Matrix &v, Matrix& Z);

public:
  Lanczpol();
  ~Lanczpol();


  
  void set_mineig(Real ie){minval=ie;}
  void set_maxmult(Integer mop){maxop=mop;}
  void set_maxiter(Integer mi){maxiter=mi;}
  void set_relprec(Real relprec){eps=relprec;}  
  void set_nchebit(Integer cheb){choicencheb=cheb;} 
  void set_nblockmult(Integer nb){choicenbmult=nb;}
  
  void enable_stop_above(Real ub){stop_above=1; upper_bound=ub;}
  void disable_stop_above(){stop_above=0;}
  
  void set_retlanvecs(Integer nl){retlanvecs=nl;} 
  
  int get_lanczosvecs(Matrix& val,Matrix& vecs) const;
  
  Real get_relprec(void){return eps;} //relative precision
  
  int get_err() const{return ierr;}
  Integer get_iter() const{return iter;}
  Integer get_nmult() const{return nmult;}

    
    int compute(const Lanczosmatrix* bigmat, 
                Matrix& eigval,        
                Matrix& eigvec,        
                Integer nreig,         
                Integer in_blocksz=0,  
                Integer maxcol=0       
                );



  
  void set_out(std::ostream* o=0,int pril=1)
  {myout=o; print_level=pril;}


  std::ostream& save(std::ostream& out) const;
  
  std::istream& restore(std::istream& in);
};


}

#endif