linalg.h

basic linear algebra classes and applications (SVD,interpolation, multivariate optimization)

  ConstMatrixColumn(const ConstMatrixColumn&);	// Not implemented and forbidden:
  void operator = (const ConstMatrixColumn&);// no cloning/assignment allowed

protected:
				// Access the i-th element of the column
  inline const REAL& get_ref(const int index) const;

public:					// Take a col of the matrix
  ConstMatrixColumn(const Matrix& m, const int col);

  const REAL operator () (const int i) const
  			{ return get_ref(i); }
};
 
inline ElementWiseConst::ElementWiseConst(const ConstMatrixColumn& mc)
  	: start_ptr(const_cast<REAL*>(mc.col_ptr)),
  	  end_ptr(const_cast<REAL*>(mc.col_ptr)+mc.nrows)
  	{ }

				// Access the i-th element of the column
inline const REAL& ConstMatrixColumn::get_ref(const int index) const
{ register const int offset = index - row_lwb;
  if( offset >= nrows || offset < 0 )
      _error("MatrixColumn index %d is out of row boundaries [%d,%d]",
	      index,q_row_lwb(),q_row_upb());
  return col_ptr[offset];
}


inline ElementWiseConst of_every(const ConstMatrixColumn& mc)
		{ return mc; }

class MatrixColumn : public ConstMatrixColumn
{
  friend class ElementWise;
  friend class LAStreamOut;
  MatrixColumn(const MatrixColumn&);	// Not implemented and forbidden:
  void operator = (const MatrixColumn&);// no cloning/assignment allowed
public:					// Take a col of the matrix
  MatrixColumn(Matrix& m, const int col) :
  	ConstMatrixColumn(m,col) {}

  REAL& operator () (const int i)
  			{ return const_cast<REAL&>(get_ref(i)); }
};

inline ElementWise::ElementWise(const MatrixColumn& mc) : ElementWiseConst(mc) {}
inline ElementWise to_every(const MatrixColumn& mc) { return mc; }

			// A view of a single row of a matrix
class ConstMatrixRow : protected Matrix::ConstReference,
		       public DimSpec
{
  friend class ElementWiseStrideConst;
  friend class LAStrideStreamIn;
  friend class LAStrideStreamOut;
  const REAL * const row_ptr;		// Pointer to the row under
  					// consideration
  const int stride;			// if ptr=@a[row,i], then
					//    ptr+stride = @a[row,i+1]
					// Since elements of a[] are stored
					// col after col, stride = nrows
  const REAL * const end_ptr;		// Points after the end of the matrix
  
  ConstMatrixRow(const ConstMatrixRow&);	// Not implemented and forbidden:
  void operator = (const ConstMatrixRow&);	// no cloning/assignment allowed

protected:
				// Access the i-th element of the row
  inline const REAL& get_ref(const int index) const;

public:					// Take a row of the matrix
  ConstMatrixRow(const Matrix& m, const int row);

  const REAL operator () (const int i) const
  			{ return get_ref(i); }
};
 
inline
ElementWiseStrideConst::ElementWiseStrideConst(const ConstMatrixRow& mr)
  	: start_ptr(const_cast<REAL*>(mr.row_ptr)),
  	  end_ptr(const_cast<REAL*>(mr.end_ptr)),
  	  stride(mr.stride)
  	{ }

                                // The return statement of_every below
                                // would implicitly call a constructor
inline ElementWiseStrideConst of_every(const ConstMatrixRow& mr)
	{ return mr; }


				// Access the i-th element of the row
inline const REAL& ConstMatrixRow::get_ref(const int index) const
{ register const int offset = index - col_lwb;
  if( offset >= ncols || offset < 0 )
      _error("MatrixRow index %d is out of column boundaries [%d,%d]",
	      index,q_col_lwb(),q_col_upb());
  return row_ptr[stride*offset];
}


class MatrixRow : public ConstMatrixRow
{
  friend class ElementWiseStride;
  friend class LAStrideStreamOut;
  MatrixRow(const MatrixRow&);		// Not implemented and forbidden:
  void operator = (const MatrixRow&);	// no cloning/assignment allowed
public:					// Take a col of the matrix
  MatrixRow(Matrix& m, const int row) :
  	ConstMatrixRow(m,row) {}

  REAL& operator () (const int i)
  			{ return const_cast<REAL&>(get_ref(i)); }
};

inline ElementWiseStride::ElementWiseStride(const MatrixRow& mr)
	: ElementWiseStrideConst(mr) {}

inline ElementWiseStride to_every(const MatrixRow& mr) { return mr; }

				// A view of the matrix' main diagonal
				// ConstMatrixDiag is described as a Matrix
				// 1:n x 1:n, where n=min(nrows,ncols)
class ConstMatrixDiag : protected Matrix::ConstReference,
		        public DimSpec
{
  friend class ElementWiseStrideConst;
  friend class LAStrideStreamIn;
  friend class LAStrideStreamOut;
  const REAL * const start_ptr;		// Pointer to the upper left corner
  const int stride;			// if ptr=@a[i,i], then
					//    ptr+stride = @a[i+1,i+1]
					// Since elements of a[] are stored
					// col after col, stride = nrows+1
  const REAL * const end_ptr;		// Points after the end of the matrix
  
  ConstMatrixDiag(const ConstMatrixDiag&);	// Not implemented and forbidden:
  void operator = (const ConstMatrixDiag&);	// no cloning/assignment allowed

protected:
				// Access the i-th element of the diag
				// Note that numbering always starts with 1,
				// regardless of col_lwb and row_lwb of the
				// original matrix
  inline const REAL& get_ref(const int index) const;

public:
  explicit ConstMatrixDiag(const Matrix& m);

  const REAL operator () (const int i) const
  			{ return get_ref(i); }
};
 
inline
ElementWiseStrideConst::ElementWiseStrideConst(const ConstMatrixDiag& md)
  	: start_ptr(const_cast<REAL*>(md.start_ptr)),
  	  end_ptr(const_cast<REAL*>(md.end_ptr)),
  	  stride(md.stride)
  	{ }

                                // The return statement of of_every below
                                // would implicitly call a constructor
inline ElementWiseStrideConst of_every(const ConstMatrixDiag& md)
	{ return md; }


				// Access the i-th element of the diagonal
inline const REAL& ConstMatrixDiag::get_ref(const int index) const
{ if( index > ncols || index < 1 )
      _error("MatrixDiag index %d is out of diag boundaries [%d,%d]",
	      index,1,ncols);
  return start_ptr[stride*(index-1)];
}


class MatrixDiag : public ConstMatrixDiag
{
  friend class ElementWiseStride;
  friend class LAStrideStreamOut;
  MatrixDiag(const MatrixDiag&);	// Not implemented and forbidden:
  void operator = (const MatrixDiag&);	// no cloning/assignment allowed
public:					// Take a col of the matrix
  explicit MatrixDiag(Matrix& m) : ConstMatrixDiag(m) {}

  REAL& operator () (const int i)
  			{ return const_cast<REAL&>(get_ref(i)); }
};

inline ElementWiseStride::ElementWiseStride(const MatrixDiag& md)
	: ElementWiseStrideConst(md) {}

inline ElementWiseStride to_every(const MatrixDiag& md) { return md; }


/*
 *------------------------------------------------------------------------
 *			Inline Matrix operations
 */

inline Matrix::Matrix(const int no_rows, const int no_cols)
	: DimSpec(no_rows,no_cols)
{
  allocate();
}

inline Matrix::Matrix(const int row_lwb, const int row_upb,
		      const int col_lwb, const int col_upb)
	: DimSpec(row_lwb,row_upb,col_lwb,col_upb)
{
  allocate();
}

inline Matrix::Matrix(const DimSpec& dimension_specs)
	: DimSpec(dimension_specs)
{
  allocate();
}

inline Matrix::Matrix(const LazyMatrix& lazy_constructor)
	: DimSpec(lazy_constructor)
{
  allocate();
  lazy_constructor.fill_in(*this);
}


				// Force a promise of a matrix
				// That is, apply a lazy_constructor
				// to the current matrix
inline Matrix& Matrix::operator = (const LazyMatrix& lazy_constructor)
{
  is_valid();
  if( !(static_cast<const DimSpec&>(lazy_constructor) == *this) )
    info(),
    _error("The matrix above is incompatible with the assigned "
    	   "Lazy matrix");
      
  lazy_constructor.fill_in(*this);
  return *this;
}
					// Copy constructor, expensive: use
					// sparingly
inline Matrix::Matrix(const Matrix& another)
	: DimSpec(another)
{
  another.is_valid();
  allocate();
  *this = another;
}

				// Resize the matrix to new dimensions
inline void Matrix::resize_to(const int row_lwb, const int row_upb,
		 const int col_lwb, const int col_upb)
{
  resize_to(DimSpec(row_lwb,row_upb,col_lwb,col_upb));
}

inline const REAL& MatrixDABase::operator () (const int rown, const int coln) const
{
  register int arown = rown-row_lwb;		// Effective indices
  register int acoln = coln-col_lwb;

  if( arown >= nrows || arown < 0 )
    _error("Row index %d is out of Matrix boundaries [%d,%d]",
	   rown,row_lwb,nrows+row_lwb-1);
  if( acoln >= ncols || acoln < 0 )
    _error("Col index %d is out of Matrix boundaries [%d,%d]",
	   coln,col_lwb,ncols+col_lwb-1);
  
  return (index[acoln])[arown];
}

inline Matrix& Matrix::clear(void)	// Clean the Matrix
{
  is_valid();
  memset(elements,0,nelems*sizeof(REAL));
  return *this;
}

inline void are_compatible(const Matrix& m1, const Matrix& m2)
{
  m1.is_valid();
  m2.is_valid();
  
  if( !(static_cast<const DimSpec&>(m1) == m2) )
    m1.info(), m2.info(), _error("The matrices above are incompatible");
}

				// Assignment
inline Matrix& Matrix::operator = (const Matrix& source)
{
  are_compatible(*this,source);
  memcpy(elements,source.elements,nelems*sizeof(REAL));
  return *this;
}

				// Apply a user-defined action to each
				// element of a collection
inline ElementWiseConstAction&
ElementWiseConst::apply(ElementWiseConstAction& functor) const
{
  register const REAL * ep = start_ptr;
  while( ep < end_ptr )
    functor(*ep++);
  return functor;
}

inline ElementWiseAction&
ElementWise::apply(ElementWiseAction& functor)
{
  register REAL * ep = start_ptr;
  while( ep < end_ptr )
    functor(*ep++);
  return functor;
}

inline ElementWiseConstAction&
ElementWiseStrideConst::apply(ElementWiseConstAction& functor) const
{
  for(register const REAL * ep=start_ptr; ep < end_ptr; ep += stride )
    functor(*ep);
  return functor;
}

inline ElementWiseAction&
ElementWiseStride::apply(ElementWiseAction& functor)
{
  for(register REAL * ep=start_ptr; ep < end_ptr; ep += stride )
    functor(*ep);
  return functor;
}

			// Create a blank vector of a given size
			// (indexation starts at 1)
inline Vector::Vector(const int n)
	: Matrix(n,1) {}
  
			// Create a general lwb:upb vector blank vector
inline Vector::Vector(const int lwb, const int upb) : Matrix(lwb,upb,1,1) {}

			// Force the promise of a vector
inline Vector& Vector::operator = (const LazyMatrix& source)
  	{ Matrix::operator =(source); return *this; }

				// Resize the vector for a specified number
				// of elements, trying to keep intact as many
				// elements of the old vector as possible.
				// If the vector is expanded, the new elements
				// will be zeroes
inline void Vector::resize_to(const int n)	{ Vector::resize_to(1,n); }

inline void Vector::resize_to(const Vector& v)
{
  Vector::resize_to(v.q_lwb(),v.q_upb());
}

				// Get access to a vector element
inline REAL& Vector::operator () (const int ind)
{
  is_valid();
  register const int aind = ind - row_lwb;
  if( aind >= nelems || aind < 0 )
    _error("Requested element %d is out of Vector boundaries [%d,%d]",
	   ind,row_lwb,nrows+row_lwb-1);
  
  return elements[aind];
}

				// Make a vector following a recipe
inline Vector::Vector(const LazyMatrix& lazy_constructor)
  : Matrix(lazy_constructor)
{
  assure(ncols == 1,
         "cannot make a vector from a promise of a full matrix");
}

				// The following is just a syntactic sugar
				// to avoid typing of_every and to_every
inline Matrix& Matrix::operator = (const REAL val)
	{ to_every(*this) = val; return *this; }
	
inline Matrix& operator -= (Matrix& m, const double val)
	{ to_every(m) -= val; return m; }
inline Matrix& operator += (Matrix& m, const double val)
	{ to_every(m) += val; return m; }
inline Matrix& operator *= (Matrix& m, const double val)
	{ to_every(m) *= val; return m; }
	
inline Matrix::Reference& operator -= (Matrix::Reference& m, const double val)
	{ to_every(m) -= val; return m; }
inline Matrix::Reference& operator += (Matrix::Reference& m, const double val)
	{ to_every(m) += val; return m; }
inline Matrix::Reference& operator *= (Matrix::Reference& m, const double val)
	{ to_every(m) *= val; return m; }


inline bool operator ==  (const Matrix& m, const REAL val)
	{ return of_every(m) == val; }
inline bool operator !=  (const Matrix& m, const REAL val)
	{ return of_every(m) != val; }
inline bool operator ==  (const Matrix::ConstReference& m, const REAL val)
	{ return of_every(m) == val; }
inline bool operator !=  (const Matrix::ConstReference& m, const REAL val)
	{ return of_every(m) != val; }

#endif