linalg.h

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

};

class ElementWiseConst
{
  friend class ElementWise;
  friend class ElementWiseStrideConst;
  friend class ElementWiseStride;

  REAL * const start_ptr;		// Pointer to the beginning of
  					// the group of elements
  REAL * const end_ptr;			// Points after the end of the group
  
  void operator=(const ElementWiseConst&);// is not implemented, ergo, is
                                        // not allowed
  ElementWiseConst(const ElementWiseConst&);// Cloning is not allowed, either
  
				// A private constructor, to make
				// sure the object can't be constructed
				// and left, while the matrix would disappear
				
				// Use of_every() friend to do
				// the actual construction
  ElementWiseConst(const Matrix& m)
  	: start_ptr(const_cast<REAL*>(m.elements)),
  	  end_ptr(const_cast<REAL*>(m.elements+m.nelems))
  	{ m.is_valid(); }

  inline ElementWiseConst(const ConstMatrixColumn& mc);
  
  				// Just a helper class to disambiguate 
  				// type conversion from ElementWiseConst
  				// to ElementWiseStrideConst and others
  struct ElementWiseConstRef
  {
    const ElementWiseConst& ref;
    ElementWiseConstRef(const ElementWiseConstRef&);
    void operator = (const ElementWiseConstRef&);
    ElementWiseConstRef(const ElementWiseConst& _ref) : ref(_ref) {}
  };

public:

					// No public constructors...

                                        // ElementWiseConst(Matrix& m) would be
                                        // implicitly called
  friend inline ElementWiseConst of_every(const Matrix& m) { return m; }
  friend inline ElementWiseConst of_every(const Matrix::ConstReference& m)
  	{ return m.ref(); }
  
  friend inline ElementWiseConst of_every(const ConstMatrixColumn& mc);
 				
  				// Every ElementWiseConst is a
  				// ElementWiseStrideConst with a stride of 1
  inline operator ElementWiseStrideConst (void);
  
				// Comparisons
				// Find out if the predicate
				// "element op val" is true for ALL
				// elements of the group
  bool	  operator ==  (const REAL val) const;	// ? all elems == val
  bool	  operator !=  (const REAL val) const;	// ? all elems != val
  bool	  operator <   (const REAL val) const;	// ? all elems <  val
  bool	  operator <=  (const REAL val) const;	// ? all elems <= val
  bool	  operator >   (const REAL val) const;	// ? all elems >  val
  bool	  operator >=  (const REAL val) const;	// ? all elems >= val
  
				// Find out if the predicate
				// "element op another.elem" is true for ALL
				// elements of the two groups
  bool	  operator ==  (const ElementWiseConst& another) const;
  bool	  operator !=  (const ElementWiseConst& another) const;
  bool	  operator <   (const ElementWiseConst& another) const;
  bool	  operator <=  (const ElementWiseConst& another) const;
  bool	  operator >   (const ElementWiseConst& another) const;
  bool	  operator >=  (const ElementWiseConst& another) const;
				
  void sure_compatible_with(const ElementWiseConst& another) const
  	{ assure(end_ptr-start_ptr == another.end_ptr-another.start_ptr,
  		"Two groups of elements to operate on have different sizes"); }
  
  				// Data reduction: reduce a collection
  				// to one value: a "norm" of the collection
  double sum(void) const;
  double sum_squares(void) const;
  double sum_abs(void) const;
  double max_abs(void) const;
  //double foldr(const double seed, Reductor& reductor);
  
  				// Reduce a difference between two collections
  				// to a single value: a "norm" of the
  				// difference
  double sum_squares(const ElementWiseConst& another) const;
  double sum_abs(const ElementWiseConst& another) const;
  double max_abs(const ElementWiseConst& another) const;
  //double foldr(const double seed, Reductor2& reductor);
  
  				// Just let the user do what he wants
  				// with each element in turn
  				// (without modifying them, of course)
  ElementWiseConstAction& apply(ElementWiseConstAction& functor) const;
};


class ElementWise : public ElementWiseConst
{
  void operator=(const ElementWise&);	// is not implemented, ergo, is
                                        // not allowed
  ElementWise(const ElementWise&);	// Cloning is not allowed, either
  
				// A private constructor, to make
				// sure the object can't be constructed
				// and left, while the matrix would disappear
				
				// Use a to_every() friend to do
				// the actual construction
  ElementWise(Matrix& m) : ElementWiseConst(m) {}
  inline ElementWise(const MatrixColumn& mc);
public:

					// No public constructors...

                                        // ElementWise(Matrix& m) would be
                                        // implicitly called
  friend inline ElementWise to_every(Matrix& m) { return m; }
  friend inline ElementWise to_every(const Matrix::Reference& m)
  	{ return m.ref(); }
//  friend inline ElementWiseConst of_every(const ConstMatrixColumn& mc);
  friend inline ElementWise to_every(const MatrixColumn& mc); // see below

   				// Every ElementWise is a
  				// ElementWiseStride with a stride of 1
  inline operator ElementWiseStride (void);
  				// The same as the conversion operator above,
  				// use when a compiler (gcc 2.8.1) fails to
  				// apply it
  inline ElementWiseStride with_stride (void);

				// group-scalar arithmetics
				// Modify every element of the group
				// according to the operation
  void operator =   (const REAL val);	// Assignment to all the elems
  void operator +=  (const double val);	// Add to elements
  void operator -=  (const double val);	// Take from elements
  void operator *=  (const double val);	// Multiply elements by a val

				// Other element-wise matrix operations
  void abs(void);			// Take an absolute value of a matrix
  void sqr(void);			// Square each element
  void sqrt(void);			// Take the square root

  void operator = (const ElementWiseConst& another);
  void operator += (const ElementWiseConst& another);
  void operator -= (const ElementWiseConst& another);
  void operator *= (const ElementWiseConst& another);
  void operator /= (const ElementWiseConst& another);

  ElementWiseAction& apply(ElementWiseAction& functor);
};

inline Vector& Vector::operator = (const ElementWiseConst& stream)
  	{ to_every(*this) = stream; return *this; }


inline double Vector::norm_1(void) const
	{ return of_every(*this).sum_abs(); }
inline double Vector::norm_2_sqr(void) const
	{ return of_every(*this).sum_squares(); }
inline double Vector::norm_inf(void) const
	{ return of_every(*this).max_abs(); }

//---------------------------------------------------------------------------
//			The same, only with a stride

class ElementWiseStrideConst
{
  friend class ElementWiseStride;
  friend class ElementWiseConst;

  REAL * const start_ptr;		// Pointer to the beginning of
  					// the group of elements
  REAL * const end_ptr;			// Points after the end of the group
  const int stride;
  
  void operator=(const ElementWiseStrideConst&);// is not implemented, ergo, is
                                        // not allowed
  ElementWiseStrideConst(const ElementWiseStrideConst&);// Cloning is not allowed, either
  
				// A private constructor, to make
				// sure the object can't be constructed
				// and left, while the matrix would disappear
				
				// Use of_every() friend to do
				// the actual construction
  ElementWiseStrideConst(const ElementWiseConst::ElementWiseConstRef& ewc)
  	: start_ptr(ewc.ref.start_ptr),
  	  end_ptr(ewc.ref.end_ptr),
  	  stride(1)
  	  {}
  	  
  inline ElementWiseStrideConst(const ConstMatrixRow& mr);
  inline ElementWiseStrideConst(const ConstMatrixDiag& md);

public:

					// No public constructors...
 
 friend inline ElementWiseStrideConst of_every(const ConstMatrixRow& mr);
 friend inline ElementWiseStrideConst of_every(const ConstMatrixDiag& md);

				// Comparisons
				// Find out if the predicate
				// "element op val" is true for ALL
				// elements of the group
  bool	  operator ==  (const REAL val) const;	// ? all elems == val
  bool	  operator !=  (const REAL val) const;	// ? all elems != val
  bool	  operator <   (const REAL val) const;	// ? all elems <  val
  bool	  operator <=  (const REAL val) const;	// ? all elems <= val
  bool	  operator >   (const REAL val) const;	// ? all elems >  val
  bool	  operator >=  (const REAL val) const;	// ? all elems >= val
  
				// Find out if the predicate
				// "element op another.elem" is true for ALL
				// elements of the two groups
  bool	  operator ==  (const ElementWiseStrideConst& another) const;
  bool	  operator !=  (const ElementWiseStrideConst& another) const;
  bool	  operator <   (const ElementWiseStrideConst& another) const;
  bool	  operator <=  (const ElementWiseStrideConst& another) const;
  bool	  operator >   (const ElementWiseStrideConst& another) const;
  bool	  operator >=  (const ElementWiseStrideConst& another) const;
				
//  void sure_compatible_with(const ElementWiseStrideConst& another) const
//  	{ assure(end_ptr-start_ptr == another.end_ptr-another.start_ptr,
//  		"Two groups of elements to operate on have different sizes"); }
  
  				// Data reduction: reduce a collection
  				// to one value: a "norm" of the collection
  double sum(void) const;
  double sum_squares(void) const;
  double sum_abs(void) const;
  double max_abs(void) const;
  //double foldr(const double seed, Reductor& reductor);
  
  				// Reduce a difference between two collections
  				// to a single value: a "norm" of the
  				// difference
  double sum_squares(const ElementWiseStrideConst& another) const;
  double sum_abs(const ElementWiseStrideConst& another) const;
  double max_abs(const ElementWiseStrideConst& another) const;
  //double foldr(const double seed, Reductor2& reductor);
  
  				// Just let the user do what he wants
  				// with each element in turn
  				// (without modifying them, of course)
  ElementWiseConstAction& apply(ElementWiseConstAction& functor) const;
};
  
class ElementWiseStride : public ElementWiseStrideConst
{
  friend class ElementWise;
  void operator=(const ElementWiseStride&);	// is not implemented, ergo, is
                                        // not allowed
  ElementWiseStride(const ElementWiseStride&);	// Cloning is not allowed, either
  
				// A private constructor, to make
				// sure the object can't be constructed
				// and left, while the matrix would disappear
				
				// Use a to_every() friend to do
				// the actual construction
  inline ElementWiseStride(const MatrixRow& mr);
  inline ElementWiseStride(const MatrixDiag& md);
  ElementWiseStride(const ElementWiseConst::ElementWiseConstRef& ewc)
  	: ElementWiseStrideConst(ewc) {}
public:

					// No public constructors...

                                // The return statement of the friends below
                                // would implicitly call a constructor
  friend inline ElementWiseStride to_every(const MatrixRow& mr);
  friend inline ElementWiseStride to_every(const MatrixDiag& md);
 
				// group-scalar arithmetics
				// Modify every element of the group
				// according to the operation
  void operator =   (const REAL val);	// Assignment to all the elems
  void operator +=  (const double val);	// Add to elements
  void operator -=  (const double val);	// Take from elements
  void operator *=  (const double val);	// Multiply elements by a val

				// Other element-wise matrix operations
  void abs(void);			// Take an absolute value of a coll
  void sqr(void);			// Square each element
  void sqrt(void);			// Take the square root

  void operator = (const ElementWiseStrideConst& another);
  void operator += (const ElementWiseStrideConst& another);
  void operator -= (const ElementWiseStrideConst& another);
  void operator *= (const ElementWiseStrideConst& another);
  void operator /= (const ElementWiseStrideConst& another);

  ElementWiseAction& apply(ElementWiseAction& functor);
};

inline ElementWiseConst::operator ElementWiseStrideConst (void)
	{ return ElementWiseConstRef(*this); }

inline ElementWise::operator ElementWiseStride (void)
	{ return ElementWiseConstRef(*this); }

inline ElementWiseStride ElementWise::with_stride (void)
	{ return ElementWiseConstRef(*this); }

inline Vector& Vector::operator = (const ElementWiseStrideConst& stream)
  	{  to_every(*this).with_stride() = stream; return *this; }


/*
 *------------------------------------------------------------------------
 *	MatrixDA: a view of a matrix that provides a direct access
 *	to its elements. A column index is built to speed up access
 */

				// A mixin class that builds a row index
				// for a matrix grid
class MatrixDABase : public DimSpec
{
  REAL * const * const index;		// index[i] = &matrix(0,i) (col index)
  MatrixDABase(const MatrixDABase&);	// Not implemented and forbidden:
  void operator = (const MatrixDABase&);// no cloning/assignment allowed
  
  static REAL * const * build_index(const Matrix& m);

protected:
  MatrixDABase(const Matrix& m)
        : DimSpec(m),
          index( build_index(m) )
          {}
  inline const REAL& operator () (const int rown, const int coln) const;
public:
  ~MatrixDABase(void);
};

class ConstMatrixDA : public Matrix::ConstReference, public MatrixDABase
{
  ConstMatrixDA(const ConstMatrixDA&);	// Not implemented and forbidden:
  void operator = (const ConstMatrixDA&);// no cloning/assignment allowed
public:
  ConstMatrixDA(const Matrix& m)
        : Matrix::ConstReference(m), MatrixDABase(m)
          {}
				// Individual element manipulations
  const REAL operator () (const int rown, const int coln) const
  			{ return MatrixDABase::operator() (rown,coln); }
  const REAL operator () (const rowcol rc) const
  		{ return operator () (rc.row,rc.col); }
  				// I wish the cast to const Matrix& worked,
  				// as given by the following operator. Alas,
  				// an explicit member function seems to
  				// be needed	
//  const Matrix& get_container(void) const
//  	{ return Matrix::ConstReference::operator const Matrix& (); }
};
 
class MatrixDA : public Matrix::Reference, public MatrixDABase
{
  MatrixDA(const MatrixDA&);		// Not implemented and forbidden:
  void operator = (const MatrixDA&);	// no cloning/assignment allowed

public:
  MatrixDA(Matrix& m)
        : Matrix::Reference(m), MatrixDABase(m)
          {}
				// Individual element manipulations
  REAL& operator () (const int rown, const int coln) const
  	{ return const_cast<REAL&>(MatrixDABase::operator()(rown,coln)); }
  REAL& operator () (const rowcol rc) const
  		{ return operator () (rc.row,rc.col); }
  	
  				// I wish the cast to Matrix& worked by itself,
  				// as given by the following operator. Alas,
  				// an explicit member function seems to
  				// be needed
#if defined (__GNUC__)
  Matrix& get_container(void) { return Matrix::Reference::operator Matrix& (); }
#else
  Matrix& get_container(void) { return *this; }
#endif
};
 
/*
 *------------------------------------------------------------------------
 *		MatrixRow, MatrixCol, MatrixDiag
 * They are special kind of views, of selected parts of the matrix
 * they are implemented as special kind of streams that walk only selected
 * parts of the matrix
 */

class ConstMatrixColumn : protected Matrix::ConstReference,
			  public DimSpec
{
  friend class ElementWiseConst;
  friend class LAStreamIn;
  friend class LAStreamOut;
  const REAL * const col_ptr;		// Pointer to the column under
  					// consideration