linalg.h

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

// This may look like C code, but it is really -*- C++ -*-
/*
 ************************************************************************
 *
 *			  Linear Algebra Package
 *
 * The present package implements all the basic algorithms dealing
 * with vectors, matrices, matrix columns, etc.
 * Matrix is a basic object in the package; vectors, symmetric matrices,
 * etc. are considered matrices of a special type.
 *
 * Matrix elements are arranged in memory in a COLUMN-wise
 * fashion (in FORTRAN's spirit). In fact, this makes it very easy to
 * pass the matrices to FORTRAN procedures implementing more
 * elaborate algorithms.
 *
 * Matrix and vector indices always start with 1, spanning up to 
 * a given limit, unless specified otherwise by a user. 
 *
 * The present package provides all facilities to completely AVOID returning
 * matrices. If one really needs to return a matrix, return
 * a LazyMatrix object instead. The conversion is completely transparent
 * to the end user, e.g. "Matrix m = haar_matrix(5);" and _is_ efficient.
 *
 * Container and views: object of a class Matrix is the only one that
 * owns the data, that is, a 2D grid of values. All other classes provide
 * different views of this grid (as a long 1D array of a matrix stretched
 * column-by-column, a view of a single matrix row, column, or a diagonal,
 * a view of a certain matrix block, etc).
 *
 * The views are designed to be as lightweight as possible: most of
 * the view functionality could be inlined. That is why views do not
 * have any virtual functions. The view classes are supposed to be
 * final, in JAVA parlance. It makes little sense to inherit from them.
 * It is unnecessary as views are designed to be flexible and efficient,
 * they provide a wide variety of access to the matrix with safety and
 * as efficient as possible. So flexible that many former Matrix methods
 * can be implemented outside of the Matrix class as they no longer need
 * a priviledged access to Matrix internals, without any loss of performance.
 *
 * Importance of Matrix streams: a sequential view/access to a matrix or
 * its parts. Many of Linear Algebra algorithms actually require only
 * sequential access to a matrix or its rows/columns, which is simpler and
 * faster than a random access. That's why LinAlg stresses streams. There are
 * two kinds of sequential views: ElementWise are transient views that
 * typically exist only for the duration of a elementwise action.
 * That's why it doesn't require a "locking" of a matrix. On the other hand,
 * LAStream and the ilk are more permanent views of the matrix, with
 * roughly the same functionality. An application itself can traverse the
 * stream at its own convenience, bookmarking certain spots and possibly
 * returning to them later. LAStream does assert a shared lock on the matrix,
 * to prevent the element matrix from moving/disappearing.
 *
 * $Id: LinAlg.h,v 4.3 1998/12/22 21:02:55 oleg Exp oleg $
 *
 ************************************************************************
 */

#ifndef __GNUC__
#pragma once
#endif
#ifndef _LinAlg_h
#define _LinAlg_h 1

#if defined(__GNUC__)
#pragma interface
#define __unused(X) X __attribute__((unused))
class ostream;			// An Opaque class
#else
#include <iostream.h>	  // When ostream is a template, can't make it opaque
#define __unused(X) X
#endif

#include "myenv.h"
#include "std.h"
#if defined(__GNUC__)
#include <values.h>
#else
#include <limits.h>
#define MAXINT INT_MAX
#endif
#include <math.h>
#include "builtin.h"
#include "minmax.h"
#if defined(index)
#undef index
#endif

/*
 *------------------------------------------------------------------------
 *			Auxiliary classes
 */

typedef float REAL;			// Scalar field of the Linear Vector
					// space

			// Dimensions specifier of a 2D object
class DimSpec
{
protected:
  int nrows;				// No. of rows
  int ncols;				// No. of columns
  int row_lwb;				// Lower bound of the row index
  int col_lwb;				// Lower bound of the col index

public:

  DimSpec(const int _nrows, const int _ncols)	// Indices start with 1
  	: nrows(_nrows), ncols(_ncols), row_lwb(1), col_lwb(1)
  	  { assert(nrows>0 && ncols>0); }
  DimSpec(const int _row_lwb, const int _row_upb,	// Or with low:upper
	  const int _col_lwb, const int _col_upb)	// boundary specif.
   	: nrows(_row_upb-_row_lwb+1), 
   	  ncols(_col_upb-_col_lwb+1),
   	  row_lwb(_row_lwb), col_lwb(_col_lwb)
  	  { assert(nrows>0 && ncols>0); }
  
				// Status inquires
  int q_row_lwb(void) const			{ return row_lwb; }
  int q_row_upb(void) const			{ return nrows+row_lwb-1; }
  int q_nrows(void) const			{ return nrows; }
  int q_col_lwb(void) const			{ return col_lwb; }
  int q_col_upb(void) const			{ return ncols+col_lwb-1; }
  int q_ncols(void) const			{ return ncols; }
  
  bool operator == (const DimSpec& another_dim_spec) const
  	{ return nrows == another_dim_spec.nrows &&
  		 ncols == another_dim_spec.ncols &&
  		 row_lwb == another_dim_spec.row_lwb &&
  		 col_lwb == another_dim_spec.col_lwb; }
 
  friend ostream& operator << (ostream& os, const DimSpec& dimspec);
};

				// A 2D index, a pair representing
				// a row,col location in some 2D structure
struct rowcol
{
  int row, col;
  rowcol(const int _row, const int _col) : row(_row), col(_col) {}
  bool operator == (const rowcol another) const
  	{ return row == another.row && col == another.col; }
  friend ostream& operator << (ostream& os, const rowcol& rc);
};

		// A Range of indices: lwb:upb
		// if lwb > upb, the range is "empty"
		// If lwb is not specified, it defaults to -IRange::INF
		// (meaning the range spans from the earliest possible
		// element in a collection the range is being applied to)
		// If upb is not specified, it defaults to IRange::INF
		// meaning the range spans through the end of the collection,
		// whatever that may be.
struct IRange
{
  const int lwb, upb;
  enum { INF = MAXINT };
  IRange(const int _lwb, const int _upb) : lwb(_lwb), upb(_upb) {}
  static IRange from(const int lwb) { return IRange(lwb,INF); }
  static IRange through(const int upb) { return IRange(-INF,upb); }
  bool q_empty(void) const	{ return lwb > upb; }
  bool q_proper(void) const
  	{ return abs(lwb) != INF && abs(upb) != INF && lwb <= upb; }
  friend ostream& operator << (ostream& os, const IRange range);
  int no_less_than(const int strict_lwb) const
   	{ if( lwb == -IRange::INF ) return strict_lwb;
	  assert(lwb >= strict_lwb); return lwb; }
  int no_more_than(const int strict_upb) const
   	{ if( upb == IRange::INF ) return strict_upb;
	  assert(upb <= strict_upb); return upb; }
};

		// The same as above, but with a given stride
struct IRangeStride
{
  const int lwb, upb, stride;
  IRangeStride(const int _lwb, const int _upb, const int _stride)
  	: lwb(_lwb), upb(_upb), stride(_stride) {}
  IRangeStride(const IRange range)
  	: lwb(range.lwb), upb(range.upb), stride(1) {}
  bool q_empty(void) const	{ return lwb > upb; }
  bool q_proper(void) const
  	{ return abs(lwb) != IRange::INF && abs(upb) != IRange::INF
		&& lwb <= upb && stride > 0; }
};


		// A cut from DimSpec in two dimensions
class DimSpecSubranged : public DimSpec
{
protected:
   int min_offset;	// An offset to a[imin,jmin] from the original DimSpec
   int max_offset;	// An offset to a[imax,jmax] from the original DimSpec
   const int original_nrows;	// The number of rows in the original DimSpec

public:
 DimSpecSubranged(const DimSpec& orig, const IRange row_range,
		  const IRange col_range)
 	: DimSpec(row_range.no_less_than(orig.q_row_lwb()),
		  row_range.no_more_than(orig.q_row_upb()),
		  col_range.no_less_than(orig.q_col_lwb()),
		  col_range.no_more_than(orig.q_col_upb())),
	  original_nrows(orig.q_nrows())
	  { min_offset = (col_lwb - orig.q_col_lwb())*original_nrows +
	  		 (row_lwb - orig.q_row_lwb());
	    max_offset = min_offset + (ncols-1)*original_nrows + nrows - 1;
	  }
		// we are sure row_lwb is within [orig.row_lwb,orig.row_upb]
	  	// and so is the new row_upb (and the same for col)

  int q_min_offset(void) const		{ return min_offset; }
  int q_max_offset(void) const		{ return max_offset; }
  int q_row_diff(void) const		{ return original_nrows - nrows; }
};
				// An interface describing an operation to apply
				// to a matrix element.
				// This interface is used by Matrix::apply()
				// and similar for-each-type iterators. They
				// call a concrete instance of this interface
				// on each element of a matrix etc. collection
				// under consideration.
				// The operation receives the (i,j) index of
				// the current element matrix and its _value_:
				// thus this interface may not modify the
				// element itself.
class ConstElementAction
{
public:
  virtual void operation(const REAL element, const int i, const int j) = 0;
};

				// The same as above, only an operation
				// is allowed to modify a matrix element
				// it is being applied to.
class ElementAction
{
public:
  virtual void operation(REAL& element, const int i, const int j) = 0;

//private:			// Those aren't implemented; but making them
				// private forbids the assignement
//  ElementAction(const ElementAction&);
//  void operator= (const ElementAction&);
};


				// An interface that performs an operation
				// on two elements of Matrices being jointly
				// traversed. The operation receives the values
				// of the two current elements in both matrices
				// and the element's index (which must be
				// the same for each matrix involved).
				// As interface recieves elements' _values_
				// it may not modify the elements.
				// This interface is being used to compute
				// some property of two matrices being jointly
				// traversed.
class ConstElementAction2
{
public:
  virtual void operation(const REAL element1, const REAL element2,
			 const int i, const int j) = 0;
};

				// Lazy matrix constructor
				// That is, a promise of a matrix rather than
				// a matrix itself. The promise is forced
				// by a Matrix constructor or assignment
				// operator, see below.
				// It's highly recommended a function never
				// return Matrix itself. Use LazyMatrix
				// instead
class LazyMatrix : public DimSpec
{
  friend class Matrix;
  virtual void fill_in(Matrix& m) const = 0;

				// Not implemented: cloning is prohibited
  //  LazyMatrix(const LazyMatrix&);
  LazyMatrix& operator = (const LazyMatrix&);
public:
  LazyMatrix(const int nrows, const int ncols)	// Indices start with 1
    : DimSpec(nrows,ncols) {}
  LazyMatrix(const int row_lwb, const int row_upb,// Or with low:upper
	     const int col_lwb, const int col_upb)// boundary specif.
    : DimSpec(row_lwb,row_upb,col_lwb,col_upb) {}
  LazyMatrix(const DimSpec& dims) : DimSpec(dims) {}
};

				// A pair specifying extreme values
				// of some domain/result
class MinMax
{
  REAL min_val, max_val;		// Invariant: max_val >= min_val
public:
  MinMax(REAL _min_val, REAL _max_val)
    : min_val(_min_val), max_val(_max_val) {}
  explicit MinMax(REAL val)
    : min_val(val), max_val(val) {}
  MinMax& operator << (REAL val)
  	{ if( val > max_val ) max_val = val;
	  else if(val < min_val) min_val = val; return *this; }
  REAL min(void) const		{ return min_val; }
  REAL max(void) const		{ return max_val; }
  double ratio(void) const	{ return max_val/min_val; }
  friend ostream& operator << (ostream& os, const MinMax& mx);
};

                // The following class watches over exclusive (write)
                // and shared (read) access to a data structure the object
                // is a member of.
                // The ref_count field tells how the object is being
                // currently accessed:
                //      = 0: the object is not engaged
                //      =-1: exclusive (write access)
                //      >0 : the object is being currently read (viewed)
                //      by one or several readers
                // The watchdog watches over obvious blunders as trying
                // to get a read access while the object is being exclusively
                // held, or trying to grab the write access or destroy the
                // object while it is engaged in some other activity.
                // Note, if an object contains a pointer to some data
                // structure, using this pointer requires a read (shared)
                // access, even if the data structure would be modified.
                // Indeed, as long as all 'readers' access the data structure
                // at the same place in memory, they certainly are manipulating
                // the same (and the real) thing. Modifying a data structure
                // _pointer_ on the other hand requires an exclusive access
                // to the object (pointer holder). Thus read/write access
                // is meant a shared/exclusive access to
                // object's data and object's layout (but not the contents
                // of the data itself)
class RWWatchDog
{
  int ref_count;                        // see explanation above

  RWWatchDog(const RWWatchDog&);        // Deliberately unimplemented:
  void operator = (const RWWatchDog&);  // no cloning allowed!

  volatile void access_violation(const char reason []);

public:
  bool q_engaged(void) const            { return ref_count != 0; }
  bool q_shared(void) const             { return ref_count > 0; }
  bool q_exclusive(void) const          { return ref_count == -1; }

  void get_exclusive(void)
        { if( q_engaged() ) access_violation("get_exclusive"); ref_count = -1; }
  void release_exclusive(void)
        { if(!q_exclusive()) access_violation("release_exclusive");
          ref_count = 0; }

  void get_shared(void)
        { if( q_exclusive() ) access_violation("get_shared"); ref_count++; }
  void release_shared(void)
        { if( !q_shared() ) access_violation("release_shared"); ref_count--; }

  RWWatchDog(void) : ref_count(0) {}
  ~RWWatchDog(void)
  	{ if( ref_count!=0 ) access_violation("destroying"); 
  	  ref_count = -1; }

  				// Dump the current status
  friend ostream& operator << (ostream& os, const RWWatchDog& wd);
};

#if 0
class MinMaxLoc : public MinMax
{
};
#endif


/*
 *------------------------------------------------------------------------
 *	Special kinds of matrices that have to be forward-declared
 */
 
class Matrix;
class MatrixColumn;
class ConstMatrixColumn;
class MatrixRow;
class MatrixDiag;
				// Create an orthogonal (2^n)*(no_cols) Haar
				// (sub)matrix, whose columns are Haar
				// functions
				// If no_cols is 0, create the complete matrix
				// with 2^n columns
class haar_matrix : public LazyMatrix
{
  void fill_in(Matrix& m) const;
public:
  haar_matrix(const int n, const int no_cols=0);
};

				// Lazy transposition
class LazyTransposedMatrix : public LazyMatrix
{
  const Matrix& proto;
  void fill_in(Matrix& m) const;
  inline LazyTransposedMatrix(const Matrix & _proto);
public:
  friend inline const LazyTransposedMatrix transposed(const Matrix & proto);
};

/*
 *------------------------------------------------------------------------