//
// Matrix.cpp
//
// This file contains function definitions for the Matrix class.
// This code is part of the Bioinformatics Template Library (BTL).
//
// Copyright (C) 1997 Birkbeck College, Malet Street, London WC1E 7HX, U.K.
// (classlib@mail.cryst.bbk.ac.uk)
//
// This library is free software; you can
// redistribute it and/or modify it under the terms of
// the GNU Library General Public License as published by the Free
// Software Foundation; either version 2 of the License, or (at your
// option) any later version.  This library is distributed in the hope
// that it will be useful, but WITHOUT ANY WARRANTY; without even the
// implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE.  See the GNU Library General Public License for more details.
// You should have received a copy of the GNU Library General Public
// License along with this library; if not, write to the Free Software
// Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
////////////////////////////////////////////////////////////////////////////


#include "Typedefn.h"
#include "MathVector.h"
#include "Matrix.h"
#include "NumericLimits.h"
#include "NumericUtilities.h"
#include <math.h>
#include <stdlib.h>
#include <float.h>
#include <iomanip.h>

//............................................................................
// Subroutine to calculate the single value decomposition of a Matrix:
//
// A = U * L * Vtrans
//
// where A is this matrix. 
// 
// Output:
//  	Matrix U is written over this Matrix, V is output, L is a diagonal
//  	matrix and is output as a Vector. 

void Matrix::SVD(Vector& L, Matrix& V)
{
    if (L.size() != ncols)
    	FATAL_ERROR("The size of the input Vector must equal the number of "
    	            "columns in this Matrix.");

    if (V.nrows != ncols || V.ncols != ncols)
    	FATAL_ERROR("The number of columns and rows of the input Matrix must " 
    	            "equal the number of columns in this Matrix.");
    
    Matrix ATA = (*this) & (*this);
    ATA.Eigens(L,V);
    
    // Find inverse of Matrix with L at its diagonal.
    for (size_type i=1; i<=ncols; i++)
    	L(i) = sqrt(L(i));
    
    // U = A * V * Linv
    // Linv is a diagonal matrix and is represented here as a Vector.
    // N.B. x /= y is often less efficient than x *= 1/y
    //
    (*this) = (*this) * V;
    for (size_type row=1; row<=nrows; row++)
    	for (size_type col=1; col<=ncols; col++)
    	    (*this)(row,col) *= 1/L(col);
}

//.............................................................................
// Subroutine to compute eigenvalues & eigenvectors of a real symmetric
// matrix, from IBM SSP manual (see p165).
// 
// 		       Ian Tickle April 1992
// 
// 		(modified by David Moss February 1993)
// 
// Eigenvalues & vectors of real symmetric matrix stored in triangular form.
// 
// Arguments:
// 
// mv = 0 to compute eigenvalues only.
// mv = 1 to compute eigenvalues & eigenvectors.
// 
// n - supplied dimension of n*n matrix.
// 
// a - supplied array of size n*(n+1)/2 containing n*n matrix in the order:
// 
// 	   1      2      3    ...
//     1    a[0]
//     2    a[1]   a[2]
//     3    a[3]   a[4]   a[5]   ...
//     ...
//     NOTE a is used as working space and is overwritten on return.
//     Eigenvalues are written into diagonal elements of a
//     i.e.  a[0]  a[2]  a[5]  for a 3*3 matrix.
// 
// r - Resultant matrix of eigenvectors stored columnwise in the same
//     order as eigenvalues.

void Matrix::Eigen(int mv, int n, long double* a, long double* r) const
{ 
    int i, il, ilq, ilr, im, imq, imr, ind, iq, j, l, ll, lm, lq, m, mm, mq;
  long double anorm, anrmx, cosx, cosx2, sincs, sinx, sinx2, thr, x, y;
  long double theshold = LDBL_EPSILON;

  if (mv)
  { for (i=1; i<n*n; i++) r[i]=0.;
    for (i=0; i<n*n; i+=n+1) r[i]=1.;
  }

/* Initial and final norms (anorm & anrmx). */
  anorm=0.;
  iq=0;
  for (i=0; i<n; i++) for (j=0; j<=i; j++)
  { if (j!=i) anorm+=a[iq]*a[iq];
    ++iq;
  }

  if (anorm>0.)
  { anorm=sqrt(2.*anorm);
    anrmx=theshold*anorm/n;

/* Compute threshold and initialise flag. */
    thr=anorm;
    do
    { thr/=n;
      do
      { ind=0;
	l=0;
	do
	{ lq=l*(l+1)/2;
	  ll=l+lq;
	  m=l+1;
	  ilq=n*l;
	  do

/* Compute sin & cos. */
	  { mq=m*(m+1)/2;
	    lm=l+mq;
	    if (fabs(a[lm])>=thr)
	    { ind=1;
	      mm=m+mq;
	      x=.5*(a[ll]-a[mm]);
	      y=-a[lm]/sqrt(a[lm]*a[lm]+x*x);
	      if (x<0.) y=-y;
	      sinx=y/sqrt(2.*(1.+(sqrt(1.-y*y))));
	      sinx2=sinx*sinx;
	      cosx=sqrt(1.-sinx2);
	      cosx2=cosx*cosx;
	      sincs=sinx*cosx;

/* Rotate l & m columns. */
	      imq=n*m;
	      for (i=0; i<n; i++)
	      { iq=i*(i+1)/2;
		if (i!=l && i!=m)
		{ if (i<m) im=i+mq;
		  else im=m+iq;
		  if (i<l) il=i+lq;
		  else il=l+iq;
		  x=a[il]*cosx-a[im]*sinx;
		  a[im]=a[il]*sinx+a[im]*cosx;
		  a[il]=x;
		}
		if (mv)
		{ ilr=ilq+i;
		  imr=imq+i;
		  x=r[ilr]*cosx-r[imr]*sinx;
		  r[imr]=r[ilr]*sinx+r[imr]*cosx;
		  r[ilr]=x;
		}
	      }

	      x=2.*a[lm]*sincs;
	      y=a[ll]*cosx2+a[mm]*sinx2-x;
	      x=a[ll]*sinx2+a[mm]*cosx2+x;
	      a[lm]=(a[ll]-a[mm])*sincs+a[lm]*(cosx2-sinx2);
	      a[ll]=y;
	      a[mm]=x;
	    }

/* Tests for completion.
   Test for m = last column. */
	  } while (++m!=n);

/* Test for l =penultimate column. */
	} while (++l!=n-1);
      } while (ind);

/* Compare threshold with final norm. */
    } while (thr>anrmx);
  }
}       // Matrix::eigen

//.............................................................................
// Sort eigenvalues & eigenvectors (if mv=1) in order of descending eigenvalue.
// Arguments are as for eigen, which must have been called.

void Matrix::Esort(int mv, int n, long double* a, long double* r) const
{
  int i, im, j, k, km, l, m;
  long double am;

  k=0;
  for (i=0; i<n-1; i++)
  { im=i;
    km=k;
    am=a[k];
    l=0;
    for (j=0; j<n; j++)
    { if (j>i && a[l]>am)
      { im=j;
	km=l;
	am=a[l];
      }
      l+=j+2;
    }
    if (im!=i)
    { a[km]=a[k];
      a[k]=am;
      if (mv)
      { l=n*i;
	m=n*im;
	for (j=0; j<n; j++)
	{ am=r[l];
	  r[l++]=r[m];
	  r[m++]=am;
	}
      }
    }
    k+=i+2;
  }
}   // Matrix::Esort

//.............................................................................
// Constructs a unit matrix of size i*i.

Matrix::Matrix(const size_type& s) : nrows(s), ncols(s)
{
    if (s<1) FATAL_ERROR("The dimensions must be greater than 0");
    
    size_type size = s*s;
    mat = new value_type[size];
    
    if (mat == 0) FATAL_ERROR("Not enough memory");
    
    size_type diagpos = s+1;
    for (size_type i=0, flag=0; i<size; i++, flag++)
    {
    	if (flag == diagpos) flag=0;
    	if (flag != 0)
    	    mat[i] = 0.0;
    	else
    	    mat[i] = 1.0; 
    }
}	

//.............................................................................
// Constructor for default 3x3 matrix and initialise elements to zero.

Matrix::Matrix() : nrows(3), ncols(3)
{
    mat = new value_type[9];
    if (mat == 0) FATAL_ERROR("Not enough memory");

    mat[0]=0.0; mat[1]=0.0; mat[2]=0.0; mat[3]=0.0; mat[4]=0.0;
    mat[5]=0.0; mat[6]=0.0; mat[7]=0.0; mat[8]=0.0;
}

//.............................................................................
// Constructor with row and column arguments. Initialises each element to zero

Matrix::Matrix(const size_type& p, const size_type& q): nrows(p),ncols(q)
{
    size_type size = p*q;

    if (p<1 || q<1) FATAL_ERROR("Both dimensions must be greater than 0");

    mat = new value_type[size];

    if (mat == 0) FATAL_ERROR("Not enough memory");

    for (size_type i = 0; i < size; i++)
    	mat[i] = 0.0;
}

//.............................................................................
// Construct from a two dimensional array of reals ( given value_type* )

Matrix::Matrix(const value_type * const array,
               const size_type& p,
               const size_type& q) : nrows(p), ncols(q)
{
    if (p<1 || q<1) FATAL_ERROR("Both dimensions must be greater than 0");

    size_type size = p*q;
    mat = new value_type[size];

    if (mat == 0) FATAL_ERROR("Not enough memory");

    for (size_type i = 0; i < size; i++)
    	mat[i] = array[i];
} 

//.............................................................................
// Constructor for 3x3 matrix with initialisation

Matrix::Matrix
    	    (const value_type& e1, const value_type& e2, const value_type& e3,
             const value_type& e4, const value_type& e5, const value_type& e6,
             const value_type& e7, const value_type& e8, const value_type& e9) 
             : nrows(3), ncols(3)
{ 
    mat = new value_type[9];

    if (mat == 0) FATAL_ERROR("Not enough memory");

    mat[0]=e1; mat[1]=e2; mat[2]=e3; mat[3]=e4; mat[4]=e5; mat[5]=e6;
    mat[6]=e7; mat[7]=e8; mat[8]=e9;
}
//..............................................................................
// Constructor for 3x3 rotation matrix for a rotation about a given axis by a
// given number of degrees. Create the following matrix:
//
//  [cosA+a1a1(1-cosA)      -a3sinA+a1a2(1-cosA) 	a2sinA+a1a2(1-cosA) ] 
//  [a3sinA+a1a2(1-cosA)    cosA+a2a2(1-cosA)   	-a1sinA+a2a3(1-cosA)] 
//  [-a2sinA+a3a1(1-cosA)   a1sinA+(1-cosA)a2a3 	cosA+(1-cosA)a3a3   ]
//
// where a is the unit vector in the direction of the axis and A is the angle
// of rotation.

Matrix::Matrix(const Vector& axis, const value_type& angle) : nrows(3), ncols(3)
{
    mat = new value_type[9];
    if (mat == 0) FATAL_ERROR("Not enough memory");

    const double DEGRAD = M_PI/180.0; 	// converts degrees to radians
    double radAngle = angle * DEGRAD;
    double cosA = cos(radAngle);
    double mCosA = 1-cosA;    	
    double sinA = sin(radAngle);
    Vector a = axis/axis.Modulus();
    
    mat[0] =  cosA  	    +   a(1)*a(1) * mCosA;
    mat[1] = -a(3) * sinA   + 	a(1)*a(2) * mCosA;
    mat[2] =  a(2) * sinA   +	a(1)*a(3) * mCosA;
    mat[3] =  a(3) * sinA   +	a(1)*a(2) * mCosA;
    mat[4] =  cosA	    +	a(2)*a(2) * mCosA;
    mat[5] = -a(1) * sinA   +	a(2)*a(3) * mCosA;
    mat[6] = -a(2) * sinA   +	a(1)*a(3) * mCosA;
    mat[7] =  a(1) * sinA   +	a(2)*a(3) * mCosA;
    mat[8] =  cosA	    +	a(3)*a(3) * mCosA;
}

//.............................................................................
// Copy constructor

Matrix::Matrix(const Matrix &m) : nrows(m.nrows), ncols(m.ncols)
{
    mat = new value_type[nrows * ncols];
    if (mat == 0) FATAL_ERROR("Not enough memory");
    for (size_type i=0; i<nrows * ncols; i++)
    	mat[i] = m.mat[i];
}

//.............................................................................
// matrix assignment

Matrix& Matrix::operator=(const Matrix& m)
{
    // REPLACE mat by stl vector OR OPTIMISE
    if (this == &m) return *this;

    delete[] mat;
    mat = new value_type[m.nrows * m.ncols];
    if (mat == 0) FATAL_ERROR("Not enough memory");

    nrows = m.nrows;
    ncols = m.ncols;
    for (size_type i=0; i<nrows * ncols; i++) 
    	mat[i] = m.mat[i];
    return *this;
}

//.............................................................................
// Read only access a matrix element given its indices e.g. x(i,j) i = row,
// j = col N.B. (i,j >= 1)
	
Matrix::value_type Matrix::operator()(const size_type& i,
                                      const size_type& j) const
{
	
#if defined(DEBUG_VERSION)
    // Bounds checks
    //
    if (i<1 || j<1 || i>nrows || j>ncols) 
     	FATAL_ERROR("Matrix indices out of range");
#endif

    return mat[(ncols*(i-1) + j - 1)];
}

//.............................................................................
// Read/write access a matrix element given its indices e.g. x(i,j) i = row,
// j = col N.B. (i,j >= 1)
	
Matrix::value_type& Matrix::operator()(const size_type& i,
                                       const size_type& j)
{
	
#if defined(DEBUG_VERSION)
    // Bounds checks
    //
    if (i<1 || j<1 || i>nrows || j>ncols) 
    	FATAL_ERROR("Matrix indices out of range");
#endif

    return mat[(ncols*(i-1) + j - 1)];
}

//.............................................................................
// Returns an iterator that points the first elment of a given row
// N.B. ( 0 <= i < nrows )"] */

Matrix::iterator Matrix::operator[](const size_type& i)
{ 

#if defined(DEBUG_VERSION)
    // Bounds checks
    //
    if (i>=nrows) 
    	FATAL_ERROR("Matrix index out of range");
#endif

    return (&mat[0] + ncols*i);
}

//.............................................................................
// Returns a const_iterator that points the first elment of a given row
// N.B. ( 0 <= i < nrows )"] */

Matrix::const_iterator Matrix::operator[](const size_type& i) const
{ 

#if defined(DEBUG_VERSION)
    // Bounds checks
    //
    if (i>=nrows) 
    	FATAL_ERROR("Matrix index out of range");
#endif

    return (&mat[0] + ncols*i);
}


//.............................................................................
// Matrix multiplication

Matrix Matrix::operator*(const Matrix& m) const
{
    if (ncols != m.nrows) 
    	FATAL_ERROR("The number of rows in the input Matrix must equal the "
    	            "number of columns in this Matrix");
 
    Matrix answer(nrows,m.ncols);
    for (size_type i=1; i<=nrows; i++)
    	for (size_type j=1; j<=m.ncols; j++)
    	    for (size_type k=1; k<=m.nrows; k++)
	    	answer(i,j) += (*this)(i,k) * m(k,j);

    return answer;
}

//.............................................................................
// Premultiplication of a vector by a matrix

Vector Matrix::operator*(const Vector &v) const
{
    if (ncols != v.size()) 
    	FATAL_ERROR("The size of the input Vector must equal the number of "
    	            "columns in this Matrix");

    Vector v1(nrows);
    for(size_type i=1, iarr=0; i<=nrows; i++, iarr=ncols*(i-1))
    	for(size_type j=1; j<=ncols; j++)
    	    // v1(i) = *this(i,j) * v(j)
    	    v1(i) += mat[iarr + j-1] * v(j);  
    return v1;

}

//.............................................................................
// Postmultiplication of a vector by a matrix

Vector operator*(const Vector &v, const Matrix &m)
{
    typedef Matrix::size_type size_type;
    size_type size = min(m.ReadNRows(),v.size());

    if (v.size() != m.ReadNRows())
     	FATAL_ERROR("The size of the input Vector must equal the number of rows"
     	            " in this Matrix");

    Vector v1(size);
    for(size_type i=1, iarr=0; i<=size; i++, iarr=i-1)
    	for(size_type j=1; j<=size; j++)
    	    // v1(i)+=v(j)*m(j,i)
    	    v1(i) += v(j) * m.mat[(m.ncols*(j-1) + iarr)]; 

    return v1;
}

//.............................................................................
// multiple matrix by a value_type

Matrix& Matrix::operator*=(const value_type& r)
{
    for (size_type i=0; i<size(); i++)
	mat[i] *= r;
    return *this;
}

//.............................................................................
// Multiple matrix by a value_type

Matrix Matrix::operator*(const value_type& r) const
{
    Matrix m1 = *this;
    m1 *= r;
    return m1;
}

//.............................................................................
// Multiple matrix by a value_type

Matrix operator*(const Matrix::value_type& r, const Matrix& m)
{
    Matrix m1(m);
    m1 *= r;
    return m1;
}

//.............................................................................

// Divide matrix by a value_type
Matrix& Matrix::operator/=(const value_type& r)
{
#if defined(DEBUG_VERSION)
    if (r == 0.0) WARNING("Divide by zero");
#endif

    for (size_type i=0; i<size(); i++)
	    mat[i] /= r;
    return *this;
}

//.............................................................................
// Divide matrix by a value_type

Matrix Matrix::operator/(const value_type& r) const
{
#if defined(DEBUG_VERSION)
    if (r == 0.0) WARNING("Divide by zero");
#endif
    Matrix m1 = *this;
    m1 /= r;
    return m1;
}

//..............................................................................
// Subtraction of a number from each element

Matrix& Matrix::operator-=(const value_type &r)
{
    value_type *vp = &mat[0];
    for (size_type i=0; i<size(); i++)
    	*vp++ -= r;
    return *this;
}    	

//..............................................................................
// Subtraction of a number from each element.

Matrix Matrix::operator-(const value_type &r) const
{
    Matrix m1 = *this;
    m1 -= r;
    return m1;
}    	

//..............................................................................
// Addition of a number to each element.

Matrix& Matrix::operator+=(const value_type &r)
{
    value_type *vp = &mat[0];
    for (size_type i=0; i<size(); i++)
    	*vp++ += r;
    return *this;
}    	

//..............................................................................
// Addition of a number to each element.

Matrix Matrix::operator+(const value_type &r) const
{
    Matrix m1 = *this;
    m1 += r;
    return m1;
}    	

//.............................................................................
// Subtraction of a vector from the rows of a matrix

Matrix& Matrix::operator-=(const Vector &v)
{

    if (v.size() != ncols)
    	FATAL_ERROR("The the size of the input Vector must equal the number of "
    	            "columns in this Matrix");

    for (size_type i=0, ipos=0; i<nrows; i++)
    	for (size_type j=0; j<ncols; j++,ipos++)
    	    mat[ipos] -= v[j];
    	    
    return *this;
}
    	    
//.............................................................................
// Subtraction of a vector from the rows of a matrix

Matrix Matrix::operator-(const Vector &v) const
{
    Matrix m1 = *this;
    m1 -= v;
    return m1;

}  //  Matrix::operator-

//.............................................................................
// Addition of a vector from the rows of a matrix

Matrix& Matrix::operator+=(const Vector &v)
{

    if (v.size() != ncols)
    	FATAL_ERROR("The the size of the input Vector must equal the number of "
    	            "columns in this Matrix");

    for (size_type i=0, ipos=0; i<nrows; i++)
    	for (size_type j=0; j<ncols; j++,ipos++)
    	    mat[ipos] += v[j];
    	    
    return *this;
}
    	    
//.............................................................................
// Addition of a vector from the rows of a matrix

Matrix Matrix::operator+(const Vector &v) const
{
    Matrix m1 = *this;
    m1 += v;
    return m1;

}

//..............................................................................
// Matrix subtraction. The the size of Matrix m must equal the size of this 
// Matrix.

Matrix& Matrix::operator-=(const Matrix &m)
{
    if (nrows != m.nrows || ncols != m.ncols)
    	FATAL_ERROR("The input Matrix must be the same size as this Matrix");

    for (size_type i=0; i<nrows*ncols; i++)
    	mat[i] -= m.mat[i];

    return *this;
}

//..............................................................................
// Matrix subtraction. The the size of Matrix m must equal the size of this
// Matrix.

Matrix Matrix::operator-(const Matrix &m) const
{
    Matrix m1 = *this;
    if (nrows != m.nrows || ncols != m.ncols)
    	FATAL_ERROR("The input Matrix must be the same size as this Matrix");
    m1 -= m;
    return m1;    
}


//..............................................................................
// Matrix addition. The the size of Matrix m must equal the size of this Matrix

Matrix& Matrix::operator+=(const Matrix &m)
{
    if (nrows != m.nrows || ncols != m.ncols)
    	FATAL_ERROR("The input Matrix must be the same size as this Matrix");
    for (size_type i=0; i<nrows*ncols; i++)
    	mat[i] += m.mat[i];
    return *this;
}


//..............................................................................
// Matrix addition. The the size of Matrix m must equal the size of this Matrix

Matrix Matrix::operator+(const Matrix &m) const
{
    Matrix m1 = *this;
    if (nrows != m.nrows || ncols != m.ncols)
    	FATAL_ERROR("The input Matrix must be the same size as this Matrix");
    m1 += m;
    return m1;    
}


//.............................................................................
// Matrix multiplication Transpose(m1) * m2

Matrix Matrix::operator&(const Matrix& m) const
{
    if (m.nrows != nrows)
    	FATAL_ERROR("The input Matrix must be the same size as this Matrix");
    
    Matrix m1(ncols,m.ncols);
    Vector temp(nrows);
    size_type i,j,k,iarr;
    for (i=0, iarr=0; i<ncols; i++, iarr=m1.ncols*i)
    {
    	for (j=0; j<nrows; j++) temp[j] = mat[ncols*j + i];
    	for (j=0; j<m.ncols; j++)
    	    for (k=0; k<nrows; k++)    	    
    	    	// m1(i,j)+= *this(k,i) * m(k,j)
    	    	//
    	    	m1.mat[iarr + j] +=temp[k] * m.mat[m.ncols*k + j];
    }   	
    return m1;

}   // operator&

//.............................................................................
// Matrix multiplication m1 * Transpose(m2)

Matrix Matrix::operator%(const Matrix& m) const
{
    if (m.ncols != ncols)
    	FATAL_ERROR("Both this and the input Matrix must have the same number "
    	            "of columns");
    
    Matrix answer(nrows,m.nrows);
    for (size_type i=1; i<=nrows; i++)
    	for (size_type j=1; j<=m.nrows; j++)
    	    for (size_type k=1; k<=ncols; k++)
    	    	answer(i,j) += (*this)(i,k) * m(j,k);

    return answer;
}

//.............................................................................
// equality operator

bool Matrix::operator==(const Matrix& m) const
{ 

    if ( nrows != m.nrows || ncols != m.ncols )
    	return false;
    	
    unsigned int error = IsNotEqualTo(m);
    if (error == 0) 
    	return true;
    else
    	return false;
}

//.............................................................................
// write out a matrix eg cout << mat;

ostream& operator<<(ostream &os, const Matrix &m)
{
	os.setf(ios::showpoint);
	os.setf(ios::fixed, ios::floatfield);
	for (Matrix::size_type i=1; i<=m.nrows; i++)
	{
		os << '\n';
		for (Matrix::size_type j=1; j<=m.ncols; j++)
			os << m(i,j) << " ";
	}
	os << '\n';
	return os;
}       // operator <<

//.............................................................................
// Matrix adjoint. Restrictions: only works for square matrices

Matrix Matrix::Adjoint() const
{
    if (nrows != ncols)
    	FATAL_ERROR("Cannot calculate adjoint of non-square matrix");
    	
    if (nrows < 3)
    	FATAL_ERROR("Cannot calculate adjoint of matrix with dimensions less "
    	            "than 3");
    
    Matrix adj(nrows,ncols);
    for (size_type col = 1; col <= ncols; col++)
    	for (size_type row = 1; row <= nrows; row++)
    	{
    	    int sign = -1;
    	    if ((row+col)%2 == 0)
    	    	sign = 1;
    	    	
    	    Matrix minor = Minor(row,col);
    	    adj(row,col) = sign * minor.Determinant();
    	}
    adj = adj.Transpose();	     
    return adj;
         
} // Matrix::Adjoint

//.............................................................................
// Compare to two Matrix's - return zero if they are identical within the
// floating point precision used and the number of significant digits in the
// maximum difference between elements otherwise.

unsigned int Matrix::IsNotEqualTo(const Matrix& m) const
{
    
    if ( nrows != m.nrows || ncols != m.ncols )
    {
    	WARNING("The input Matrix must have the same dimensions as "
              "this Matrix");
    	return NumericLimits<value_type>::digits10();
    }

    unsigned int error, maxError=0;
    for (size_type i=0; i<size(); i++)
    {
    	error = NumericUtilities<value_type>::AreNotEqual(mat[i], m.mat[i]);
    	if (error > maxError) maxError = error;
    }
    	     
    return maxError;
}

//.............................................................................
// D S Moss May 1991
// 
// Inverts a symmetric positive definite matrix.
// Employs Cholesky decomposition in step (1).
// (1)M=LL(tr)   (2)compute L(-1) (3)M(-1)=L(-1)(tr)L(-1)
// 
// mat - a pointer to a column of pointers to the rows of a square or lower 
// triangular matrix containing the lower triangle of M and overwritten with L
// then with L(-1) and finally with the lower triangle of M(-1)
// 

void Matrix::CholeskyInverse()
{
    value_type a;
    size_type i,j,k;

    if (nrows != ncols)
    	FATAL_ERROR("This Matrix must be square");

    // Compute L where M=LL(tr) using Cholesky decomposition.
    for (i=1; i<=nrows; i++)
    	for (j=1; j<=i; j++)
    	{
	    a=(*this)(i,j);
	    for (k=1; k<j; k++)
	    	a-=(*this)(i,k)*(*this)(j,k);
	    if (i==j)
	    	if (a > NumericLimits<value_type>::min() * 100)
	    	    (*this)(i,i)=sqrt(a);
	    	else
	    	{
	    	    WARNING("A principal minor is singular\n");
	    	    (*this)(i,i)=0.0;
	    	}
	    else
	    	(*this)(i,j)=a/(*this)(j,j);
    	}

    // compute L(-1)
    for (j=1; j<=nrows; j++)
    {
    	(*this)(j,j)=1.0/(*this)(j,j);
    	for (i=j+1; i<=nrows; i++)
    	{
    	    a=0.0;
    	    for (k=j; k<i; k++)
    	    	a-=(*this)(i,k)*(*this)(k,j);
    	    (*this)(i,j)=a/(*this)(i,i);
    	}
    }

    // compute M(-1) = L(-1)(tr)L(-1)
    for (j=1; j<=nrows; j++)
    	for (i=j; i<=nrows; i++)
    	{
    	    a=0.0;
    	    for (k=i; k<=nrows; k++)
    	    	a+=(*this)(k,i)*(*this)(k,j);
    	    (*this)(i,j)=(*this)(j,i)=a;
    	}

} // Matrix::Cholinv

//.............................................................................
// Find column mean of *this

Vector Matrix::ColumnMean() const
{
    Vector mean(ncols);
    for (size_type i=0; i<nrows; i++)
    	for (size_type j=0; j<ncols; j++)
    	    mean[j] += mat[ncols*i + j];

    mean /= nrows;
    
    return mean;
}

//.............................................................................
// Matrix determinant of small matrices

Matrix::value_type Matrix::Determinant() const
{
    if (nrows != ncols)
    	FATAL_ERROR("Cannot calculate determinant of non-square matrix");
    	
    if (nrows == 1) return mat[0];
    	 
    if (nrows == 2) return mat[0] * mat[3] - mat[1] * mat[2];
    
    value_type det = 0.0;
    for (size_type i = 1; i <= nrows; i++)
    {
    	Matrix minor = Minor(i,1);
    	if (i%2 == 0)
            det -= mat[ncols * (i-1)] * minor.Determinant();
        else
            det += mat[ncols * (i-1)] * minor.Determinant();
    }
    return det;
}

//.............................................................................
// Eigenvalues & vectors of real symmetric matrix *this.
// 		This function calls Eigen and Esort
// Arguments:
// 
// eval - returned Vector of size n containing eigenvalues in decreasing
//        magnitude.
// 
// evec - returned matrix of size n^2 containing n*n matrix of eigenvectors
//        stored columnwise in the same order as eigenvalues.

void Matrix::Eigens(Vector &eval, Matrix &evec) const
{
    if (nrows != ncols) FATAL_ERROR("This Matrix must be square");
    
    if (eval.size() != nrows) 
    	FATAL_ERROR("The size of the input Vector must equal the number of rows"
    	            " in this Matrix");

    if (evec.nrows != nrows || evec.ncols != nrows)
    	FATAL_ERROR("The number of rows and columns of the input Matrix must "
    	            "equal the number of rows in this Matrix");
        
    size_type i, j, k;
    long double *a = new long double[nrows*(nrows+1)/2];
    long double *r = new long double[nrows*nrows];

    if (a == 0 || r == 0) FATAL_ERROR("Not enough memory");

    // Copy lower triangle of this Matrix to an array.
    for(i=1, k=0; i<=nrows; i++)
    	for(j=1; j<=i; j++)
	    a[k++] = (*this)(i,j);


    Eigen(1, nrows, a, r);
    Esort(1, nrows, a, r);
    for (j=0, k=0; j<nrows; j++)
    { 
    	eval[j]=a[k];
      	k+=j+2;
    }

    delete[] a;

    // Store eigen vectors column-wise in the evec Matrix.
    for(i=1, k=0; i<=nrows; i++)
    	for(j=1; j<=nrows; j++)
    	    evec(j,i) = r[k++];
    delete[] r;
}   // Matrix::Eigens

//.............................................................................
// The inverse of a matrix. N.B. This Matrix is overwritten by its inverse

void Matrix::Inverse()
{
    Vector vecL(ncols);
    Matrix matV(ncols,ncols);
    
    // Calculuate SVD
    //	 
    SVD(vecL,matV);
    	
    // Find maximum singular value.
    //
    value_type vecLmax = 0.0;
    for (size_type i=1; i<=ncols; i++)
    	if (vecL(i) > vecLmax) vecLmax = vecL(i);

    // Set very small values to zero in the matLinv Matrix to reduce rounding
    // error problems.
    //
    Matrix matLinv(ncols,ncols);  
    value_type vecLmin = vecLmax * NumericLimits<value_type>::epsilon();
    {for (size_type i=1; i<=ncols; i++)
    	if (vecL(i) < vecLmin)
      	    matLinv(i,i) = 0.0;
    	else
      	    matLinv(i,i) = 1.0/vecL(i);}

    // Replace this Matrix with its inverse
    *this = matV * matLinv % (*this); 	
}

//.............................................................................
// Find the largest difference between elements of matrices of the same 
// dimensions.

Matrix::value_type Matrix::MaxDifference(const Matrix& m) const
{
    if (m.ncols != ncols || m.nrows != nrows)
    {
    	WARNING("The input Matrix should have the same dimensions as "
              "this Matrix.");
    	return 99999.9;
    }
    value_type diff, maxdiff=0.0;
    for (size_type i=0; i<size(); i++)
    {
    	diff = fabs(mat[i] - m.mat[i]);
    	if (diff > maxdiff) maxdiff = diff;
    }
    return maxdiff;
}

//.............................................................................

Matrix Matrix::Minor(size_type i, size_type j) const
{
    if (nrows != ncols) 
    	FATAL_ERROR("Cannot calculate Minor of non-square matrix.");
    
    if ( i<1 || i>nrows || j<1 || j>ncols)
    	FATAL_ERROR("One or both indexes are out of range");
    	
    Matrix minor(nrows-1,ncols-1);
    for (size_type row = 1; row <= nrows; row++) 
    	for (size_type col = 1; col <= ncols; col++)
    	{
    	    if (row < i)
    	    {
    	    	if (col < j)
    	    	    minor(row,col) = mat[(ncols*(row-1) + col - 1)];
    	    	else if (col > j)
    	    	    minor(row,col-1) = mat[(ncols*(row-1) + col - 1)];
    	    }    	    	
    	    else if (row > i)
    	    {
    	    	if (col < j)
    	    	    minor(row-1,col) = mat[(ncols*(row-1) + col - 1)];
    	    	else if (col > j)
    	    	    minor(row-1,col-1) = mat[(ncols*(row-1) + col - 1)];
    	    }
    	}
    return minor;
}    	    	
    	        	       	
//.............................................................................
// Copy a row of this Matrix into a Vector
    	
Vector Matrix::ReadRow(const size_type& rowIndx) const
{
    Vector vec(ncols);

    if (rowIndx>nrows || rowIndx<1)
    	FATAL_ERROR("Matrix index out of range");
    
    size_type rowarr = ncols * (rowIndx-1);
    for (size_type i=0; i<ncols; i++)
    	vec[i] = mat[rowarr + i];
 
    return vec;
}

//.............................................................................
// Calculate inverse of square matrix = adjoint/determinant.

Matrix Matrix::SquareInverse() const
{

    // Work out determinant
    value_type det = Determinant();
    
    if (fabs(det) < NumericLimits<value_type>::epsilon())
    	FATAL_ERROR("Singular matrix: Inverse is not defined");
      
    // find Adjoint of this Matrix.
    Matrix invmat = Adjoint();

    // divide thru by the determinant if non-zero
    invmat /= det;
    
    return invmat;
}

//.............................................................................

Matrix Matrix::Transpose() const
{
    Matrix trans(ncols,nrows);
    for (size_type row = 1; row <= nrows; row++)
    	for (size_type col = 1; col <= ncols; col++)
    	    trans(col,row) = mat[ncols*(row-1) + col - 1];
    return trans;
}
    	        	       	
//.............................................................................
// Write a row of this Matrix from a Vector
    	
void Matrix::WriteRow(const Vector& rowVector, const size_type& rowIndx)
{

    if (rowIndx>nrows || rowIndx<1)
    	FATAL_ERROR("Matrix index out of range");

    if (rowVector.size() != ncols)
    	FATAL_ERROR("The size of the input Vector must be the same as the "
    	            "number of columns in this matrix");
    
    size_type rowarr = ncols * (rowIndx-1);
    for (size_type i=0; i<ncols; i++)
    	mat[rowarr + i] = rowVector[i];
}
    	        	       	
//.............................................................................
// Write a row of this Matrix from an array
    	
void Matrix::WriteRow(const value_type* rowArr, const size_type& rowIndx)
{

    if (rowIndx>nrows || rowIndx<1)
    	FATAL_ERROR("Matrix index out of range");

    size_type rowarr = ncols * (rowIndx-1);
    for (size_type i=0; i<ncols; i++)
    	mat[rowarr + i] = rowArr[i];
}