linalg.h

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

 *	Data classes: Matrix and (its particular case) Vector
 *
 * Note that view classes (below) may "borrow" matrix data area
 * for a while; to make sure they "return" the data, we use a limited form
 * of reference counting. That is, the destruction/reallocation of a data
 * class are only possible when the reference counter is 1
 * (that is, there are no oustanding views). Otherwise it is a grave error.
 *
 */

class Matrix : public DimSpec
{
public:
  friend class ElementWiseConst;
  friend class ElementWiseStrideConst;
  friend class LAStreamIn;
  friend class LAStreamOut;
  friend class LAStrideStreamIn;
  friend class LAStrideStreamOut;
  friend class LABlockStreamIn;
  friend class LABlockStreamOut;
  friend class MatrixDABase;
  class ConstReference;			// Forward decl of nested classes
  class Reference;
  friend class ConstReference;
  
  friend class Vector;
  friend class ConstMatrixColumn;
  friend class ConstMatrixRow;
  friend class ConstMatrixDiag;

private:			// Private part
  int valid_code;			// Validation code
  enum { MATRIX_val_code = 5757 };	// Matrix validation code value
  RWWatchDog ref_counter;

protected:			// May be used in derived classes
  const char * name;				// Name for the matrix
  int nelems;				// Total no of elements, nrows*ncols
  REAL * elements;			// Elements themselves

  void allocate(void);

  friend void haar_matrix::fill_in(Matrix& m) const;
  friend void LazyTransposedMatrix::fill_in(Matrix& m) const;

public:			// Public interface

				// Constructors and destructors
					// Make a blank matrix
  Matrix(const int nrows, const int ncols);	// Indices start with 1
  Matrix(const int row_lwb, const int row_upb,	// Or with low:upper
	 const int col_lwb, const int col_upb);	// boundary specif.
  Matrix(const DimSpec& dimension_specs);
  inline Matrix(const Matrix&  another); // A real copy constructor, expensive


					// Construct a matrix applying a spec
					// operation to two prototypes
					// Example: Matrix A(10,12), B(12,5);
					// Matrix C(A,Matrix::Mult,B);
//  enum MATRIX_CREATORS_2op { Mult, 		// A*B
//			     TransposeMult, 	// A'*B
//			     InvMult, 		// A^(-1)*B
//			     AtBA }; 		// A'*B*A
//  Matrix(const Matrix& A, const MATRIX_CREATORS_2op op, const Matrix& B);
//  Matrix(const Vector& x, const Vector& y);	// x'*y (diad) matrix

  Matrix(const LazyMatrix& lazy_constructor);//Make a matrix using given recipe
  explicit Matrix(const char * file_name);	// Read the matrix from the file
					// (not yet implemented!)

  ~Matrix(void);			// Destructor

  void set_name(const char * name);	// Set a new matrix name

					// Erase the old matrix and create a
					// new one according to new boundaries
  void resize_to(const int nrows, const int ncols);	// Indexation @ 1
  void resize_to(const int row_lwb, const int row_upb,	// Or with low:upper
		 const int col_lwb, const int col_upb);	// boundary specif.
  void resize_to(const DimSpec& dimension_specs);	// Like other matrix m


  void is_valid(void) const
  { if( valid_code != MATRIX_val_code ) invalid_matrix(); }

				// Status inquires
  int q_no_elems(void) const		{ return nelems; }
  const char * q_name(void) const	{ return name; }

  class ConstReference
  {
    friend class Reference;
    Matrix& m;
    ConstReference(const ConstReference&);	// Not implemented and forbidden
    void operator = (const ConstReference&);	// reference isn't transferable
  public:
    ConstReference(const Matrix& _m)
    	: m(const_cast<Matrix&>(_m))
    	{ m.is_valid(); m.ref_counter.get_shared(); }
   ~ConstReference(void) { m.is_valid(); m.ref_counter.release_shared(); }
    /*explicit*/ operator const Matrix& (void) const { return m; }
    const Matrix& ref(void) const { return m; }
  };

  class Reference : public ConstReference
  {
    Reference(const Reference&);	// Not implemented and forbidden
    void operator = (const Reference&);	// reference isn't transferable
  public:
    Reference(Matrix& _m) : ConstReference(_m) {}
    /*explicit*/ operator Matrix& (void) const { return m; }
    Matrix& ref(void) const { return m; }
  };

					// Apply a user-defined action
  Matrix& apply(ElementAction& action);		// to each matrix element
  const Matrix& apply(ConstElementAction& action) const; // to each matrix element

  					// Although the following method does
  					// not need a priviledged access to a
  					// Matrix, it is required to be a member
  					// of the class
  inline Matrix& operator = (const REAL val);
					// Invert the matrix returning the
					// determinant if desired
					// determinant = 0 if the matrix is
					// singular
					// If determ_ptr=dummy_determinant_ref
					// and the matrix *is*
					// singular, throw up
  static double dummy_determinant_ref;
  Matrix& invert(double &determ_ref = dummy_determinant_ref);

				// Element-wise operations on two matrices
  inline Matrix& operator = (const Matrix& source);	// Assignment
  inline Matrix& operator = (const ConstReference& ref)
  		{ return *this = ref.ref(); }
  Matrix& operator = (const LazyMatrix& source);// Force the promise of a matrix
  Matrix& clear(void);			// Clear the matrix (fill out with 0)
 

					// Comparisons
  bool operator == (const Matrix& im2) const;
  friend inline void are_compatible(const Matrix& im1, const Matrix& im2);


				// True matrix operations
				// (on a matrix as a whole)  
  Matrix& operator *= (const Matrix& source);	// Inplace multiplication
						// possible only for square src
  inline Matrix& operator *= (const ConstReference& ref)
  		{ return *this *= ref.ref(); }

  
  void mult(const Matrix& A, const Matrix& B);  // Compute A*B
 
				// Compute matrix norms
  double row_norm(void) const;		// MAX{ SUM{ |M(i,j)|, over j}, over i}
  double norm_inf(void) const		// Alias, shows the norm is induced
         { return row_norm(); }		// 	by the vector inf-norm
  double col_norm(void) const;		// MAX{ SUM{ |M(i,j)|, over i}, over j}
  double norm_1(void) const		// Alias, shows the norm is induced
         { return col_norm(); }		// 	by the vector 1-norm
  double e2_norm(void) const;		// SUM{ m(i,j)^2 }, Note it's square
					// of the Frobenius rather than 2. norm

  friend double e2_norm(const Matrix& m1, const Matrix& m2);
					// e2_norm(m1-m2)

  double determinant(void) const;	// Matrix must be a square one

  double asymmetry_index(void) const;	// How far is the matrix from being
					// symmetrical (0 means complete symm)
					// (not yet implemented)

				// Make matrices of a special kind
  Matrix& unit_matrix(void);		// Matrix needn't be a square
  Matrix& hilbert_matrix(void);		// Hilb[i,j] = 1/(i+j-1); i,j=1..max

				// I/O: write, read, print 
  					// Write to a file
					// "| command name" is OK as a file
					// name
  void write(const char * file_name,const char * title = "") const;
  void info(void) const;		// Print the info about the Matrix
  void print(const char title []) const;// Print the Matrix as a table
  volatile void invalid_matrix(void) const;   // The matrix in question is not valid
  					// ==> crash the program
};

void compare(const Matrix& im1, const Matrix& im2, const char title[]);


/*
 *------------------------------------------------------------------------
 *		   Vector as a n*1 matrix (that is, a col-matrix)
 */
class Vector : public Matrix
{
public:
  Vector(const int n);		// Create a blank vector of a given size
				// (indexation starts at 1)
				// Create a general lwb:upb vector blank vector
  Vector(const int lwb, const int upb);
//  Vector(const Vector& another);	// a copy constructor (to be inferred)

					// Make a vector and assign init vals
  Vector(const int lwb, const int upb,  // lower and upper bounds
	 double iv1, ...    		// DOUBLE numbers. The last arg of
	 ); 				// the list must be string "END"
					// E.g: Vector f(1,2,0.0,1.5,"END");

  Vector(const LazyMatrix& lazy_constructor);//Make a vector using given recipe

				// Resize the vector (keeping as many old
				// elements as possible), expand by zeros
  void resize_to(const int n);			// Indexation @ 1
  void resize_to(const int lwb, const int upb); 	// lwb:upb specif
  void resize_to(const Vector& v);			// like other vector

  				// Unlike Matrix, Vector permits a direct
  				// access to its elements without further ado
  				// Still, streams/of_every() are more efficient
  				// for uniform access
  inline REAL& operator () (const int index);
  inline REAL  operator () (const int index) const
	{ return (const_cast<Vector&>(*this))(index); }

				// Listed below are specific vector operations
				// (unlike n*1 matrices)

					// Status inquires
  int q_lwb(void) const		{ return q_row_lwb(); }
  int q_upb(void) const		{ return q_row_upb(); }

					// Compute the scalar product
  friend double operator * (const Vector& v1, const Vector& v2);

					// "Inplace" multiplication
					// target = A*target
					// A needn't be a square one (the
					// target will be resized to fit)
  Vector& operator *= (const Matrix& A);

					// Vector norms
  inline double norm_1(void) const;   		// SUM{ |v[i]| }
  inline double norm_2_sqr(void) const;	       	// SUM{ v[i]^2 }
  inline double norm_inf(void) const;		// MAX{ |v[i]| }

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

  					// Although the following methods do
  					// not need a priviledged access to
  					// Vector, they are required to be
  					// members of the class
  inline Vector& operator = (const ElementWiseConst& stream);
  inline Vector& operator = (const ElementWiseStrideConst& stream);
  inline Vector& operator = (const REAL val)
  	{ Matrix::operator =(val); return *this; }
};

				// Service functions (useful in the
				// verification code). They print some detail
				// info if the validation condition fails
void verify_element_value(const Matrix& m,const REAL val);
void verify_matrix_identity(const Matrix& m1, const Matrix& m2);
inline void verify_element_value(const Matrix::ConstReference& m,const REAL val)
	{ verify_element_value(m.ref(),val); }
inline void verify_matrix_identity(const Matrix::ConstReference& m1,
				   const Matrix& m2)
	{ verify_matrix_identity(m1.ref(),m2); }
inline void verify_matrix_identity(const Matrix& m1,
				   const Matrix::ConstReference& m2)
	{ verify_matrix_identity(m1,m2.ref()); }
inline void verify_matrix_identity(const Matrix::ConstReference& m1,
				   const Matrix::ConstReference& m2)
	{ verify_matrix_identity(m1.ref(),m2.ref()); }

				// Multiply by the diagonal
				// of another matrix
Matrix& operator *= (Matrix& m, const ConstMatrixDiag& diag);

				// Compute the product of two matrices
class LazyMatrixProduct : public LazyMatrix
{
  const Matrix& A;
  const Matrix& B;
  void fill_in(Matrix& m) const { m.mult(A,B); }
  LazyMatrixProduct(const Matrix& _A, const Matrix& _B)
	: LazyMatrix(_A.q_row_lwb(),_A.q_row_upb(),
		     _B.q_col_lwb(),_B.q_col_upb()),
	  A(_A), B(_B) {}
  
public:
  friend const LazyMatrixProduct operator * (const Matrix& A, const Matrix& B)
    { return LazyMatrixProduct(A,B); }
};

				// Make a zero matrix
class LazyZeroMatrix : public LazyMatrix
{
  void fill_in(Matrix& m) const { __unused(Matrix& m1) = m; }
  LazyZeroMatrix(const Matrix& proto)
    : LazyMatrix(static_cast<const DimSpec&>(proto)) {}
public:
  friend const LazyZeroMatrix zero(const Matrix& proto)
    { proto.is_valid(); return proto; }
};

				// Make a unit matrix
class LazyUnitMatrix : public LazyMatrix
{
  void fill_in(Matrix& m) const { m.unit_matrix(); }
  LazyUnitMatrix(const Matrix& proto)
    : LazyMatrix(static_cast<const DimSpec&>(proto)) {}
public:
  friend const LazyUnitMatrix unit(const Matrix& proto)
    { proto.is_valid(); return proto; }
};

				// Make an inverse matrix
class LazyInverseMatrix : public LazyMatrix
{
  const Matrix& proto;
  void fill_in(Matrix& m) const { m = proto; m.invert(); }
  LazyInverseMatrix(const Matrix& _proto)
    : LazyMatrix(static_cast<const DimSpec&>(_proto)),
      proto(_proto) {}
public:
  friend const LazyInverseMatrix inverse(const Matrix& proto)
    { proto.is_valid(); return proto; }
};


				// Make a transposed matrix
inline LazyTransposedMatrix::LazyTransposedMatrix(const Matrix & _proto)
  	: LazyMatrix(_proto.q_col_lwb(),_proto.q_col_upb(),
  		     _proto.q_row_lwb(),_proto.q_row_upb()),
  	  proto(_proto) {}
inline const LazyTransposedMatrix transposed(const Matrix & proto)
  	{ proto.is_valid(); return proto; }

/*
 *------------------------------------------------------------------------
 *	Matrix Views: provide access to elements of a Matrix
 *	(or groups of elements: rows, columns, diagonals, etc)
 * Although views can be created and exist for some time, one should
 * not attempt to pass around views (especially return them as a result
 * of the function) when the data object itself (a matrix) could be
 * disposed of. Thanks to the reference counting, this situation would
 * be detected, and the program would crash.
 */
 
/*
 *------------------------------------------------------------------------
 */

/*
 * Functions/procedures (generally grouped together under a ElementWise
 * class umbrella) that perform a variety of operation on each element
 * of some structure in turn.
 *
 * ElementWiseConst and ElementWise classes are largely a syntactic
 * sugar; you can't explicitly construct an object of this class, on stack
 * or on heap. The only way the class can be used is within a pattern
 *      to_every(matrix) = 1;
 *	if( of_every(matrix) > 1 ) ...
 * etc. Besides, this is the only way it makes sense.
 *
 * Note that the ElementWise classes are friends of Matrix, so they are
 * privy of Matrix internals and can use them to make processing faster.
 */


				// A functor describing an operation to apply 
				// to every element of a collection in
				// turn
				// This operation receives only the elements's
				// value, and thus it may not modify the
				// element under consideration.
class ElementWiseConstAction
{
				// Those aren't implemented; but making them
				// private forbids the assignement
  ElementWiseConstAction(const ElementWiseConstAction&);
  ElementWiseConstAction& operator= (const ElementWiseConstAction&);
public:
  ElementWiseConstAction(void) {}
  virtual void operator () (const REAL element) = 0;
};

				// The same as above but this functor is given
				// a reference to an element, which it can
				// modify
class ElementWiseAction
{
				// Those aren't implemented; but making them
				// private forbids the assignement
  ElementWiseAction(const ElementWiseAction&);
  ElementWiseAction& operator= (const ElementWiseAction&);
public:
  ElementWiseAction(void) {}
  virtual void operator () (REAL& element) = 0;