sparsmat.hxx

Go to the documentation of this file.

#ifndef CH_MATRIX_CLASSES__SPARSMAT_HXX
#define CH_MATRIX_CLASSES__SPARSMAT_HXX

#ifndef CH_MATRIX_CLASSES__SYMMAT_HXX
#include "symmat.hxx"
#endif

namespace CH_Matrix_Classes {

//everything involving a "Sparsesym" is implemented in "sparssym.cxx/hxx"

#define SPARSE_ZERO_TOL 1e-60


class Sparsemat: protected Memarrayuser
{
    friend class Indexmatrix;
    friend class Matrix;
    friend class Symmatrix;
    friend class Sparsesym;
    
private:

    Mtype mtype;   
    Integer nr,    
            nc;    

    //column by column representation
    Indexmatrix colinfo;  
    Indexmatrix colindex; 
    Matrix colval;        

    //row by row representation (see column by column)
    Indexmatrix rowinfo; 
    Indexmatrix rowindex;
    Matrix rowval;       

        
    Real tol;          

    bool is_init;      

    void init_to_zero();          
    Integer find_column(Integer i) const;
    
    Integer find_row(Integer i) const;

public:

  //----------------------------------------------------
  //----  constructors, destructor, and initialization
  //----------------------------------------------------


  inline Sparsemat();                      
  inline Sparsemat(const Sparsemat& A,Real d=1.);       
  inline Sparsemat(Integer nr,Integer nc);  
  inline Sparsemat(Integer nr,Integer nc,Integer nz,
                   const Integer *ini,const Integer *inj,const Real* va);
  inline Sparsemat(Integer nr,Integer nc,Integer nz,
                   const Indexmatrix& ini,const Indexmatrix& inj, const Matrix& va);

#if (CONICBUNDLE_DEBUG>=1)
  void set_init(bool i){is_init=i;}
  bool get_init() const {return is_init;} 
#else 
  void set_init(bool /* i */){}
  bool get_init() const {return true;}
#endif 

  inline Sparsemat& init(const Sparsemat& A,Real d=1.);       
  inline Sparsemat& init(const Matrix& A,Real d=1.);          
  inline Sparsemat& init(const Indexmatrix& A,Real d=1.);     
  inline Sparsemat& init(const Symmatrix&,Real d=1.);     
  inline Sparsemat& init(const Sparsesym&,Real d=1.);    
  inline Sparsemat& init(Integer nr,Integer nc);  //zero-matrix of size nr*nc
  Sparsemat& init(Integer nr,Integer nc,Integer nz,
                  const Integer *ini,const Integer *inj,const Real* va);
  Sparsemat& init(Integer nr,Integer nc,Integer nz,
                  const Indexmatrix& ini,const Indexmatrix& inj, const Matrix& va);

  void set_tol(Real t){tol=t;}



  inline Sparsemat(const Matrix& A, Real d=1.);          
  inline Sparsemat(const Indexmatrix& A,Real d=1.);
  inline Sparsemat(const Symmatrix& A,Real d=1.); 
  inline Sparsemat(const Sparsesym& A,Real d=1.);


  //----------------------------------------------------
  //----  size and type information
  //----------------------------------------------------


  void dim(Integer& in_nr, Integer& in_nc) const {in_nr=nr; in_nc=nc;}

  Integer dim() const {return nr*nc;}

  Integer rowdim() const {return nr;}

  Integer coldim() const {return nc;}

  Integer nonzeros() const {return colval.dim();}

  Integer col_nonzeros(Integer i,Integer* startind=0) const;

  Integer row_nonzeros(Integer i,Integer* startind=0) const;

  Mtype get_mtype() const {return mtype;}

  


    
  //--------------------------------
  //----  Indexing and Submatrices
  //--------------------------------


    
  Real operator()(Integer i,Integer j) const;     

  inline Real operator()(Integer i) const;           

  inline Real operator[](Integer i) const; 

  Sparsemat col(Integer i) const;

  Sparsemat row(Integer i) const;
  Sparsemat cols(const Indexmatrix& ind) const;
  Sparsemat rows(const Indexmatrix& ind) const;
  
  Sparsemat& delete_rows(const Indexmatrix& ind);
  Sparsemat& delete_cols(const Indexmatrix& ind);

  Sparsemat& insert_row(Integer i,const Sparsemat& v); 
  Sparsemat& insert_col(Integer i,const Sparsemat& v); 

  Sparsemat& concat_right(const Sparsemat& A);
  Sparsemat& concat_below(const Sparsemat& A);

  const Indexmatrix& get_colinfo() const {return colinfo;}
  const Indexmatrix& get_colindex() const {return colindex;}
  const Matrix& get_colval() const {return colval;}

  const Indexmatrix& get_rowinfo() const {return rowinfo;}
  const Indexmatrix& get_rowindex() const {return rowindex;}
  const Matrix& get_rowval() const {return rowval;}
  
  void get_edge_rep(Indexmatrix& I, Indexmatrix& J, Matrix& val) const;
  int get_edge(Integer i,Integer& indi,Integer& indj, Real& val) const;
  int contains_support(const Sparsemat& A) const;





  friend inline Sparsemat concat_right(const Sparsemat& A,const Sparsemat& B);
  friend inline Sparsemat concat_below(const Sparsemat& A,const Sparsemat& B);

  friend void swap(Sparsemat& A, Sparsemat& B);


  //------------------------------
  //----  BLAS-like Routines
  //------------------------------


  Sparsemat& xeya(const Sparsemat& A,Real d=1.);
  Sparsemat& xeya(const Matrix& A,Real d=1.);
  Sparsemat& xeya(const Indexmatrix& A,Real d=1.);



  friend Sparsemat& xbpeya(Sparsemat& x,const Sparsemat& y,Real alpha,Real beta,int ytrans);
  
  friend Matrix& genmult(const Sparsemat& A,const Matrix& B,Matrix &C,
                         Real alpha,Real beta, int atrans,int btrans);
      //C=beta*C+alpha A B

  friend Matrix& genmult(const Sparsemat& A,const Matrix& B,Integer colB,Matrix &C,Integer colC,Real alpha,Real beta, int atrans,int btrans);
      //C=beta*C+alpha A B

  friend Matrix& genmult(const Matrix& A,const Sparsemat& B,Matrix &C,
                         Real alpha,Real beta, int atrans,int btrans);
      //C=beta*C+alpha A B
  
  friend Matrix& genmult(const Sparsemat& A,const Sparsemat& B,Matrix &C,
                         Real alpha,Real beta, int atrans,int btrans);
      //C=beta*C+alpha A B


  //------------------------------
  //----  usual operators
  //------------------------------


  inline Sparsemat& operator=(const Sparsemat& A);
  inline Sparsemat& operator+=(const Sparsemat& A);
  inline Sparsemat& operator-=(const Sparsemat& A);
  inline Sparsemat operator-() const;
  
  inline Sparsemat& operator*=(Real d);
  inline Sparsemat& operator/=(Real d);
  
  Sparsemat& operator=(const Matrix& A);
  Sparsemat& operator=(const Indexmatrix& A);

  Sparsemat& transpose();          
   

  
  friend Sparsemat operator*(const Sparsemat& A,const Sparsemat& B);
  friend inline Sparsemat operator+(const Sparsemat &A,const Sparsemat& B);
  friend inline Sparsemat operator-(const Sparsemat &A,const Sparsemat& B);

  friend inline Sparsemat operator*(const Sparsemat& A,Real d);
  friend inline Sparsemat operator*(Real d,const Sparsemat& A);
  friend inline Sparsemat operator/(const Sparsemat& A,Real d);
  
  friend inline Matrix operator*(const Sparsemat& A,const Matrix& B);
  friend inline Matrix operator*(const Matrix& A,const Sparsemat& B);
  friend inline Matrix operator+(const Sparsemat& A,const Matrix& B);
  friend inline Matrix operator+(const Matrix& A,const Sparsemat& B);
  friend inline Matrix operator-(const Sparsemat& A,const Matrix& B);
  friend inline Matrix operator-(const Matrix& A,const Sparsemat& B);

  friend inline Sparsemat transpose(const Sparsemat& A);


  //------------------------------------------
  //----  Connections to other Matrix Classes
  //------------------------------------------
  

  Sparsemat& xeya(const Sparsesym& A,Real d=1.);
  inline Sparsemat& operator=(const Symmatrix& A);
  inline Sparsemat& operator=(const Sparsesym& A);



  friend Matrix& genmult(const Symmatrix& A,const Sparsemat& B,Matrix& C,
                         Real alpha,Real beta,int btrans);
      //returns C=beta*C+alpha*A*B, where A and B may be transposed

  friend Matrix& genmult(const Sparsemat& A,const Symmatrix& B,Matrix& C,
                         Real alpha,Real beta, int atrans);

  friend Matrix& genmult(const Sparsesym& A,const Sparsemat& B,Matrix &C,
                         Real alpha,Real beta, int btrans);

  friend Matrix& genmult(const Sparsemat& A,const Sparsesym& B,Matrix &C,
                         Real alpha,Real beta,int atrans);

  friend Symmatrix& rankadd(const Sparsemat& A,Symmatrix& C,
                            Real alpha,Real beta,int trans);
      // returns C=beta*C+alpha* A*A^T, where A may be transposed
    
  friend Symmatrix& scaledrankadd(const Sparsemat& A,const Matrix& D,Symmatrix& C,
                            Real alpha,Real beta,int trans);
      // returns C=beta*C+alpha* A*D*A^T, where A may be transposed
    
  friend Symmatrix& rank2add(const Sparsemat& A,const Matrix& B,Symmatrix& C,
                             Real alpha,Real beta,int trans);
      // returns C=beta*C+alpha* sym(A*B^T) [or sym(A^T*B)]

  friend inline Matrix operator*(const Symmatrix& A,const Sparsemat& B);
  friend inline Matrix operator*(const Sparsemat& A,const Symmatrix& B);


  //------------------------------
  //----  Elementwise Operations
  //------------------------------


  friend Sparsemat abs(const Sparsemat& A);                 


  //----------------------------
  //----  Numerical Methods
  //----------------------------


  Sparsemat& scale_rows(const Matrix& vec);    //A=diag(vec)*A
  Sparsemat& scale_cols(const Matrix& vec);    //A=A*diag(vec)




  friend Real trace(const Sparsemat& A);               //=sum(diag(A))
  friend Real ip(const Sparsemat& A, const Sparsemat& B); //=trace(B^t*A)
  friend Real ip(const Sparsemat& A, const Matrix& B); //=trace(B^t*A)
  friend inline Real ip(const Matrix& A, const Sparsemat& B);
  friend Real colip(const Sparsemat& A,Integer j,const Matrix* scaling);
  friend Real rowip(const Sparsemat& A,Integer i,const Matrix* scaling);
  friend Sparsemat colsip(const Sparsemat& A,const Matrix* scaling);
  friend Sparsemat rowsip(const Sparsemat& A,const Matrix* scaling);
  friend inline Real norm2(const Sparsemat& A);               //=sqrt(ip(A,A));

  friend Matrix sumrows(const Sparsemat& A);   //=(1 1 1 ... 1)*A
  friend Matrix sumcols(const Sparsemat& A);   //=A*(1 1 ... 1)^t
  friend inline Real sum(const Sparsemat& A);  //=(1 1 ... 1)*A*(1 1 ... 1)^t


  //---------------------------------------------
  //----  Comparisons / Max / Min / sort / find
  //---------------------------------------------


  friend int equal(const Sparsemat& A, const Sparsemat& B,Real eqtol);



  //--------------------------------
  //----  Input / Output
  //--------------------------------


    void display(std::ostream& out, 
                 int precision=0,   
                 int width=0,       
                 int screenwidth=0  
                 ) const;



    friend std::ostream& operator<<(std::ostream& o,const Sparsemat &v);
    friend std::istream& operator>>(std::istream& i,Sparsemat &v);


};



// **************************************************************************
//                make non inline friends available outside
// **************************************************************************

  void swap(Sparsemat& A, Sparsemat& B); 


  
  Sparsemat& xbpeya(Sparsemat& x,const Sparsemat& y,Real alpha=1.,Real beta=0.,int ytrans=0);

Matrix& genmult(const Sparsemat& A,const Matrix& B,Matrix &C,
                  Real alpha=1.,Real beta=0., int atrans=0,int btrans=0);
  Matrix& genmult(const Matrix& A,const Sparsemat& B,Matrix &C,
                  Real alpha=1.,Real beta=0., int atrans=0,int btrans=0);
  Matrix& genmult(const Sparsemat& A,const Sparsemat& B,Matrix &C,
                  Real alpha=1.,Real beta=0., int atrans=0,int btrans=0);

  Sparsemat operator*(const Sparsemat& A,const Sparsemat& B);

  
  Matrix& genmult(const Sparsesym& A,const Sparsemat& B,Matrix &C,
                         Real alpha=1.,Real beta=0., int btrans=0);

  Matrix& genmult(const Sparsemat& A,const Sparsesym& B,Matrix &C,
                         Real alpha=1.,Real beta=0.,int atrans=0);

  Sparsemat abs(const Sparsemat& A);
  
  Real trace(const Sparsemat& A);               //=sum(diag(A))

  Real ip(const Sparsemat& A, const Sparsemat& B); //=trace(B^t*A)

  Real ip(const Sparsemat& A, const Matrix& B); //=trace(B^t*A)

  Real colip(const Sparsemat& A,Integer j,const Matrix* scaling=0);

  Real rowip(const Sparsemat& A,Integer i,const Matrix* scaling=0);

  Sparsemat colsip(const Sparsemat& A, const Matrix* colscaling=0);

  Sparsemat rowsip(const Sparsemat& A, const Matrix* rowscaling=0);

  Matrix sumrows(const Sparsemat& A);   //=(1 1 1 ... 1)*A

  Matrix sumcols(const Sparsemat& A);   //=A*(1 1 ... 1)^t  

  std::ostream& operator<<(std::ostream& o,const Sparsemat &v);

  std::istream& operator>>(std::istream& i,Sparsemat &v);


// **************************************************************************
//                   implementation of inline functions
// **************************************************************************

inline void Sparsemat::init_to_zero()
{
 mtype=MTsparse;
 nr=nc=0;
 tol=SPARSE_ZERO_TOL;
 colinfo.init(0,3,Integer(0)); colindex.init(0,1,Integer(0)); colval.init(0,1,0.);
 rowinfo.init(0,3,Integer(0)); rowindex.init(0,1,Integer(0)); rowval.init(0,1,0.);
#if (CONICBUNDLE_DEBUG>=1)
 is_init=true;
#endif
}

//initialize

inline Sparsemat& Sparsemat::init(const Sparsemat& A,Real d)
{ return xeya(A,d);}
inline Sparsemat& Sparsemat::init(const Matrix& A,Real d)
{ return xeya(A,d);}
inline Sparsemat& Sparsemat::init(const Indexmatrix& A,Real d)
{ return xeya(A,d);}


inline Sparsemat& Sparsemat::init(Integer r,Integer c)
{
chk_range(r,c,-1,-1);
 colinfo.init(0,3,Integer(0)); colindex.init(0,1,Integer(0)); colval.init(0,1,0.);
 rowinfo.init(0,3,Integer(0)); rowindex.init(0,1,Integer(0)); rowval.init(0,1,0.);
nr=r; nc=c;
chk_set_init(*this,1);
return *this;
}

inline Sparsemat::Sparsemat()
{ init_to_zero(); chk_set_init(*this,1);}

  inline Sparsemat::Sparsemat(const Sparsemat& A,double d):Memarrayuser()
{ init_to_zero(); xeya(A,d);}
inline Sparsemat::Sparsemat(const Matrix& A,double d)
{ init_to_zero(); xeya(A,d);}
inline Sparsemat::Sparsemat(const Indexmatrix& A,double d)
{ init_to_zero(); xeya(A,d);}

inline Sparsemat::Sparsemat(Integer r,Integer c)
{ init_to_zero(); init(r,c); }

inline Sparsemat::Sparsemat(Integer in_nr,Integer in_nc,Integer nz,
       const Integer *ini,const Integer *inj,const Real* va)
{ init_to_zero(); init(in_nr,in_nc,nz,ini,inj,va);}

inline Sparsemat::Sparsemat(Integer in_nr,Integer in_nc,Integer nz,
              const Indexmatrix& ini,const Indexmatrix& inj, const Matrix& va)
{ init_to_zero(); init(in_nr,in_nc,nz,ini,inj,va);}

inline Real Sparsemat::operator()(Integer i) const
{return (*this)(i%nr,i/nr);}           
inline Real Sparsemat::operator[](Integer i) const
{return (*this)(i%nr,i/nr);}
    
inline Sparsemat& Sparsemat::transpose()
{chk_init(*this);Integer h=nr; nr=nc;nc=h;swap(colinfo,rowinfo);
 swap(colindex,rowindex);swap(colval,rowval);return *this;}          

inline Sparsemat concat_right(const Sparsemat& A,const Sparsemat& B)
{Sparsemat C(A);C.concat_right(B);return C;}

inline Sparsemat concat_below(const Sparsemat& A,const Sparsemat& B)
{Sparsemat C(A);C.concat_below(B);return C;}

inline Real norm2(const Sparsemat& A)
{
 chk_init(A);
 return norm2(A.colval);
}


inline Sparsemat& Sparsemat::operator=(const Sparsemat &A)
{ return xeya(A); }
inline Sparsemat& Sparsemat::operator+=(const Sparsemat &A)
{ return xbpeya(*this,A,1.,1.); }
inline Sparsemat& Sparsemat::operator-=(const Sparsemat &A)
{ return xbpeya(*this,A,-1.,1.); }
inline Sparsemat& Sparsemat::operator*=(Real d)
{chk_init(*this); colval*=d; rowval*=d; return *this;}
inline Sparsemat& Sparsemat::operator/=(Real d)
{chk_init(*this);colval/=d;rowval/=d;return *this;}
inline Sparsemat Sparsemat::operator-() const
{return Sparsemat(*this,-1.); }



inline Sparsemat operator+(const Sparsemat &A,const Sparsemat& B)
{ Sparsemat C(A); return xbpeya(C,B,1.,1.); }

inline Sparsemat operator-(const Sparsemat &A,const Sparsemat& B)
{ Sparsemat C(A); return xbpeya(C,B,-1.,1.); }

inline Sparsemat operator*(const Sparsemat& A,Real d)
{return Sparsemat(A,d);}

inline Sparsemat operator*(Real d,const Sparsemat& A) 
{return Sparsemat(A,d);}

inline Sparsemat operator/(const Sparsemat& A,Real d) 
{return Sparsemat(A,1./d);}

inline Sparsemat& Sparsemat::operator=(const Matrix &A)
{return xeya(A);}
inline Sparsemat& Sparsemat::operator=(const Indexmatrix &A)
{return xeya(A);}


inline Matrix operator*(const Sparsemat& A,const Matrix& B) 
    {Matrix C;return genmult(A,B,C);}

inline Matrix operator*(const Matrix& A,const Sparsemat& B) 
    {Matrix C;return genmult(A,B,C);}

inline Matrix operator+(const Sparsemat& A,const Matrix& B)
    {Matrix C(B);return C.xpeya(A);}

inline Matrix operator+(const Matrix& A,const Sparsemat& B)
    {Matrix C(A);return C.xpeya(B);}

inline Matrix operator-(const Sparsemat& A,const Matrix& B)
    {Matrix C(B,-1.);return C.xpeya(A);}

inline Matrix operator-(const Matrix& A,const Sparsemat& B)
    {Matrix C(A);return C.xpeya(B,-1.);}


inline Sparsemat::Sparsemat(const Symmatrix& A,Real d)
{ init_to_zero(); xeya(Matrix(A),d);}
inline Sparsemat& Sparsemat::init(const Symmatrix &A,Real d)
{return xeya(Matrix(A),d);}
inline Sparsemat& Sparsemat::operator=(const Symmatrix &A)
{return xeya(Matrix(A));}

inline Real ip(const Matrix& A, const Sparsemat& B)
{return ip(B,A);}

inline Real sum(const Sparsemat& A)
{return sum(A.colval);}

inline Sparsemat transpose(const Sparsemat& A)
{Sparsemat B(A); B.transpose(); return B;}


inline Matrix::Matrix(const Sparsemat& A, Real d)        
{init_to_zero(); xeya(A,d);}

inline Matrix& Matrix::init(const Sparsemat& A, Real d)  
{return xeya(A,d);}

inline Matrix& Matrix::operator=(const Sparsemat & A)       
{return xeya(A);}

inline Matrix& Matrix::operator*=(const Sparsemat &A)
{ Matrix C; return *this=genmult(*this,A,C);}

inline Matrix& Matrix::operator+=(const Sparsemat &A)       
{return xpeya(A);}

inline Matrix& Matrix::operator-=(const Sparsemat &A)       
{return xpeya(A,-1.);}



  int equal(const Sparsemat& A, const Sparsemat& B,Real eqtol=1e-10);



}

//in order to define the inlines associated with the other matrix
//classes, the hxx-files of these are included here, as well. 

#ifndef CH_MATRIX_CLASSES__SPARSSYM_HXX
#include "sparssym.hxx"
#endif

#endif