matrix1.cc

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

// This may look like C code, but it is really -*- C++ -*-
/*
 ************************************************************************
 *
 *			Linear Algebra Package
 *
 *		Basic linear algebra operations, level 1
 *		      Element-wise operations
 *
 * $Id: matrix1.cc,v 4.3 1998/12/19 03:14:20 oleg Exp oleg $
 *
 ************************************************************************
 */

#ifdef __GNUC__
#pragma implementation "LinAlg.h"
#endif

#include "LAStreams.h"
#include <math.h>
#include "iostream.h"

/*
 *------------------------------------------------------------------------
 *			Constructors and destructors
 */

void Matrix::allocate(void)
{
  valid_code = MATRIX_val_code;

  assure( (nelems = nrows * ncols) > 0,
  	"The number of matrix cols and rows has got to be positive");

  assure( !ref_counter.q_engaged(),
  	"An attempt to allocate Matrix data which are in use" );
  
  name = "";

  assert( (elements = (REAL *)calloc(nelems,sizeof(REAL))) != 0 );

}


Matrix::~Matrix(void)		// Dispose of the Matrix struct
{
  is_valid();
  free(elements);
  if( name[0] != '\0' )
    delete const_cast<char*>(name);
  valid_code = 0;
}

				// Set a new Matrix name
void Matrix::set_name(const char * new_name)
{
  if( name != 0 && name[0] != '\0' )	// Dispose of the previous matrix name
    delete const_cast<char*>(name);

  if( new_name == 0 || new_name[0] == '\0' )
    name = "";				// Matrix is anonymous now
  else
    name = new char[strlen(new_name)+1],
      strcpy(const_cast<char *>(name), new_name);
}

				// Erase the old matrix and create a
				// new one according to new boundaries
				// with indexation starting at 1
void Matrix::resize_to(const int nrows, const int ncols)
{
  is_valid();
  if( nrows == Matrix::nrows && ncols == Matrix::ncols )
    return;

#if 0
  if( ncols != 1 )
    free(index);
#endif
  assert( !ref_counter.q_engaged() );
  free(elements);

  const char * const old_name = name;
  Matrix::nrows = nrows;
  Matrix::ncols = ncols;
  allocate();
  name = old_name;
}
				// Erase the old matrix and create a
				// new one according to new boundaries
void Matrix::resize_to(const DimSpec& dimension_specs)
{
  is_valid();
  Matrix::row_lwb = dimension_specs.q_row_lwb();
  Matrix::col_lwb = dimension_specs.q_col_lwb();
  resize_to(dimension_specs.q_nrows(),dimension_specs.q_ncols());
}

#if 0
					// Routing constructor module
Matrix::Matrix(const Matrix& A, const MATRIX_CREATORS_2op op, const Matrix& B)
{
  A.is_valid();
  B.is_valid();
  switch(op)
  {
    case Mult:
         _AmultB(A,B);
	 break;

    case TransposeMult:
	 _AtmultB(A,B);
	 break;

    default:
	 _error("Operation %d is not yet implemented",op);
  }
}
#endif


			// Build a column index to facilitate direct
			// access to matrix elements
REAL * const * MatrixDABase::build_index(const Matrix& m)
{
  m.is_valid();
  if( m.ncols == 1 )		// Only one col - index is dummy actually
    return const_cast<REAL * const *>(&m.elements);

  REAL ** const index = (REAL **)calloc(m.ncols,sizeof(REAL *));
  assert( index != 0 );
  register REAL * col_p = &m.elements[0];
  for(register REAL** ip = index; ip<index+m.ncols; col_p += m.nrows)
    *ip++ = col_p;
  return index;
}

MatrixDABase::~MatrixDABase(void)
{
  if( ncols != 1 )
    free((void *)index);
}

/*
 *------------------------------------------------------------------------
 * 		    Making a matrix of a special kind	
 */

				// Make a unit matrix
				// (Matrix needn't be a square one)
				// The matrix is traversed in the
				// natural (that is, col by col) order
Matrix& Matrix::unit_matrix(void)
{
  is_valid();
  register REAL *ep = elements;

  for(register int j=0; j < ncols; j++)
    for(register int i=0; i < nrows; i++)
        *ep++ = ( i==j ? 1.0 : 0.0 );

  return *this;
}

				// Make a Hilbert matrix
				// Hilb[i,j] = 1/(i+j-1), i,j=1...max, OR
				// Hilb[i,j] = 1/(i+j+1), i,j=0...max-1
				// (Matrix needn't be a square one)
				// The matrix is traversed in the
				// natural (that is, col by col) order
Matrix& Matrix::hilbert_matrix(void)
{
  is_valid();
  register REAL *ep = elements;

  for(register int j=0; j < ncols; j++)
    for(register int i=0; i < nrows; i++)
        *ep++ = 1./(i+j+1);

  return *this;
}

			// Create an orthonormal (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
			// E.g., the complete Haar matrix of the second order
			// is
			// column 1: [ 1  1  1  1]/2
			// column 2: [ 1  1 -1 -1]/2
			// column 3: [ 1 -1  0  0]/sqrt(2)
			// column 4: [ 0  0  1 -1]/sqrt(2)
			// Matrix m is assumed to be zero originally
void haar_matrix::fill_in(Matrix& m) const
{
  m.is_valid();
  assert( m.ncols <= m.nrows && m.ncols > 0 );
  register REAL * cp = m.elements;
     const REAL * const m_end = m.elements + m.nelems;

  double norm_factor = 1/sqrt((double)m.nrows);

  for(register int i=0; i<m.nrows; i++)	// First column is always 1
    *cp++ = norm_factor;	// (up to a normalization)

				// The other functions are kind of steps:
				// stretch of 1 followed by the equally
				// long stretch of -1
				// The functions can be grouped in families
				// according to their order (step size),
				// differing only in the location of the step
  int step_length = m.nrows/2;
  while( cp < m_end && step_length > 0 )
  {
    for(register int step_position=0; cp < m_end && step_position < m.nrows;
	step_position += 2*step_length, cp += m.nrows)
    {
      register REAL * ccp = cp + step_position;
      for(register int i=0; i<step_length; i++)
	*ccp++ = norm_factor;
      for(register int j=0; j<step_length; j++)
	*ccp++ = -norm_factor;
    }
    step_length /= 2;
    norm_factor *= sqrt(2.0);
  }
  assert( step_length != 0 || cp == m_end );
  assert( m.nrows != m.ncols || step_length == 0 );
}

haar_matrix::haar_matrix(const int order, const int no_cols)
    : LazyMatrix(1<<order,no_cols == 0 ? 1<<order : no_cols)
{
  assert(order > 0 && no_cols >= 0);
}

/*
 *------------------------------------------------------------------------
 * 			Collection-scalar arithmetics
 *	walk through the collection and do something with each
 *			element in turn, and the scalar
 */

				// For every element, do `elem OP value`
#define COMPUTED_VAL_ASSIGNMENT(OP,VALTYPE)				\
									\
void ElementWise::operator OP (const VALTYPE val)			\
{									\
  register REAL * ep = start_ptr;					\
  while( ep < end_ptr )							\
    *ep++ OP val;							\
}									\

COMPUTED_VAL_ASSIGNMENT(=,REAL)
COMPUTED_VAL_ASSIGNMENT(+=,double)
COMPUTED_VAL_ASSIGNMENT(-=,double)
COMPUTED_VAL_ASSIGNMENT(*=,double)

#undef COMPUTED_VAL_ASSIGNMENT


				// is "element OP val" true for all
				// elements in the collection?

#define COMPARISON_WITH_SCALAR(OP)					\
									\
bool ElementWiseConst::operator OP (const REAL val) const		\
{									\
  register const REAL * ep = start_ptr;					\
  while( ep < end_ptr )							\
    if( !(*ep++ OP val) )						\
      return false;							\
									\
  return true;								\
}									\


COMPARISON_WITH_SCALAR(==)
COMPARISON_WITH_SCALAR(!=)
COMPARISON_WITH_SCALAR(<)
COMPARISON_WITH_SCALAR(<=)
COMPARISON_WITH_SCALAR(>)
COMPARISON_WITH_SCALAR(>=)

#undef COMPARISON_WITH_SCALAR

/*
 *------------------------------------------------------------------------
 *	Apply algebraic functions to all elements of a collection
 */

				// Take an absolute value of a matrix
void ElementWise::abs(void)
{
  for(register REAL * ep=start_ptr; ep < end_ptr; ep++)
    *ep = ::abs(*ep);
}

				// Square each element
void ElementWise::sqr(void)
{
  for(register REAL * ep=start_ptr; ep < end_ptr; ep++)
    *ep = *ep * *ep;
}

				// Take the square root of all the elements
void ElementWise::sqrt(void)
{
  for(register REAL * ep=start_ptr; ep < end_ptr; ep++)
    if( *ep >= 0 )
      *ep = ::sqrt(*ep);
    else
      _error("%d-th element, %g, is negative. Can't take the square root",
	     (ep-start_ptr), *ep );
}

				// Apply a user-defined action to each matrix
				// element. The matrix is traversed in the
				// natural (that is, col by col) order
Matrix& Matrix::apply(ElementAction& action)
{
  is_valid();
  register REAL * ep = elements;
  for(register int j=col_lwb; j<col_lwb+ncols; j++)
    for(register int i=row_lwb; i<row_lwb+nrows; i++)
      action.operation(*ep++,i,j);
  assert( ep == elements+nelems );

  return *this;
}

const Matrix& Matrix::apply(ConstElementAction& action) const
{
  is_valid();
  register const REAL * ep = elements;
  for(register int j=col_lwb; j<col_lwb+ncols; j++)
    for(register int i=row_lwb; i<row_lwb+nrows; i++)
      action.operation(*ep++,i,j);
  assert( ep == elements+nelems );

  return *this;
}


/*
 *------------------------------------------------------------------------
 * 		Element-wise operations on two groups of elements
 */

				// Check to see if two matrices are identical
bool Matrix::operator == (const Matrix& m2) const
{
  are_compatible(*this,m2);
  return (memcmp(elements,m2.elements,nelems*sizeof(REAL)) == 0);
}

				// For every element, do `elem OP another.elem`
#define TWO_GROUP_COMP(OP)					\
									\
bool ElementWiseConst::operator OP (const ElementWiseConst& another) const \
{									\
  sure_compatible_with(another);					\
  register const REAL * sp = another.start_ptr;				\
  register const REAL * tp = start_ptr;					\
  while( tp < end_ptr )							\
    if( !(*tp++ OP *sp++) )						\
      return false;							\
									\
  return true;								\
}									\

TWO_GROUP_COMP(==)
TWO_GROUP_COMP(!=)
TWO_GROUP_COMP(<)
TWO_GROUP_COMP(<=)
TWO_GROUP_COMP(>)
TWO_GROUP_COMP(>=)

#undef TWO_GROUP_COMP

				// For every element, do `elem OP another.elem`
#define TWO_GROUP_OP(OP)					\
									\
void ElementWise::operator OP (const ElementWiseConst& another)		\
{									\
  sure_compatible_with(another);					\