matrix.hxx

Go to the documentation of this file.

#ifndef CH_MATRIX_CLASSES__MATRIX_HXX
#define CH_MATRIX_CLASSES__MATRIX_HXX

#include <math.h>
#ifndef CH_MATRIX_CLASSES__INDEXMAT_HXX
#include "indexmat.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"

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


class Realrange
{
public:
  Real from;
  Real to;
  Real step;
  Real tol;
  Realrange(Real _from,Real _to,Real _step=1,Real _tol=1e-8):
    from(_from),to(_to),step(_step),tol(_tol){}
  ~Realrange(){} 
};


// **************************************************************************
//                               Matrix
// **************************************************************************


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

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

  //--- auxilliary storage space for repeated use of information, e.g., in qr
  enum AuxTask {
    none, 
    qr 
  };
  AuxTask aux_task;  
  Integer aux_dim;   
  Real *aux_m;      
  
  bool is_init;      
  
  inline void init_to_zero();   

public:

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



  inline Matrix();                              
  inline Matrix(const Matrix&,Real d=1.,int atrans=0);       
  inline Matrix(const Realrange&);              
  inline Matrix(Integer nr,Integer nc);         
  inline Matrix(Integer nr,Integer nc,Real d);  
  inline Matrix(Integer nr,Integer nc,const Real* dp,Integer incr=1,Real d=1.); 
  inline Matrix(const std::vector<Real>& vec);
  ~Matrix(){memarray->free(m);if (aux_m) memarray->free(aux_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 Matrix& init(const Matrix &A,Real d=1.,int atrans=0);                   
  inline Matrix& init(const Indexmatrix& A,Real d=1.);
  inline Matrix& init(const Sparsemat& A, Real d=1.); //*this=A*d;
  inline Matrix& init(const Symmatrix& S,Real d=1.);
  inline Matrix& init(const Sparsesym&,Real d=1.);
  Matrix& init(const Realrange&);
  inline Matrix& init(Integer nr,Integer nc,Real d);                    
  inline Matrix& init(Integer nr,Integer nc,const Real *dp,Integer incr=1,Real d=1.); 
  inline Matrix& init(const std::vector<Real>& vec);
  inline Matrix& init_diag(int nr,Real d=1.);
  inline Matrix& init_diag(const Matrix& vec,Real d=1.);
  inline Matrix& init_diag(const Indexmatrix& vec,Real d=1.);
    
  void newsize(Integer nr,Integer nc);  //resize matrix without initialization




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


  //----------------------------------------------------
  //----  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 Real& operator()(Integer i,Integer j);      

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

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

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

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

  Matrix operator()(const Indexmatrix& A) const;
 
  
  inline Real& operator[](Integer i);      //{return (*this)(i);}
  inline Real operator[](Integer i) const; //{return (*this)(i);}
  
  Matrix col(Integer i) const;          //returns column i as column vector
  Matrix row(Integer i) const;          //returns row i as row vector
  Matrix cols(const Indexmatrix &vec) const;  //returns cols as indexed by vec
  Matrix rows(const Indexmatrix &vec) const;  //returns rows as indexed by vec
  Matrix& swap_rowsij(Integer i,Integer j);
  Matrix& swap_colsij(Integer i,Integer j);
  Matrix& pivot_permute_rows(const Indexmatrix & piv,bool inverse=false);
  Matrix& pivot_permute_cols(const Indexmatrix & piv,bool inverse=false);

  Matrix& triu(Integer d=0);    // (*this)(i,j)=0 for j<i+d
  Matrix& tril(Integer d=0);    // (*this)(i,j)=0 for j>i+d
  
  Matrix& subassign(const Indexmatrix& vecrow,const Indexmatrix& veccol,
                    const Matrix& A);
     //(*this)(vecrow(i),veccol(j))=A(i,j) for all i,j
  Matrix& subassign(const Indexmatrix& vec,const Matrix& A);
     //(*this)(vec(i))=A(i);
     //*this, vec, and A may be rect matrices but will be used as vectors
  
  Matrix& delete_rows(const Indexmatrix& ind,bool sorted_increasingly=false);
  Matrix& delete_cols(const Indexmatrix& ind,bool sorted_increasingly=false);

  Matrix& insert_row(Integer i,const Matrix& v); // 0<=i<=nr, v vector
  Matrix& insert_col(Integer i,const Matrix& v); // 0<=i<=nc, v vector
  
  inline Matrix& reduce_length(Integer n);  
     //interpret *this as a column vector and cut at n 

  Matrix& concat_right(const Matrix& A,int Atrans=0);
  Matrix& concat_below(const Matrix& A);
  Matrix& concat_right(Real d);   //only for column vectors (nr==1)
  Matrix& concat_below(Real d);   //only for row vectors (nc==1)
  Matrix& enlarge_right(Integer addnc); 
  Matrix& enlarge_below(Integer addnr);   
  Matrix& enlarge_right(Integer addnc,Real d); 
  Matrix& enlarge_below(Integer addnr,Real d);   
  Matrix& enlarge_right(Integer addnc,const Real* dp,Real d=1.); 
  Matrix& enlarge_below(Integer addnr,const Real *dp,Real d=1.);   
  
  Real* get_store() {return m;}       //use cautiously, do not use delete!
  const Real* get_store() const {return m;}   //use cautiously





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

  friend Matrix triu(const Matrix& A,Integer i);
  friend Matrix tril(const Matrix& A,Integer i);


  friend Matrix concat_right(const Matrix& A,const Matrix& B);
  friend Matrix concat_below(const Matrix& A,const Matrix& B);

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


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


  Matrix& xeya(const Matrix& A,Real d=1.,int atrans=0);   //*this=d*A
  Matrix& xpeya(const Matrix& A,Real d=1.);  //*this+=d*A;
  Matrix& xeya(const Indexmatrix& A,Real d=1.);   //*this=d*A
  Matrix& xpeya(const Indexmatrix& A,Real d=1.);  //*this+=d*A;
    


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

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


  inline Matrix& operator=(const Matrix &A);
  inline Matrix& operator*=(const Matrix &s);
  inline Matrix& operator+=(const Matrix &v);
  inline Matrix& operator-=(const Matrix &v);
  inline Matrix& operator%=(const Matrix &A);  
  inline Matrix& operator/=(const Matrix &A);  
  inline Matrix operator-() const;
  
  inline Matrix& operator*=(Real d);
  inline Matrix& operator/=(Real d);
  inline Matrix& operator+=(Real d);
  inline Matrix& operator-=(Real d);
  
  Matrix& transpose();            //transposes itself
  


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

  friend Matrix transpose(const Matrix& A);


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


  
  Matrix& xeya(const Symmatrix& A,Real d=1.);  //*this=A*d;
  Matrix& xpeya(const Symmatrix& A,Real d=1.); //*this+=A*d;
  inline Matrix& operator=(const Symmatrix& S);
  inline Matrix& operator*=(const Symmatrix& S);
  inline Matrix& operator+=(const Symmatrix& S);
  inline Matrix& operator-=(const Symmatrix& S);

  Matrix& xeya(const Sparsesym& A,Real d=1.);  //*this=A*d;
  Matrix& xpeya(const Sparsesym& A,Real d=1.); //*this+=A*d;
  inline Matrix& operator=(const Sparsesym &);
  inline Matrix& operator*=(const Sparsesym& S);
  inline Matrix& operator+=(const Sparsesym& S);
  inline Matrix& operator-=(const Sparsesym& S);
    
  Matrix& xeya(const Sparsemat& A,Real d=1.);  //*this=A*d;
  Matrix& xpeya(const Sparsemat& A,Real d=1.); //*this+=A*d;
  inline Matrix& operator=(const Sparsemat & A);      
  inline Matrix& operator*=(const Sparsemat &A);      
  inline Matrix& operator+=(const Sparsemat &A);      
  inline Matrix& operator-=(const Sparsemat &A);      



  friend std::vector<double>& assign(std::vector<double>& vec,const Matrix& A);


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

  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);

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


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


  Matrix& rand(Integer nr,Integer nc,CH_Tools::GB_rand* random_generator=0);
  Matrix& rand_normal(Integer nr,Integer nc,Real mean=0.,Real variance=1.,int generator_type=0);
  Matrix& shuffle(CH_Tools::GB_rand* random_generator=0);
  Matrix& inv(void);                        //reciprocal value componentwise
  Matrix& sqrt(void);
  Matrix& sqr(void);
  Matrix& sign(Real tol=1e-12);
  Matrix& floor(void);
  Matrix& ceil(void);
  Matrix& rint(void);
  Matrix& round(void);
  Matrix& abs(void);        



  friend Matrix rand(Integer nr,Integer nc,CH_Tools::GB_rand* random_generator);
  friend Matrix inv(const Matrix& A);
  friend Matrix sqrt(const Matrix& A);
  friend Matrix sqr(const Matrix& A);
  friend Matrix sign(const Matrix& A,Real tol);
  friend Matrix floor(const Matrix& A);
  friend Matrix ceil(const Matrix& A);
  friend Matrix rint(const Matrix& A);
  friend Matrix round(const Matrix& A);
  friend Matrix abs(const Matrix& A);                 


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

  bool contains_nan();

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

  //----- triangular routines
  int triu_solve(Matrix& rhs,Real tol=1e-10);   //only elements i>=j are used
  int tril_solve(Matrix& rhs,Real tol=1e-10);   //only elements i<=j are used


  //----- QR-Factorization
  int QR_factor(Real tol=1e-10);                            //factorization stored in *this
  int QR_factor(Matrix& Q,Real tol=1e-10);                   //afterwards *this is R
  int QR_factor(Matrix& Q,Matrix& R,Real tol) const;    //*this is unchanged

  int QR_factor(Indexmatrix& piv,Real tol=1e-10);
  int QR_factor_relpiv(Indexmatrix& piv,Real tol=1e-10);
  int QR_factor(Matrix& Q,Indexmatrix& piv,Real tol=1e-10);
  int QR_factor(Matrix& Q,Matrix& R,Indexmatrix& piv,Real tol) const;
  int Qt_times(Matrix& A,Integer r) const;
  int Q_times(Matrix& A,Integer r) const;
  int times_Q(Matrix& A,Integer r) const;
  int QR_solve(Matrix& rhs,Real tol=1e-10);  
  
  int QR_concat_right(const Matrix& A,Indexmatrix& piv,Integer r,Real tol=1e-10);


  //----- Least squares  
  int ls(Matrix & rhs, Real tol);
  
  int nnls(Matrix& rhs,Matrix* dual=0,Real tol=1e-10) const;



  friend Real trace(const Matrix& A);               //=sum(diag(A))
  friend Real ip(const Matrix& A, const Matrix& B); //=trace(B^t*A)
  friend Real ip_min_max(const Matrix& A, const Matrix& B,Real &minval,Real& maxval); //=trace(B^t*A), minval=min(B_ij*A_ij), maxval=max(B_ij*A_ij)
  friend Real colip(const Matrix& A,Integer j,const Matrix* scaling);
  friend Real rowip(const Matrix& A,Integer i,const Matrix* scaling);
  friend Matrix colsip(const Matrix& A);
  friend Matrix rowsip(const Matrix& A);
  friend Real norm2(const Matrix& A);               //=sqrt(ip(A,A));
  friend Real normDsquared(const Matrix& A,const Matrix& d,int atrans,int dinv);       

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

  //----- Householder rotations
  friend Matrix house(const Matrix &A,Integer i,Integer j,Real tol);
  friend int rowhouse(Matrix &A,const Matrix& v,Integer i,Integer j);
  friend int colhouse(Matrix &A,const Matrix& v,Integer i,Integer j);   

  friend int QR_factor(const Matrix& A,Matrix& Q,Matrix &R,Real tol);
  friend int QR_factor(const Matrix& A,Matrix& Q,Matrix &R,Indexmatrix& piv,Real tol);


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


  Indexmatrix find(Real tol=1e-10) const;   //finds nonzeros
  Indexmatrix find_number(Real num=0.,Real tol=1e-10) const; 



  friend Matrix operator<(const Matrix &A,const Matrix &B);
  friend Matrix operator>(const Matrix &A,const Matrix &B);
  friend Matrix operator<=(const Matrix &A,const Matrix &B);
  friend Matrix operator>=(const Matrix &A,const Matrix &B);
  friend Matrix operator==(const Matrix &A,const Matrix &B);
  friend Matrix operator!=(const Matrix &A,const Matrix &B);
    friend Matrix operator<(const Matrix &A,Real d);
  friend Matrix operator>(const Matrix &A,Real d);
  friend Matrix operator<=(const Matrix &A,Real d);
  friend Matrix operator>=(const Matrix &A,Real d);
  friend Matrix operator==(const Matrix &A,Real d);
  friend Matrix operator!=(const Matrix &A,Real d);
  friend Matrix operator<(Real d,const Matrix &A);
  friend Matrix operator>(Real d,const Matrix &A);
  friend Matrix operator<=(Real d,const Matrix &A);
  friend Matrix operator>=(Real d,const Matrix &A);
  friend Matrix operator==(Real d,const Matrix &A);
  friend Matrix operator!=(Real d,const Matrix &A);

  friend bool equal(const Matrix& A,const Matrix& B);

 friend Matrix minrows(const Matrix& A);   //min of each column (over the rows)
  friend Matrix mincols(const Matrix& A);   //min of each row (over the columns)
  friend Real min(const Matrix& A,Integer *iindex,Integer *jindex);
  friend Matrix maxrows(const Matrix& A);   //similar
  friend Matrix maxcols(const Matrix& A);   //similar
  friend Real max(const Matrix& A,Integer *iindex,Integer *jindex);
  
  friend Indexmatrix sortindex(const Matrix& vec,bool nondecreasing);
  friend void sortindex(const Matrix& vec,Indexmatrix &ind,bool nondecreasing);
  
  friend Indexmatrix find(const Matrix& A,Real tol);
  friend Indexmatrix find_number(const Matrix& A,Real num,Real tol);


  //--------------------------------
  //----  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 Matrix &v);
  friend std::istream& operator>>(std::istream& i,Matrix &v);
  
  
};


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

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

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

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

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

Matrix transpose(const Matrix& A);

std::vector<double>& assign(std::vector<double>& vec,const Matrix& A);

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

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

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

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

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

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

Matrix abs(const Matrix& A);

Real trace(const Matrix& A);               

Real normDsquared(const Matrix& A,const Matrix& d,int atrans=0,int dinv=0);      

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

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

Real sum(const Matrix& A);         //=(1 1 ... 1)*A*(1 1 ... 1)^t

Matrix house(const Matrix &x,Integer i=0,Integer j=0,Real tol=1e-10);

int rowhouse(Matrix &A,const Matrix& v,Integer i=0,Integer j=0);

int colhouse(Matrix &A,const Matrix& v,Integer i=0,Integer j=0);

Matrix operator<(const Matrix &A,const Matrix &B);

Matrix operator<=(const Matrix &A,const Matrix &B);

Matrix operator==(const Matrix &A,const Matrix &B);

Matrix operator!=(const Matrix &A,const Matrix &B);

Matrix operator<(const Matrix &A,Real d);

Matrix operator>(const Matrix &A,Real d);

Matrix operator<=(const Matrix &A,Real d);

Matrix operator>=(const Matrix &A,Real d);

Matrix operator==(const Matrix &A,Real d);

Matrix operator!=(const Matrix &A,Real d);

Matrix minrows(const Matrix& A);   //min of each column (over the rows)

Matrix mincols(const Matrix& A);   //min of each row (over the columns)

Real min(const Matrix& A,Integer *iindex=0,Integer *jindex=0);

Matrix maxrows(const Matrix& A);   //similar

Matrix maxcols(const Matrix& A);   //similar

Real max(const Matrix& A,Integer *iindex=0,Integer *jindex=0);

Indexmatrix sortindex(const Matrix& vec,bool nondecreasing=true);

void sortindex(const Matrix& vec,Indexmatrix &ind, bool nondecreasing=true);

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

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

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

inline void Matrix::init_to_zero()
{
 nr=nc=0;mem_dim=0;m=0;
 aux_task=none;
 aux_dim=0;
 aux_m=0;
 chk_set_init(*this,1);
}

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

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

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

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

inline Matrix& Matrix::init(const std::vector<Real>& 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 Matrix& Matrix::init_diag(int n,Real d)
{
  init(n,n,0.);
  mat_xea(n,m,n+1,d);
  return *this;
}

inline Matrix& Matrix::init_diag(const Matrix& vec,Real d)
{
  Integer n=vec.dim();
  init(n,n,0.);
  if(d==1.)
    mat_xey(n,m,n+1,vec.get_store(),1);
  else
    mat_xeya(n,m,n+1,vec.get_store(),1,d);
  return *this;
}

inline Matrix& Matrix::init_diag(const Indexmatrix& vec,Real d)
{
  Integer n=vec.dim();
  init(n,n,0.);
  if(d==1.)
    for(Integer i=0;i<n;i++) m[i*(n+1)]=Real(vec.m[i]);
  else
    for(Integer i=0;i<n;i++) m[i*(n+1)]=vec.m[i]*d;
  return *this;
}


inline Matrix::Matrix():Memarrayuser()
{
 init_to_zero();
}

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

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

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

inline Matrix::Matrix(const Realrange& range):Memarrayuser()
{
 init_to_zero();
 init(range);
}

inline Matrix::Matrix(Integer inr,Integer inc,Real d):Memarrayuser()
{
 init_to_zero();
 init(inr,inc,d);
}

inline Matrix::Matrix(Integer inr,Integer inc,const Real *p,Integer incr,Real d):Memarrayuser()
{
 init_to_zero();
 init(inr,inc,p,incr,d);
}

inline Matrix::Matrix(const std::vector<Real>& vec):Memarrayuser()
{
 init_to_zero();
 init(vec);
}

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

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

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

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

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

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

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



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


inline Real colip(const Matrix& A,Integer j,const Matrix* scaling=0)
{
 chk_init(A);
 chk_range(j,j,A.nc,A.nc);
 if (scaling==0)
   return mat_ip(A.nr,A.m+j*A.nr);
 assert(scaling->dim()==A.nr);
 const Real* mp=A.m+j*A.nr;
 const Real* const mend=mp+A.nr;
 const Real* sp=scaling->m;
 Real sum=0;
 for(;mp!=mend;) {
    const Real d=(*mp++);
    sum+=d*d*(*sp++);
 }
 return sum;
}

inline Real rowip(const Matrix& A,Integer i,const Matrix* scaling=0)
{
 chk_init(A);
 chk_range(i,i,A.nr,A.nr);
 if (scaling==0)
   return mat_ip(A.nc,A.m+i,A.nr);
 assert(scaling->dim()==A.nc);
 const Real* mp=A.m+i;
 const Real* const mend=mp+A.nc*A.nr;
 const Real* sp=scaling->m;
 Real sum=0;
 for(;mp!=mend;mp+=A.nr) {
    const Real d=(*mp);
    sum+=d*d*(*sp++);
 }
 return sum;
}

inline Matrix colsip(const Matrix& A)
{
 chk_init(A);
 Matrix tmp(1,A.nc); chk_set_init(tmp,1);
 for (Integer j=0;j<A.nc;j++)
   tmp[j]=colip(A,j);
 return tmp;
}

inline Matrix rowsip(const Matrix& A)
{
 chk_init(A);
 Matrix tmp(A.nr,1,0.);
 Real* mp =A.m;
 for (Integer j=0;j<A.nc;j++){
   Real *tp=tmp.m;
   for (Integer i=0;i<A.nr;i++,mp++)
     (*tp++)+=(*mp)*(*mp);
 }
 return tmp;
}

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

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

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

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

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

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

Matrix& Matrix::operator/=(const Matrix &A)
{ chk_add(*this,A); mat_xinvhadey(nr*nc,m,A.m); return *this; }

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

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

inline Matrix& Matrix::operator/=(Real d)
{ chk_init(*this); mat_xmultea(nr*nc,m,1./d); return *this; }

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

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

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

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

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

inline Matrix operator%(const Matrix &A,const Matrix &B) 
    {Matrix C(A); return C%=B;}

inline Matrix operator/(const Matrix &A,const Matrix &B) 
    {Matrix C(A); return C/=B;}

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

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

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

inline Matrix operator+(const Matrix &A,Real d)           
    {Matrix B(A); return B+=d;}

inline Matrix operator+(Real d,const Matrix &A)           
    {Matrix B(A); return B+=d;}

inline Matrix operator-(const Matrix &A,Real d)       
    {Matrix B(A); return B-=d;}

inline Matrix operator-(Real d,const Matrix &A)       
    {Matrix B(A,-1.); return B+=d;}

inline int Matrix::QR_factor(Matrix& Q,Matrix& R,Real tol=1e-10) const  //*this is unchanged
{R=*this;return R.QR_factor(Q,tol);}
inline int Matrix::QR_factor(Matrix& Q,Matrix& R,Indexmatrix& piv,Real tol=1e-10) const
{R=*this;return R.QR_factor(Q,piv,tol);}
inline int Matrix::ls(Matrix & rhs, Real tol=1e-10)
    { return this->QR_solve(rhs,tol);}

inline int QR_factor(const Matrix& A,Matrix& Q,Matrix &R,Real tol)
{return A.QR_factor(Q,R,tol);}

inline int QR_factor(const Matrix& A,Matrix& Q,Matrix &R,Indexmatrix& piv,Real tol)
{return A.QR_factor(Q,R,piv,tol);}

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

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

inline Matrix concat_right(const Matrix& A,const Matrix& B)
          {Matrix C(A.dim()+B.dim(),1);C=A;C.concat_right(B);return C;}

inline Matrix concat_below(const Matrix& A,const Matrix& B)
          {Matrix C(A.dim()+B.dim(),1);C=A;C.concat_below(B);return C;}
    
inline void swap(Matrix &A,Matrix &B)
{ 
 Real *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 Matrix rand(Integer rows,Integer cols,CH_Tools::GB_rand* random_generator=0)
{Matrix A; return A.rand(rows,cols,random_generator);}

inline Matrix inv(const Matrix& A)
          {Matrix B(A); return B.inv();}

inline Matrix sqrt(const Matrix& A)
          {Matrix B(A); return B.sqrt();}

inline Matrix sqr(const Matrix& A)
          {Matrix B(A); return B.sqr();}

inline Matrix sign(const Matrix& A,Real tol=1e-12)
          {Matrix B(A); return B.sign(tol);}

inline Matrix floor(const Matrix& A)
          {Matrix B(A); return B.floor();}

inline Matrix ceil(const Matrix& A)
          {Matrix B(A); return B.floor();}

inline Matrix rint(const Matrix& A)
          {Matrix B(A); return B.floor();}

inline Matrix round(const Matrix& A)
          {Matrix B(A); return B.floor();}


inline Matrix operator>(const Matrix &A,const Matrix &B)
{return B<A;}

inline Matrix operator>=(const Matrix &A,const Matrix &B)
{return B<=A;}
inline Matrix operator<(Real d,const Matrix &A)
{return A>d;}
inline Matrix operator>(Real d,const Matrix &A)
{return A<d;}

inline Matrix operator<=(Real d,const Matrix &A)
{return A>=d;}

inline Matrix operator>=(Real d,const Matrix &A)
{return A<=d;}

inline Matrix operator==(Real d,const Matrix &A)
{return A==d;}

inline Matrix operator!=(Real d,const Matrix &A)
{return A!=d;}

inline bool equal(const Matrix& A,const Matrix& 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 Matrix& vec,bool nondecreasing)
{Indexmatrix ind;sortindex(vec,ind,nondecreasing);return ind;}


inline Indexmatrix find(const Matrix& A,Real tol=1e-10)
          {return A.find(tol);}
inline Indexmatrix find_number(const Matrix& A,Real num=0.,Real tol=1e-10)
          {return A.find_number(num,tol);}

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

}


//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__SYMMAT_HXX
#include "symmat.hxx"
#endif

#endif