sparssym.hxx

Go to the documentation of this file.

#ifndef CH_MATRIX_CLASSES__SPARSSYM_HXX
#define CH_MATRIX_CLASSES__SPARSSYM_HXX

#ifndef CH_MATRIX_CLASSES__SPARSMAT_HXX
#include "sparsmat.hxx"
#endif

namespace CH_Matrix_Classes {

// **************************************************************************
//                            class definition
// **************************************************************************


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

  Mtype mtype;   
  Integer nr;    
  
  //column by column representation of lower triangle
  Indexmatrix colinfo;  
  Indexmatrix colindex; 
  Matrix colval;   
  Indexmatrix suppind;  
  Indexmatrix suppcol;  
        
    Real tol;          

    bool is_init;      

    inline void init_to_zero();            

    void update_support();     

public:

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


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

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

  inline Sparsesym& init(const Sparsesym&,Real d=1.);   
  inline Sparsesym& init(const Matrix&,Real d=1.);       
  inline Sparsesym& init(const Indexmatrix&,Real d=1.);  
  inline Sparsesym& init(const Symmatrix&,Real d=1.); 
  inline Sparsesym& init(const Sparsemat&,Real d=1.);    
  inline Sparsesym& init(Integer nr);     
  Sparsesym& init(Integer nr,Integer nz,
                  const Integer *ini,const Integer *inj,const Real* va);
  Sparsesym& init(Integer nr,Integer nz,
                  const Indexmatrix& ini,const Indexmatrix& inj, const Matrix& va);
  
  //initialize to the same support as A but with constant value d; the same support will be generated even for d=0. 
  Sparsesym& init_support(const Sparsesym& A,Real d=0.);

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




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



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


  void dim(Integer& r, Integer& c) const {r=nr; c=nr;}

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

  Integer rowdim() const {return nr;}

  Integer coldim() const {return nr;}

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

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

    
  Real operator()(Integer i,Integer j) const;
  Real operator()(Integer i) const;    //vector of stacked columns

   Real operator[](Integer i) const;    //vector of stacked columns

  const Indexmatrix& get_colinfo() const {return colinfo;}

  const Indexmatrix& get_colindex() const {return colindex;}
  const Matrix& get_colval() const {return colval;}
  const Indexmatrix& get_suppind() const {return suppind;}
  const Indexmatrix& get_suppcol() const {return suppcol;}

  void get_edge_rep(Indexmatrix& I, Indexmatrix& J, Matrix& val) const;

  int contains_support(const Sparsesym& A) const;

  int check_support(Integer i,Integer j) const;
      




    

  friend Matrix diag(const Sparsesym& A);      //=(A(1,1),A(2,2),...)^t
  friend Sparsesym sparseDiag(const Matrix& A,Real tol);

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

  

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


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

  Sparsesym& support_xbpeya(const Sparsesym& y,Real alpha=1.,Real beta=0.);



  friend Sparsesym& xbpeya(Sparsesym& x,const Sparsesym& y,Real alpha,Real beta);
    
  friend Sparsesym& xeyapzb(Sparsesym& x,const Sparsesym& y,const Sparsesym& z,Real alpha,Real beta);
  
  friend Sparsesym& support_rankadd(const Matrix& A,Sparsesym& C,
                                    Real alpha,Real beta,int trans);

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

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

  

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


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

  Sparsesym& transpose(){return *this;}            

   

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


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

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

  inline Sparsesym& operator=(const Symmatrix &A);
  inline Symmatrix operator+(const Symmatrix &A) const;
  inline Symmatrix operator-(const Symmatrix &A) const;   

  inline Sparsesym& operator=(const Sparsemat &A);
  Sparsemat sparsemult(const Matrix& A) const;



  friend inline Symmatrix operator+(const Sparsesym &A,const Symmatrix &B);
  friend inline Symmatrix operator+(const Symmatrix &A,const Sparsesym &B);
  friend inline Symmatrix operator-(const Sparsesym& A,const Symmatrix &B);
  friend inline Symmatrix operator-(const Symmatrix &A,const Sparsesym &B);
  friend Real ip(const Symmatrix& A, const Sparsesym& B); //=trace(B^t*A)
  friend inline Real ip(const Sparsesym& A, const Symmatrix& B);

  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 inline Matrix operator*(const Sparsesym& A,const Sparsemat& B);
  friend inline Matrix operator*(const Sparsemat& A,const Sparsesym& B);


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


  friend Sparsesym abs(const Sparsesym& A);                 


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


  friend Real trace(const Sparsesym& A);               //=sum(diag(A))
  friend Real ip(const Sparsesym& A, const Sparsesym& B); //=trace(B^t*A)
  friend Real ip(const Matrix& A, const Sparsesym& B); //=trace(B^t*A)
  friend inline Real ip(const Sparsesym& A, const Matrix& B); 
  friend Real norm2(const Sparsesym& A); //=sqrt(ip(A,A))

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



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


   friend int equal(const Sparsesym& A, const Sparsesym& 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 Sparsesym &v);
  friend std::istream& operator>>(std::istream& i,Sparsesym &v);


};


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

  
  Matrix diag(const Sparsesym& A);      //=(A(1,1),A(2,2),...)^t

  Sparsesym sparseDiag(const Matrix& A,Real tol=SPARSE_ZERO_TOL);

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


  Sparsesym& xeyapzb(Sparsesym& x,const Sparsesym& y,
                     const Sparsesym& z,Real alpha=1.,Real beta=1.);

  Sparsesym& support_rankadd(const Matrix& A,Sparsesym& C,
                             Real alpha=1.,Real beta=0.,int trans=0);

  
  Real ip(const Symmatrix& A, const Sparsesym& B); //=trace(B^t*A)
  
  Sparsesym abs(const Sparsesym& A);
  
  Real trace(const Sparsesym& A);               //=sum(diag(A))
  
  Real ip(const Sparsesym& A, const Sparsesym& B); //=trace(B^t*A)
  
  Real ip(const Matrix& A, const Sparsesym& B); //=trace(B^t*A)
  
  Real norm2(const Sparsesym& A); //=sqrt(ip(A,A))
  
  Matrix sumrows(const Sparsesym& A);   //=(1 1 1 ... 1)*A
  
  Real sum(const Sparsesym& A);         //=(1 1 ... 1)*A*(1 1 ... 1)^t

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

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

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

inline void Sparsesym::init_to_zero()
{
 mtype=MTsparse;
 nr=0;
 tol=SPARSE_ZERO_TOL;
#if (CONICBUNDLE_DEBUG>=1)
 is_init=1;
#endif
}

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

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

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

inline Sparsesym& Sparsesym::init(Integer r)
{
nr=r;
colinfo.init(0,0,Integer(0)); colindex.init(0,0,Integer(0)); colval.init(0,0,0.);
chk_set_init(*this,1);
return *this;
}

inline Sparsesym::Sparsesym()
{ init_to_zero(); chk_set_init(*this,1);}
inline Sparsesym::Sparsesym(Integer r)
{ init_to_zero(); init(r);}
inline Sparsesym::Sparsesym(const Sparsesym& A,Real d):Memarrayuser()
{ init_to_zero(); xeya(A,d);}
inline Sparsesym::Sparsesym(const Matrix& A,Real d)
{ init_to_zero(); xeya(A,d);}
inline Sparsesym::Sparsesym(const Indexmatrix& A,Real d)
{ init_to_zero(); xeya(A,d);}

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

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

inline Real Sparsesym::operator()(Integer i) const
{ return (*this)(i%nr,i/nr); }
inline Real Sparsesym::operator[](Integer i) const
{ return (*this)(i%nr,i/nr); }


inline Sparsesym& xbpeya(Sparsesym& x,const Sparsesym& y,Real alpha=1.,Real beta=0.)
{
  if (beta==0.) return x.init(y,alpha); 
  Sparsesym B; xeyapzb(B,x,y,beta,alpha); swap(x,B); return x;
}

inline Sparsesym& Sparsesym::operator=(const Sparsesym &A)
{ return xeya(A);} 
inline Sparsesym& Sparsesym::operator+=(const Sparsesym &A)
{ Sparsesym B; xeyapzb(B,*this,A); swap(*this,B); return *this;}
inline Sparsesym& Sparsesym::operator-=(const Sparsesym &A)
{ Sparsesym B; xeyapzb(B,*this,A,1.,-1.); swap(*this,B); return *this;}
inline Sparsesym& Sparsesym::operator*=(Real d)
{chk_init(*this); colval*=d; return *this;}
inline Sparsesym& Sparsesym::operator/=(Real d)
{chk_init(*this);colval/=d;return *this;}
inline Sparsesym Sparsesym::operator-() const
{ return Sparsesym(*this,-1.); }

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


inline Sparsesym operator+(const Sparsesym& A,const Sparsesym& B) 
{Sparsesym C; return xeyapzb(C,A,B);}

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

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

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

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

inline Matrix operator*(const Sparsesym& A,const Matrix& B)
{ Matrix C; return genmult(A,B,C,1.,0.,0); }

inline Matrix operator*(const Matrix& A,const Sparsesym& B)
{ Matrix C; return genmult(A,B,C,1.,0.,0); }

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

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

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

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

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

inline Sparsesym transpose(const Sparsesym& A)
{return Sparsesym(A);}

inline Matrix sumcols(const Sparsesym& A)    //=A*(1 1 1 ... 1)^t
{Matrix s(sumrows(A)); return s.transpose();}


inline Matrix::Matrix(const Sparsesym& A,Real d)

{init_to_zero(); xeya(A,d); }

inline Matrix& Matrix::init(const Sparsesym& A,Real d)

{return xeya(A,d); }

inline Matrix& Matrix::operator=(const Sparsesym& A)

{return xeya(A); }

inline Matrix& Matrix::operator*=(const Sparsesym& A)

{ Matrix C; return *this=genmult(*this,A,C,1.,0.,0);}

inline Matrix& Matrix::operator+=(const Sparsesym& A)

{return xpeya(A); }

inline Matrix& Matrix::operator-=(const Sparsesym& A)

{return xpeya(A,-1.); }


inline Sparsesym::Sparsesym(const Symmatrix& A,Real d)
{ init_to_zero(); xeya(A,d);}
inline Sparsesym& Sparsesym::init(const Symmatrix& A,Real d)
{ return xeya(A,d);}
inline Sparsesym& Sparsesym::operator=(const Symmatrix& A)
{ return xeya(A);}
inline Symmatrix Sparsesym::operator+(const Symmatrix &A) const
{Symmatrix B(A); B+=*this; return B;}
inline Symmatrix Sparsesym::operator-(const Symmatrix &A) const
{Symmatrix B(A); B-=*this; return B;}

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

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

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

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

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

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

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

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

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

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

inline Symmatrix::Symmatrix(const Sparsesym& A,Real d)
    { init_to_zero(); xeya(A,d);}
inline Symmatrix& Symmatrix::init(const Sparsesym& A,Real d)          
    { return xeya(A,d);}
inline Symmatrix& Symmatrix::operator=(const Sparsesym& A)                
    { return xeya(A); }
inline Symmatrix& Symmatrix::operator+=(const Sparsesym& A)               
    { return xpeya(A);}
inline Symmatrix& Symmatrix::operator-=(const Sparsesym& A)           
    { return xpeya(A,-1.);}


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



}


#endif