matrix.h

/********************************************/
/*	matrix.h 23rd February 2005				*/
/*	(c) Danny Wilson.						*/
/*	www.danielwilson.me.uk					*/
/********************************************/

#ifndef _MATRIX_H_
#define _MATRIX_H_

#include <stdlib.h>
#include <stdio.h>
#include <vector.h>

/****************************************************************/
/*						myutils::Matrix							*/
/*																*/
/*	Matrix is a C++ style container whose memory storage is		*/
/*	designed so that elements can easily be viewed at debug		*/
/*	time in MSVC++ and to be compatible with some C code in		*/
/*	which matrices are stored as one-dimensional arrays, where	*/
/*	element (i,j) would be accessed as M[i*n+j].				*/
/*																*/
/*	Element (i,j) can be accessed in one of three ways:			*/
/*		M[i][j]				clearest syntax						*/
/*		M.element[i][j]		useful for viewing during debug		*/
/*		M.array[i*n+j]		compatible with C arrays			*/
/*																*/
/****************************************************************/

namespace myutils
{
template <typename T>
class safeArray {
public:
	T *element;
	int lo,hi;
public:
	safeArray(T *set_element, const int set_lo, const int set_hi) : element(set_element), lo(set_lo), hi(set_hi) {};
	inline T& operator[](int pos){
		if(pos<lo) error("safeArray::operator[](int pos): pos<min");
		if(pos>=hi) error("safeArray::operator[](int pos): pos>=max");
		return element[pos];
	}
};

template <typename T>
class Matrix
{
public:
	/*Preserve public access for back-compatibility*/
	T *array;
	T **element;

protected:
	int protected_nrows;
	int protected_ncols;
	int initialized;

public:
	/*Default constructor*/	Matrix()
	{
		initialized=0;
		initialize(0,0);
	}
	/*Constructor*/			Matrix(int nrows, int ncols)
	{
		initialize(nrows,ncols);
	}
	/*Constructor*/			Matrix(int nrows, int ncols, T value)
	{
		initialize(nrows,ncols);
		int i,j;
		for(i=0;i<nrows;i++)
			for(j=0;j<ncols;j++)
				element[i][j]=value;
	}
	/*Destructor*/			~Matrix()
	{
		delete[] array;
		delete[] element;
	}
	Matrix<T>& initialize(int nrows, int ncols)
	{
		int i;
		array = new T[nrows*ncols];
		if (!array) error("array allocation failure in Matrix::initialize()");

		element = new T*[nrows];
		if (!element) error("element allocation failure in Matrix::initialize()");
		for(i=0;i<nrows;i++) element[i] = &(array[i*ncols+0]);

		protected_nrows=nrows;
		protected_ncols=ncols;
		initialized=1;
		return *this;
	}
	/*All current data is lost when the Matrix is resized*/
	Matrix<T>& resize(int nrows, int ncols)
	{
		int i;
		if (!initialized) return initialize(nrows,ncols);
		if((nrows==protected_nrows)&&(ncols==protected_ncols))return *this;
		
		delete[] array;
		delete[] element;

		array = new T[nrows*ncols];
		if (!array) error("array allocation failure in Matrix::resize()");
		
		element = new T*[nrows];
		if (!element) error("element allocation failure in Matrix::resize()");
		for(i=0;i<nrows;i++) element[i] = &(array[i*ncols+0]);

		protected_nrows=nrows;
		protected_ncols=ncols;
		return *this;
	}
	int nrows(){return protected_nrows;}
	int ncols(){return protected_ncols;}
	int nrows() const {return protected_nrows;}
	int ncols() const {return protected_ncols;}
/*	void error(char* error_text)
	{
		printf("Run-time error in Matrix::");
		printf("%s%\n", error_text);
		printf("Exiting to system...\n");
		exit(13);
	}*/
	/*Copy constructor*/	Matrix(const Matrix<T> &mat)
	/*	Copy constructor for the following cases:
			Matrix mat2(mat);
			Matrix mat2=mat;
		and when Matrix is returned from a function	*/
	{
		initialize(mat.protected_nrows,mat.protected_ncols);
		int i;
		for(i=0;i<protected_nrows*protected_ncols;i++)
			array[i] = mat.array[i];
	}
	/*Assignment operator*/	Matrix<T>& operator=(const Matrix<T>& mat)
	{
		//if(this==mat)return *this;
		resize(mat.nrows(),mat.ncols());
		int i;
		for(i=0;i<protected_nrows*protected_ncols;i++)
			array[i] = mat.array[i];
		return *this;
	}
#ifdef _MYUTILS_DEBUG
	/*DEBUG Subscript operator*/inline safeArray< T > operator[](int pos){
		if(pos<0) error("Matrix::operator[](int row): row<0");
		if(pos>=protected_nrows) error("Matrix::operator[](int row): row>=nrows()");
		//return element[pos];
		return safeArray< T >(element[pos],0,protected_ncols);
	};
#else
	/*Subscript operator*/inline T* operator[](int pos){return element[pos];};
#endif

	/*Matrix multiplication*/
	Matrix<T> operator*(const Matrix<T>& mat)
	{
		if(ncols()!=mat.nrows()) error("Matrix multiplication: matrices are not conformable");
		Matrix<T> result(nrows(),mat.ncols(),0.0);
		int i,j,k;
		for(i=0;i<nrows();i++)
			for(j=0;j<mat.ncols();j++)
				for(k=0;k<ncols();k++)
					result[i][j] += element[i][k] * mat.element[k][j];
		return result;
	}
	/*apply a function to every element of the matrix*/
	Matrix<T> map(T (* f)(T))
	{
		Matrix<T> result(protected_nrows,protected_ncols);
		int i,j;
		for(i=0;i<protected_nrows;i++)
			for(j=0;j<protected_ncols;j++)
				result[i][j] = f(element[i][j]);
		return result;
	}

	/* Numerical Recipes in C++ routine for inverting a square real matrix */
	Matrix<T> invert() {
		if(protected_nrows!=protected_ncols) error("Matrix inversion: must be a symmetric matrix");
		Matrix<T> a = *this;
		Matrix<T> b(protected_nrows,protected_ncols,0);
		int i;
		for(i=0;i<protected_nrows;i++) b[i][i] = 1;

		int icol,irow,j,k,l,ll;
		double big,dum,pivinv;

		int n=a.nrows();
		int m=b.ncols();
		myutils::Vector<int> indxc(n);
		myutils::Vector<int> indxr(n);
		myutils::Vector<int> ipiv(n);
		for(j=0;j<n;j++) ipiv[j]=0;
		for(i=0;i<n;i++)
		{
			big=0.0;
			for(j=0;j<n;j++)
				if(ipiv[j]!=1)
					for(k=0;k<n;k++)
					{
						if(ipiv[k]==0)
						{
							if(fabs(a[j][k])>=big)
							{
								big=fabs(a[j][k]);
								irow=j;
								icol=k;
							}
						}
					}
			++(ipiv[icol]);
			if(irow!=icol)
			{
				for(l=0;l<n;l++) SWAP(a[irow][l],a[icol][l]);
				for(l=0;l<m;l++) SWAP(b[irow][l],b[icol][l]);
			}
			indxr[i]=irow;
			indxc[i]=icol;
			if(a[icol][icol] == 0.0) error("Matrix inversion: Singular Matrix");
			pivinv=1.0/a[icol][icol];
			a[icol][icol]=1.0;
			for(l=0;l<n;l++) a[icol][l] *= pivinv;
			for(l=0;l<m;l++) b[icol][l] *= pivinv;
			for(ll=0;ll<n;ll++)
				if (ll != icol)
				{
					dum=a[ll][icol];
					a[ll][icol]=0.0;
					for(l=0;l<n;l++) a[ll][l] -= a[icol][l]*dum;
					for(l=0;l<m;l++) b[ll][l] -= b[icol][l]*dum;
				}
		}
		for(l=n-1;l>=0;l--)
		{
			if(indxr[l]!=indxc[l])
				for(k=0;k<n;k++)
					SWAP(a[k][indxr[l]],a[k][indxc[l]]);
		}
		return a;
	}
};

};		// namespace myutils

#endif	// _MATRIX_H_