indexmat.hxx

Go to the documentation of this file.

#ifndef CH_MATRIX_CLASSES__INDEXMAT_HXX
#define CH_MATRIX_CLASSES__INDEXMAT_HXX

#include <iostream>
#include <iomanip>
#include <vector>
#include <math.h>
#ifndef CH_MATRIX_CLASSES__MEMARRAY_HXX
#include "memarray.hxx"
#endif 
#ifndef CH_MATRIX_CLASSES__MYMATH_HXX
#include "mymath.hxx"
#endif


namespace CH_Matrix_Classes {


//everything involving a "Sparsesym" is implemented in "sparssym.cxx/hxx"
//everything else involving a "Sparsemat" is implemented in "sparsmat.cxx/hxx"
//everything else involving a "Symmatrix" is implemented in "symmat.cxx/hxx"
//everything else involving a "Matrix" is implemented in "matrix.cxx/hxx"



// **************************************************************************
//                               Range
// **************************************************************************

//specifies a range of numbers: {i: i=from+k*step for some k\in N_0 and i<=to}

class Range
{
public:
  Integer from;
  Integer to;
  Integer step;

  Range(Integer _from,Integer _to,Integer _step=1):from(_from),to(_to),step(_step){}
  ~Range(){} 
};




// **************************************************************************
//                              Indexmatrix
// **************************************************************************

class Matrix;
class Symmatrix;
class Sparsemat;
class Sparsesym;
class Indexmatrix;

// **************************************************************************
//                make non inline friends available outside
// **************************************************************************
// since some of them have default values, all of them are declared in advance


  Indexmatrix diag(const Indexmatrix& A);      //=(A(1,1),A(2,2),...)^t
  void swap(Indexmatrix &A,Indexmatrix &B);

  Indexmatrix& xbpeya(Indexmatrix& x,const Indexmatrix& y,Integer alpha=1,Integer beta=0,int ytrans=0);
  
  Indexmatrix& xeyapzb(Indexmatrix& x,const Indexmatrix& y,const Indexmatrix& z,Integer alpha=1,Integer beta=1);
  
  Indexmatrix& genmult(const Indexmatrix& A,const Indexmatrix& B,Indexmatrix& C,
                         Integer alpha=1,Integer beta=0,int atrans=0,int btrans=0);
  Indexmatrix transpose(const Indexmatrix& A);
  Indexmatrix abs(const Indexmatrix& A);                 
  Integer trace(const Indexmatrix& A);               //=sum(diag(A))
  Indexmatrix sumrows(const Indexmatrix& A);   
  Indexmatrix sumcols(const Indexmatrix& A);
  Integer sum(const Indexmatrix& A);     


  Indexmatrix operator<(const Indexmatrix &A,const Indexmatrix &B);
  Indexmatrix operator<=(const Indexmatrix &A,const Indexmatrix &B);
  Indexmatrix operator==(const Indexmatrix &A,const Indexmatrix &B);
  Indexmatrix operator!=(const Indexmatrix &A,const Indexmatrix &B);
  Indexmatrix operator<(const Indexmatrix &A,Integer d);
  Indexmatrix operator>(const Indexmatrix &A,Integer d);
  Indexmatrix operator<=(const Indexmatrix &A,Integer d);
  Indexmatrix operator>=(const Indexmatrix &A,Integer d);
  Indexmatrix operator==(const Indexmatrix &A,Integer d);
  Indexmatrix operator!=(const Indexmatrix &A,Integer d);
   
  Indexmatrix minrows(const Indexmatrix& A);   
  Indexmatrix mincols(const Indexmatrix& A);   
  Integer min(const Indexmatrix& A,Integer *iindex=0,Integer *jindex=0);
  Indexmatrix maxrows(const Indexmatrix& A);   
  Indexmatrix maxcols(const Indexmatrix& A);   
  Integer max(const Indexmatrix& A,Integer *iindex=0,Integer *jindex=0);
  
  inline Indexmatrix sortindex(const Indexmatrix& vec,bool nondecreasing=true);
  void sortindex(const Indexmatrix& vec,Indexmatrix &ind, bool nondecreasing=true);

  std::ostream& operator<<(std::ostream& o,const Indexmatrix &A);
  std::istream& operator>>(std::istream& i,Indexmatrix &A);
  



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

        
private:
  static const Mtype mtype;    
    Integer mem_dim;   
    Integer nr,        
            nc;        
    Integer *m;        

    bool is_init;      

    inline void init_to_zero();   
    
public:

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


  inline Indexmatrix();                  

  inline Indexmatrix(const Indexmatrix& A,Integer d=1);

  inline Indexmatrix(const Range&);           

  inline Indexmatrix(Integer nr,Integer nc);         

  inline Indexmatrix(Integer nr,Integer nc,Integer d);  

  inline Indexmatrix(Integer nr,Integer nc,const Integer* dp,Integer incr=1); 

  inline Indexmatrix(const std::vector<Integer>& vec);

  ~Indexmatrix(){memarray->free(m);}
  
#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 Indexmatrix& init(const Indexmatrix &A,Integer d=1);    

  Indexmatrix& init(const Range&);

  inline Indexmatrix& init(Integer nr,Integer nc,Integer d);                    

  inline Indexmatrix& init(Integer nr,Integer nc,const Integer *dp,Integer incr=1);

  inline Indexmatrix& init(const std::vector<Integer>& vec);
    
  void newsize(Integer nr,Integer nc);  //resize matrix without initialization



  Indexmatrix(const Matrix&);          
  Indexmatrix(const Symmatrix&);      
  Indexmatrix(const Sparsemat&);      
  Indexmatrix(const Sparsesym&);      


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


  void dim(Integer& _nr, Integer& _nc) const {_nr=nr; _nc=nc;}

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

  Integer rowdim() const {return nr;}

  Integer coldim() const {return nc;}

  Mtype get_mtype() const {return mtype;}


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


  
  inline Integer& operator()(Integer i,Integer j);      

  inline Integer& operator()(Integer i);      //index to vector of stacked columns

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

  inline Integer operator()(Integer i) const; //index to vector of stacked columns

  Indexmatrix operator()(const Indexmatrix& vecrow,const Indexmatrix& veccol) const;

    Indexmatrix operator()(const Indexmatrix& A) const;
      //*this matrix is interpreted as vector, the resulting matrix has
      //has the same shape as A
  
  inline Integer& operator[](Integer i);      //{return (*this)(i);}

  inline Integer operator[](Integer i) const; //{return (*this)(i);}
  
  Indexmatrix col(Integer i) const;          //returns column i as column vector
  Indexmatrix row(Integer i) const;          //returns row i as row vector
  Indexmatrix cols(const Indexmatrix &vec) const;  //returns cols as indexed by vec
  
  Indexmatrix rows(const Indexmatrix &vec) const;  //returns rows as indexed by vec


  Indexmatrix& triu(Integer d=0);    // (*this)(i,j)=0 for i<j+d

  Indexmatrix& tril(Integer d=0);    // (*this)(i,j)=0 for i>j+d

  Indexmatrix& subassign(const Indexmatrix& vecrow,const Indexmatrix& veccol,
                    const Indexmatrix& A);
     //(*this)(vecrow(i),veccol(j))=A(i,j) for all i,j

  Indexmatrix& subassign(const Indexmatrix& vec,const Indexmatrix& A);

  Indexmatrix& delete_rows(const Indexmatrix& ind,bool sorted_increasingly=false);
  Indexmatrix& delete_cols(const Indexmatrix& ind,bool sorted_increasingly=false);

  Indexmatrix& insert_row(Integer i,const Indexmatrix& v); // 0<=i<=nr, v vector
  Indexmatrix& insert_col(Integer i,const Indexmatrix& v); // 0<=i<=nc, v vector

  inline Indexmatrix& reduce_length(Integer n);  

  Indexmatrix& concat_right(const Indexmatrix& A);
  Indexmatrix& concat_below(const Indexmatrix& A);
  Indexmatrix& concat_below(Integer d);   //only for column vectors (nr==1)
  Indexmatrix& concat_right(Integer d);   //only for row vectors (nc==1)
  Indexmatrix& enlarge_right(Integer addnc); 
  Indexmatrix& enlarge_below(Integer addnr);   
  Indexmatrix& enlarge_right(Integer addnc,Integer d); 
  Indexmatrix& enlarge_below(Integer addnr,Integer d);   
  Indexmatrix& enlarge_right(Integer addnc,const Integer* dp,Integer d=1); 
  Indexmatrix& enlarge_below(Integer addnr,const Integer *dp,Integer d=1);   


  Integer* get_store() {return m;}       //use cautiously, do not use delete!
  const Integer* get_store() const {return m;}   //use cautiously





  friend Indexmatrix diag(const Indexmatrix& A);      //=(A(1,1),A(2,2),...)^t

  friend Indexmatrix triu(const Indexmatrix& A,Integer d);

  friend Indexmatrix tril(const Indexmatrix& A,Integer d);

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

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


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


  Indexmatrix& xeya(const Indexmatrix& A,Integer d=1); 
  Indexmatrix& xpeya(const Indexmatrix& A,Integer d=1);
    


  friend Indexmatrix& xbpeya(Indexmatrix& x,const Indexmatrix& y,Integer alpha,Integer beta,int ytrans);
  
  friend Indexmatrix& xeyapzb(Indexmatrix& x,const Indexmatrix& y,const Indexmatrix& z,Integer alpha,Integer beta);
  
  friend Indexmatrix& genmult(const Indexmatrix& A,const Indexmatrix& B,Indexmatrix& C,
                         Integer alpha,Integer beta,int atrans,int btrans);
    
    

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


  inline Indexmatrix& operator=(const Indexmatrix &A);
  inline Indexmatrix& operator*=(const Indexmatrix &s);
  inline Indexmatrix& operator+=(const Indexmatrix &v);
  inline Indexmatrix& operator-=(const Indexmatrix &v);
  inline Indexmatrix& operator%=(const Indexmatrix &A);   //Hadamard product!!
  inline Indexmatrix operator-() const;
  
  inline Indexmatrix& operator*=(Integer d);
  inline Indexmatrix& operator/=(Integer d);
  inline Indexmatrix& operator%=(Integer d); 
  inline Indexmatrix& operator+=(Integer d);
  inline Indexmatrix& operator-=(Integer d);
  
  Indexmatrix& transpose();
  


  friend inline Indexmatrix operator*(const Indexmatrix &A,const Indexmatrix& B);
  friend inline Indexmatrix operator+(const Indexmatrix &A,const Indexmatrix& B);
  friend inline Indexmatrix operator-(const Indexmatrix &A,const Indexmatrix& B);
  friend inline Indexmatrix operator%(const Indexmatrix &A,const Indexmatrix& B);
  friend inline Indexmatrix operator*(const Indexmatrix &A,Integer d);
  friend inline Indexmatrix operator*(Integer d,const Indexmatrix &A);
  friend inline Indexmatrix operator/(const Indexmatrix& A,Integer d);
  friend inline Indexmatrix operator%(const Indexmatrix& A,Integer d);
  
  friend inline Indexmatrix operator+(const Indexmatrix& A,Integer d);
  friend inline Indexmatrix operator+(Integer d,const Indexmatrix& A);
  friend inline Indexmatrix operator-(const Indexmatrix& A,Integer d);
  friend inline Indexmatrix operator-(Integer d,const Indexmatrix& A);

  friend Indexmatrix transpose(const Indexmatrix& A);


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


  friend inline std::vector<int>& assign(std::vector<int>& vec,
                                            const Indexmatrix& A); 
  friend inline std::vector<long>& assign(std::vector<long>& vec,
                                             const Indexmatrix& A); 



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


  Indexmatrix& rand(Integer nr,Integer nc,Integer lowerb,Integer upperb,CH_Tools::GB_rand* random_generator=0);
  Indexmatrix& shuffle(CH_Tools::GB_rand* random_generator=0);
  Indexmatrix& sign(void);
  Indexmatrix& abs(void);        



  friend Indexmatrix rand(Integer nr,Integer nc,Integer lb,Integer ub,CH_Tools::GB_rand* random_generator);
  friend inline Indexmatrix sign(const Indexmatrix& A);
  friend Indexmatrix abs(const Indexmatrix& A);                 


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


  Indexmatrix& scale_rows(const Indexmatrix& vec);    
  Indexmatrix& scale_cols(const Indexmatrix& vec); 



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

  
  friend Indexmatrix sumrows(const Indexmatrix& A);   
  friend Indexmatrix sumcols(const Indexmatrix& A);
  friend Integer sum(const Indexmatrix& A);     



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


  Indexmatrix find() const;   //finds nonzeros
  Indexmatrix find_number(Integer num=0) const; 



  friend Indexmatrix operator<(const Indexmatrix &A,const Indexmatrix &B);
  friend inline Indexmatrix operator>(const Indexmatrix &A,const Indexmatrix &B);
  friend Indexmatrix operator<=(const Indexmatrix &A,const Indexmatrix &B);
  friend inline Indexmatrix operator>=(const Indexmatrix &A,const Indexmatrix &B);
  friend Indexmatrix operator==(const Indexmatrix &A,const Indexmatrix &B);
  friend Indexmatrix operator!=(const Indexmatrix &A,const Indexmatrix &B);
  friend Indexmatrix operator<(const Indexmatrix &A,Integer d);
  friend Indexmatrix operator>(const Indexmatrix &A,Integer d);
  friend Indexmatrix operator<=(const Indexmatrix &A,Integer d);
  friend Indexmatrix operator>=(const Indexmatrix &A,Integer d);
  friend Indexmatrix operator==(const Indexmatrix &A,Integer d);
  friend Indexmatrix operator!=(const Indexmatrix &A,Integer d);
  friend inline Indexmatrix operator<(Integer d,const Indexmatrix &A);
  friend inline Indexmatrix operator>(Integer d,const Indexmatrix &A);
  friend inline Indexmatrix operator<=(Integer d,const Indexmatrix &A);
  friend inline Indexmatrix operator>=(Integer d,const Indexmatrix &A);
  friend inline Indexmatrix operator==(Integer d,const Indexmatrix &A);
  friend inline Indexmatrix operator!=(Integer d,const Indexmatrix &A);

  friend inline bool equal(const Indexmatrix& A,const Indexmatrix& b);

  friend Indexmatrix minrows(const Indexmatrix& A);   
  friend Indexmatrix mincols(const Indexmatrix& A);   
  friend Integer min(const Indexmatrix& A,Integer *iindex,Integer *jindex);
  friend Indexmatrix maxrows(const Indexmatrix& A);   
  friend Indexmatrix maxcols(const Indexmatrix& A);   
  friend Integer max(const Indexmatrix& A,Integer *iindex,Integer *jindex);
  
  friend inline Indexmatrix sortindex(const Indexmatrix& vec,bool nondecreasing);
  friend void sortindex(const Indexmatrix& vec,Indexmatrix &ind, bool nondecreasing);
  
  friend inline Indexmatrix find(const Indexmatrix& A);
  friend Indexmatrix find_number(const Indexmatrix& A,Integer num);


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


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

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




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

};




//@cond

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

inline void Indexmatrix::init_to_zero()
{
 nr=nc=0;mem_dim=0;m=0;
 chk_set_init(*this,1);
}

inline Indexmatrix& Indexmatrix::init(Integer inr,Integer inc,Integer d)
{
 newsize(inr,inc);
 mat_xea(nr*nc,m,d);
 chk_set_init(*this,1);
 return *this;
}

inline Indexmatrix& Indexmatrix::init(Integer inr,Integer inc,const Integer* d,Integer incr)
{
 newsize(inr,inc);
 if (incr==1) mat_xey(nr*nc,m,d);
 else mat_xey(nr*nc,m,Integer(1),d,incr);
 chk_set_init(*this,1);
 return *this;
}

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

inline Indexmatrix& Indexmatrix::init(const std::vector<Integer>& vec)
{
  newsize(Integer(vec.size()),1); chk_set_init(*this,1);
  for(Integer i=0;i<nr;i++) m[i]=vec[(unsigned long)(i)];
  return *this;
}

inline Indexmatrix::Indexmatrix()
{
 init_to_zero();
}

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

inline Indexmatrix::Indexmatrix(Integer inr,Integer inc)
{
 init_to_zero();
 newsize(inr,inc);
}

inline Indexmatrix::Indexmatrix(const Range& range)
{
 init_to_zero();
 init(range);
}

inline Indexmatrix::Indexmatrix(Integer inr,Integer inc,Integer d)
{
 init_to_zero();
 init(inr,inc,d);
}

inline Indexmatrix::Indexmatrix(Integer inr,Integer inc,const Integer *d,Integer incr)
{
 init_to_zero();
 init(inr,inc,d,incr);
}

inline Indexmatrix::Indexmatrix(const std::vector<Integer>& vec)
{
 init_to_zero();
 init(vec);
}

inline Integer& Indexmatrix::operator()(Integer i,Integer j)
{
 chk_range(i,j,nr,nc);
 return m[j*nr+i];
}

inline Integer& Indexmatrix::operator()(Integer i) 
{
 chk_range(i,0,nr*nc,1);
 return m[i];
}

inline Integer Indexmatrix::operator()(Integer i,Integer j) const
{
 chk_range(i,j,nr,nc);
 return m[j*nr+i];
}

inline Integer Indexmatrix::operator()(Integer i) const
{
 chk_range(i,0,nr*nc,1);
 return m[i];
}

inline Integer& Indexmatrix::operator[](Integer i)
{return (*this)(i);}

inline Integer Indexmatrix::operator[](Integer i) const 
{return (*this)(i);}

inline Indexmatrix& Indexmatrix::reduce_length(Integer n) 
{ nr=min(nr*nc,max(Integer(0),n)); nc=1; return *this;}


inline Integer ip(const Indexmatrix& A, const Indexmatrix& B)
{
 chk_add(A,B);
 return mat_ip(A.nc*A.nr,A.m,B.m);
}

inline Real norm2(const Indexmatrix& A)
{
 chk_init(A);
 return ::sqrt(double(mat_ip(A.nc*A.nr,A.m,A.m)));
}

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

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

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

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

inline Indexmatrix& Indexmatrix::operator%=(const Indexmatrix &A)
{ chk_add(*this,A); mat_xhadey(nr*nc,m,A.m); return *this; }

inline Indexmatrix Indexmatrix::operator-() const
{ return Indexmatrix(*this,-1); }

inline Indexmatrix& Indexmatrix::operator*=(Integer d)
{ chk_init(*this); mat_xmultea(nr*nc,m,d); return *this; }

inline Indexmatrix& Indexmatrix::operator/=(Integer d)
{ chk_init(*this); mat_xdivea(nr*nc,m,d); return *this; }

inline Indexmatrix& Indexmatrix::operator%=(Integer d)
{ chk_init(*this); mat_xmodea(nr*nc,m,d); return *this; }

inline Indexmatrix& Indexmatrix::operator+=(Integer d)
{ chk_init(*this); mat_xpea(nr*nc,m,d); return *this; }

inline Indexmatrix& Indexmatrix::operator-=(Integer d)
{ chk_init(*this); mat_xpea(nr*nc,m,-d); return *this; }


inline Indexmatrix operator*(const Indexmatrix &A,const Indexmatrix &B) 
    {Indexmatrix C; return genmult(A,B,C);}
inline Indexmatrix operator+(const Indexmatrix &A,const Indexmatrix &B)
    {Indexmatrix C; return xeyapzb(C,A,B,1,1);}
inline Indexmatrix operator-(const Indexmatrix &A,const Indexmatrix &B)
    {Indexmatrix C; return xeyapzb(C,A,B,1,-1);}
inline Indexmatrix operator%(const Indexmatrix &A,const Indexmatrix &B) 
    {Indexmatrix C(A); return C%=B;}
inline Indexmatrix operator*(const Indexmatrix &A,Integer d)          
    {return Indexmatrix(A,d);}
inline Indexmatrix operator*(Integer d,const Indexmatrix &A)
    {return Indexmatrix(A,d);}
inline Indexmatrix operator/(const Indexmatrix &A,Integer d)
    {Indexmatrix B(A); return B/=d;}
inline Indexmatrix operator%(const Indexmatrix &A,Integer d)
    {Indexmatrix B(A); return B%=d;}
inline Indexmatrix operator+(const Indexmatrix &A,Integer d)           
    {Indexmatrix B(A); return B+=d;}
inline Indexmatrix operator+(Integer d,const Indexmatrix &A)           
    {Indexmatrix B(A); return B+=d;}
inline Indexmatrix operator-(const Indexmatrix &A,Integer d)       
    {Indexmatrix B(A); return B-=d;}
inline Indexmatrix operator-(Integer d,const Indexmatrix &A)       
    {Indexmatrix B(A,-1); return B+=d;}

inline Indexmatrix triu(const Indexmatrix& A,Integer i=0)
          {Indexmatrix B(A); B.triu(i); return B;}
inline Indexmatrix tril(const Indexmatrix& A,Integer i=0)
          {Indexmatrix B(A); B.tril(i); return B;}

inline Indexmatrix concat_right(const Indexmatrix& A,const Indexmatrix& B)
          {Indexmatrix C(A.dim()+B.dim(),1);C=A;C.concat_right(B);return C;}
inline Indexmatrix concat_below(const Indexmatrix& A,const Indexmatrix& B)
          {Indexmatrix C(A.dim()+B.dim(),1);C=A;C.concat_below(B);return C;}
    
inline void swap(Indexmatrix &A,Indexmatrix &B)
{ 
 Integer *hm=A.m;A.m=B.m;B.m=hm;
 Integer hi=A.nr;A.nr=B.nr;B.nr=hi;
 hi=A.nc;A.nc=B.nc;B.nc=hi;
 hi=A.mem_dim;A.mem_dim=B.mem_dim;B.mem_dim=hi;
#if (CONICBUNDLE_DEBUG>=1)
 bool hb=A.is_init;A.is_init=B.is_init;B.is_init=hb;
#endif
}

inline Indexmatrix rand(Integer rows,Integer cols,Integer lb,Integer ub,CH_Tools::GB_rand* random_generator=0)
{Indexmatrix A; return A.rand(rows,cols,lb,ub,random_generator);}
inline Indexmatrix sign(const Indexmatrix& A)
          {Indexmatrix B(A); return B.sign();}

inline Indexmatrix operator>(const Indexmatrix &A,const Indexmatrix &B)
{return B<A;}
inline Indexmatrix operator>=(const Indexmatrix &A,const Indexmatrix &B)
{return B<=A;}
inline Indexmatrix operator<(Integer d,const Indexmatrix &A)
{return A>d;}
inline Indexmatrix operator>(Integer d,const Indexmatrix &A)
{return A<d;}
inline Indexmatrix operator<=(Integer d,const Indexmatrix &A)
{return A>=d;}
inline Indexmatrix operator>=(Integer d,const Indexmatrix &A)
{return A<=d;}
inline Indexmatrix operator==(Integer d,const Indexmatrix &A)
{return A==d;}
inline Indexmatrix operator!=(Integer d,const Indexmatrix &A)
{return A!=d;}

inline bool equal(const Indexmatrix& A,const Indexmatrix& B)
{return ((A.nc==B.nc)&&(A.nr==B.nr)&&(mat_equal(A.nr*A.nc,A.get_store(),B.get_store())));}

inline Indexmatrix sortindex(const Indexmatrix& vec,bool nondecreasing)
{Indexmatrix ind;sortindex(vec,ind,nondecreasing);return ind;}

inline Indexmatrix find(const Indexmatrix& A)
          {return A.find();}
inline Indexmatrix find_number(const Indexmatrix& A,Integer num=0)
          {return A.find_number(num);}

inline std::vector<int>& assign(std::vector<int>& vec,const Indexmatrix& A)
{
  chk_init(A);
  vec.resize((unsigned long)(A.dim()));
  for(Integer i=0;i<A.dim();i++) vec[(unsigned long)(i)]=int(A(i));
  return vec;
}

inline std::vector<long>& assign(std::vector<long>& vec,const Indexmatrix& A)
{
  chk_init(A);
  vec.resize((unsigned long)(A.dim()));
  for(Integer i=0;i<A.dim();i++) vec[(unsigned long)(i)]=long(A(i));
  return vec;
}

  //@endcond

}


#ifndef CH_MATRIX_CLASSES__MATRIX_HXX
#include "matrix.hxx"
#endif


#endif