/*********************************************************************
 * This file is part of the cpplibs suite.
 *
 * Copyright (C) 2001-2002 Topi Mäenpää
 * Copyright (C) 2001 Jaakko Viertola
 * All rights reserved.
 *
 * This program is free software. You can redistribute and/or modify
 * it under the terms of the free software licence found in the
 * accompanying file "COPYING". The licence terms must always be
 * redistributed with this source file. The above copyright notice
 * must be reproduced in all modified and unmodified copies of this
 * source file.
 *
 * $Revision: 1.26 $
 *********************************************************************/

#ifndef _MATRIX_H
#define _MATRIX_H

#include "Exception.h"
#include "List.h"
#include "xml/XMLParser.h"
#include "Util.h"
#include "SmartPtr.h"
#include <iostream>
#include <strings.h>

/**
 * A handy macro for looping through all matrix items.
 * @param r the row index loop variable
 * @param c the column index loop variable
 * @param m the matrix to be scanned
 **/
#define AllItems(r,c,m) for(r=0;r<(m).getRows();r++) for (c=0;c<(m).getColumns();c++)

namespace util
{
	/**
	 * MatrixException is thrown when dimensions do not match or other
	 * errors occur when handling matrices.
	 **/
	class MatrixException : public Exception
	{
	public:
		MatrixException(std::string message) : Exception(message) {}
	};

	/**
	 * Matrix is used to present a two-dimensionan array of any data type.
	 * The Matrix class performs the usual matrix operations including
	 * addition, subtraction, multiplication, inversion, transpose etc.
	 **/	
	template <class T> class Matrix : virtual public Object
	{
	public:
		/**
		 * Create a size-by-size matrix. If size is omitted, a 1-by-1
		 * matrix will be created. Note that the initial contents of a
		 * a matrix are always set to zero. This is made by using the
		 * = operator with an integer argument of zero. For this reason,
		 * the = operator for an integer argument must be defined for the
		 * content class of the matrix. Note also that if a matrix is to
		 * hold objects of any class type, and if matrix atrithmetic is
		 * performed on the matrix, then the =, +, +=, -, -=, *, *=,
		 * / and /= operators of the content class must be defined for
		 * arguments of the same class type.<br>
		 * Example:<br>
		 * <pre>
		 * class MyClass
		 * {
		 *   ...
		 *
		 * public:
		 *   MyClass& operator = (const MyClass& o1);
		 *   MyClass& operator = (int value);
		 *   void operator += (const MyClass& o1);
		 *   void operator -= (const MyClass& o1);
		 *   void operator *= (const MyClass& o1);
		 *   void operator /= (const MyClass& o1);
		 *   friend MyClass operator + (const MyClass& o1, const MyClass& o2);
		 *   friend MyClass operator - (const MyClass& o1, const MyClass& o2);
		 *   friend MyClass operator * (const MyClass& o1, const MyClass& o2);
		 *   friend MyClass operator / (const MyClass& o1, const MyClass& o2);
		 *
		 *   ...
		 *  };
		 *
		 * ...
		 *
		 * Matrix&lt;MyClass&gt; mat; //invokes operator = (int) with a zero argument
		 * MyClass mc(123);
		 * mat(0,0) = mc;
		 * mat *= mc;
		 * </pre>
		 *
		 * @param size the size of the matrix
		 * @param clear if true, all matrix entries are initially set to
		 * zero. Otherwise the contents will be undefined.
		 **/
		Matrix(int size=1, bool clear=true);
		/**
		 * Create a rows-by-columns size matrix.
		 * @param rows the number of rows
		 * @param columns the number of columns
		 * @param clear if true, all matrix entries are initially set to
		 * zero. Otherwise the contents will be undefined.
		 **/
		Matrix(int rows, int columns, bool clear=true);
		/**
		 * Create a matrix that is an exact copy of another matrix.
		 * (Copy constructor)
		 * @param mat the matrix to be cloned
		 **/
		Matrix(const Matrix& other);
		/**
		 * Create a matrix that is holds the contents of another matrix
		 * typecasted to the type of this matrix. If the matrix to be cloned
		 * contains items of any other than the elementary data types, then
		 * typecast operator for the wanted type must be defined.<br>
		 * Example:<br>
		 * <pre>
		 * Matrix&lt;int&gt; a(2,2);
		 * a = 1;
		 * Matrix&lt;float&gt; b(a);
		 * </pre><br>
		 * This will result in b(0,0) holding 1.0f.
		 * @param other the matrix whose contents are to be typecasted and
		 *              copied to this matrix
		 **/
		template <class U> Matrix(const Matrix<U>& other);

		/**
		 * Create a matrix that reads its data from the given array. The
		 * size of the resulting matrix will be rows x columns. Among
		 * other things, this constructor can be used in efficiently
		 * transforming a List into a Matrix:
		 * <pre>
		 * List&lt;int&gt; lst(5);
		 * lst.addElements(5,1,2,3,4,5);
		 * Matrix&lt;int&gt; rowVector(1, 5, lst.releaseData(), true);
		 * </pre>
		 *
		 * Note that in the example, <i>lst</i> will be empty after the
		 * releaseData() call. Once can create a column matrix by
		 * exchanging the number of rows and columns.
		 *
		 * @param rows the number of rows
		 * @param columns the number of columns
		 * @param T the data array
		 * @param release if true, matrix will take the ownership of the
		 * pointer. If false, the contents will be copied.
		 **/
		Matrix(int rows, int columns, T* data, bool release=false);

		/**
		 * Create a matrix with the given number of rows and columns.
		 * Matrix contents are given as a variable-length parameter list
		 * in horizontal raster-scan order. Pay attention that the
		 * elements really are of the type you think them to be.
		 *
		 * @param rows the number of rows
		 * @param columns the number of columns
		 * @param firstElement the first element in matrix
		 * @param ... rest of the matrix data
		 *
		 * @see List::addElements(int,...)
		 **/
		Matrix(int rows, int columns, T firstElement, ...);

		/**
		 * Create a row or column matrix. Take the data for the single
		 * column from <i>vector</i>.
		 *
		 * @param vector row/column data
		 * @param row if true, create a row matrix. Otherwise create a
		 * column matrix.
		 **/
		explicit Matrix(const List<T>& vector, bool row=true);

		virtual ~Matrix();


		/**
		 * Get the data array for this matrix. The array holds the contens
		 * of the matrix from its upper left corner to the lower right
		 * one. The data is stored in a horizontal raster-scan manner.
		 * @return the data array
		 **/
		T* getData() { return _pData; }
		/**
		 * Get the data array. The actual contents of the matrix are stored in
		 * a one-dimensional array of size rows*columns. The data is ordered
		 * so that rows can be read sequentially one after another.
		 * @return a pointer to the data array
		 **/
		const T* getData() const { return _pData; }

		/**
		 * Take the ownership of the data in the matrix. The size of the
		 * matrix will be set to zero, and the internal pointer to NULL. A
		 * pointer to the data array is returned.
		 **/
		T* releaseData()
		{
			_iRows = _iColumns = 0;
			T* tmp = _pData;
			_pData = NULL;
			return tmp;
		}
		
		/**
		 * Set the size of this matrix. The initial contents of the newly
		 * allocated data array are undefined.
		 * @param rows the number of rows
		 * @param columns the number of columns
		 **/
		void setSize(int rows, int columns);

		/**
		 * Get the number of rows in a matrix.
		 **/
		int getRows() const { return _iRows; }
		/**
		 * Get the number of colums in a matrix.
		 **/
		int getColumns() const { return _iColumns; }

		/**
		 * Get the row vector at index.
		 * @param index the row index
		 * @exception MatrixException if index exceeds matrix dimensions
		 **/
		List<T> getRow(int index) const throw (MatrixException&);

		/**
		 * Get the column vector at index.
		 * @param index the column index
		 * @exception MatrixException if index exceeds matrix dimensions
		 **/
		List<T> getColumn(int index) const throw (MatrixException&);

		/**
		 * Set the row vector to index.
		 * @param index the row index
		 * @exception MatrixException if index exceeds matrix dimensions
		 * @exception MatrixException if sizes differ
		 **/
		void setRow(int index,List<T>& row) throw (MatrixException&);

		/**
		 * Set the column vector to index.
		 * @param index the column index
		 * @exception MatrixException if index exceeds matrix dimensions
		 * @exception MatrixException if sizes differ
		 **/
		void setColumn(int index,List<T>& col) throw (MatrixException&);

		/**
		 * Remove the row from the matrix at index.
		 * @param index the row index.
		 * @exception MatrixException if index exceeds matrix dimensions
		 **/
		void removeRow(int index) throw (MatrixException&);

		/**
		 * Remove the column from the matrix at index.
		 * @param index the column index.
		 * @exception MatrixException if index exceeds matrix dimensions
		 **/
		void removeColumn(int index) throw (MatrixException&);

		/**
		 * Flip the matrix in horizontal direction (columns).
		 **/
		Matrix& fliplr();
		/**
		 * Flip the matrix in vertical direction (rows).
		 **/
		Matrix& flipud();

		/**
		 * Transpose the matrix.
		 *
		 * @return a newly created matrix that is a transpose of the current one
		 **/
		Matrix getTranspose() const;
		/**
		 * Transpose the matrix in place.
		 **/
		Matrix& transpose()
		{
			if (_iRows == 1 || _iColumns == 1)
				{
					int tmp = _iRows;
					_iRows = _iColumns;
					_iColumns = tmp;
					return *this;
				}
			return operator=(getTranspose());
		}

		/**
		 * Multiply two matrices item-wise. Return the result.
		 **/
		Matrix getDotProduct(const Matrix& other) const throw (MatrixException&)
		{
			Matrix tmp(*this);
			tmp.dotProduct(other);
			return tmp;
		}

		/**
		 * Multiply two matrices item-wise. The result is stored in this
		 * matrix.
		 **/
		Matrix& dotProduct(const Matrix& other) throw (MatrixException&);

		/**
		 * Determinant of the matrix.
		 * @return the value of determinat.
		 **/
		double determinant(void) const throw (MatrixException&);
		/**
		 * Makes the inversion for the matrix.
		 *
		 * @return The inverted matrix.
		 **/
		Matrix<double> inverse(void) const throw (MatrixException&);
		/**
		 * Takes the Cofactor from the matrix at place row,col.
		 *
		 * @return The value of cofactor.
		 **/
		double cofactor(int row,int col) const throw (MatrixException&);
		/**
		 * Tests if matrix is singular.
		 **/
		bool isSingular(void);
		/**
		 * Tests if matrix is diagonal.
		 **/
		bool isDiagonal(void);
		/**
		 * Tests if matrix is symmetric.
		 **/
		bool isSymmetric(void);
		/**
		 * Tests if matrix is skew summetric.
		 **/
		bool isSkewSymmetric(void);
		/**
		 * Makes this matrix to identitymatrix.
		 * (has ones in diagonal and zeros other places).
		 * The size of matrix will be minimum of rows or columns,
		 * of course square matrix.
		 **/
		void diagonal(void);
		
		/**
		 * Calculate a sum over all matrix entries.
		 * @return the sum of all matrix entries
		 **/
		T sum(void) const;

		/**
		 * Get the maximum value.
		 **/
		T max(void) const { T m = _pData[0]; for (int i=_iRows*_iColumns;i--;) if (_pData[i] > m) m = _pData[i]; return m; }
		
		/**
		 * Get the minimum value.
		 **/
		T min(void) const { T m = _pData[0]; for (int i=_iRows*_iColumns;i--;) if (_pData[i] < m) m = _pData[i]; return m; }

		/**
		 * Add a matrix to this matrix.
		 * @param other the matrix to be added
		 * @exception MatrixException if the dimensions of
		 *            the two matrices do not match
		 **/
		void operator+= (const Matrix& other) throw (MatrixException&);
		/**
		 * Subtract a matrix from this matrix.
		 * @param other the matrix to be subtracted
		 * @exception MatrixException if the dimensions of
		 *            the two matrices do not match
		 **/
		void operator-= (const Matrix& other) throw (MatrixException&);
		/**
		 * Multiply this matrix with another matrix.
		 * @param other the matrix to multiply this matrix with
		 * @exception MatrixException if the dimensions of
		 *            the two matrices do not match
		 **/
		void operator*= (const Matrix& other) throw (MatrixException&);

		/**
		 * Add a constant value to this matrix.
		 * @param value the value to be added
		 **/
		void operator+= (T value);
		/**
		 * Subtract a constant value to this matrix.
		 * @param value the value to be subtracted
		 **/
		void operator-= (T value);
		/**
		 * Multiply this matrix by a constant value.
		 * @param value the value this matrix is to be multiplied with
		 **/
		void operator*= (T value);
		/**
		 * Divide this matrix by a constant value.
		 * @param value the value this matrix is to be divided with
		 **/
		void operator/= (T value);

		/**
		 * Set the contents of this matrix to those of another matrix.
		 * @param other the matrix to be cloned
		 **/
		Matrix& operator= (const Matrix& other);
		/**
		 * Set the contents of this matrix to those of another matrix.
		 * @param other the matrix to be cloned
		 **/
		template <class U> Matrix& operator= (const Matrix<U>& other);
		/**
		 * Set all entries in this matrix to the given value.
		 * @param value the value of each matrix entry
		 **/
		Matrix& operator= (T value);

		/**
		 * Set the contents of this matrix to those stored in the
		 * given array. The size of the array must be equal to rows*columns,
		 * and rows must be stored sequentially one after another.
		 * @param contents the new contents of the matrix in an array
		 **/
		Matrix& operator= (const T* contents);

		/**
		 * Get an item from the matrix.
		 * @param row the row index
		 * @param column the column index
		 * @return the matrix item at (row,column)
		 **/
		T& operator() (int row, int column) { return _pData[row*_iColumns+column]; }
		/**
		 * Get an item from the matrix.
		 * @param row the row index
		 * @param column the column index
		 * @return the matrix item at (row,column)
		 **/
		T operator() (int row, int column) const { return _pData[row*_iColumns+column]; }
		/**
		 * Get a sub-matrix from this matrix. Take care that the dimensions
		 * of the matrix are not exceeded.
		 * @param row the row of the upper left column of the sub-matrix
		 * @param column the column of the upper left column of the sub-matrix
		 * @param rows the number of rows to include
		 * @param columns the number of columns to include
		 **/
		Matrix operator() (int row, int column, int rows, int columns) const;

		/**
		 * Perform a matrix addition operation.
		 * @exception MatrixException if the dimensions of
		 *            m1 and m2 do not match
		 **/
		template <class U> friend Matrix<U> operator+ (const Matrix<U>& m1, const Matrix<U>& m2) throw (MatrixException&);
		/**
		 * Perform a matrix subtraction operation.
		 * @exception MatrixException if the dimensions of
		 *            m1 and m2 do not match
		 **/
		template <class U> friend Matrix<U> operator- (const Matrix<U>& m1, const Matrix<U>& m2) throw (MatrixException&);
		/**
		 * Perform a matrix multiplication operation.
		 * @exception MatrixException if the dimensions of
		 *            m1 and m2 do not match
		 **/
		template <class U> friend Matrix<U> operator* (const Matrix<U>& m1, const Matrix<U>& m2) throw (MatrixException&);
		/**
		 * Compare two matrices.
		 **/
		template <class U> friend bool operator== (const Matrix<U>& m1, const Matrix<U>& m2);
		/**
		 * Perform a constant value addition operation.
		 **/
		template <class U> friend inline Matrix<U> operator+ (const Matrix<U>& mat, U value);
		/**
		 * Perform a constant value subtraction operation.
		 **/
		template <class U> friend inline Matrix<U> operator- (const Matrix<U>& mat, U value);
		/**
		 * Perform a constant value multiplication operation.
		 **/
		template <class U> friend inline Matrix<U> operator* (const Matrix<U>& mat, U value);
		/**
		 * Perform a constant value division operation.
		 **/
		template <class U> friend inline Matrix<U> operator/ (const Matrix<U>& mat, U value);

		/**
		 * Cast the contents of a matrix to another type. If the contents
		 * of the original matrix are not of any elementary type, then
		 * an appropriate typecast operator for the content class must
		 * be defined.
		 **/
		template <class U> operator Matrix<U>();

		/**
		 * Write a matrix to an output stream.
		 **/
		template <class U> friend std::ostream& operator<< (std::ostream& sout, const Matrix<U>& m);
		/**
		 * Read a written matrix from an input stream.
		 **/
		template <class U> friend std::istream& operator>> (std::istream& sin, Matrix<U>& m);

		/**
		 * Make pivoting starting at the given row.
		 *
		 * @param row the row to start pivoting at.
		 **/
		int pivot(int row) throw (MatrixException&);
	protected:
		/**
		 * Allocate memory for a matrix. The protected variables <i>rows</i>
		 * and <i>columns</i> must be set prior to calling this method.
		 * @param clear if true, all matrix entries are set to zero
		 **/
		inline void allocate(bool clear=true);
		/**
		 * Make a copy of an existing matrix.
		 **/
		template <class U> void copy(const Matrix<U>& other);

		/**
		 * The matrix dimensions.
		 **/
		int _iRows,_iColumns;
		/**
		 * The data array.
		 **/
		T* _pData;
	};

	typedef Matrix<unsigned char> UnsignedCharMatrix;
	typedef Matrix<char> CharMatrix;
	typedef Matrix<unsigned> UnsignedIntegerMatrix;
	typedef Matrix<short> ShortMatrix;
	typedef Matrix<int> IntegerMatrix;
	typedef Matrix<long> LongMatrix;
	typedef Matrix<float> FloatMatrix;
	typedef Matrix<double> DoubleMatrix;


	/*******************************************************************
	 *                       IMPLEMENTATION                            *
	 *******************************************************************/

	template <class T> Matrix<T>::Matrix(int size, bool clear) :
		_iRows(size), _iColumns(size), _pData(NULL)
	{
		allocate(clear);
	}

	template <class T> Matrix<T>::Matrix(int r, int c, bool clear) :
		_iRows(r), _iColumns(c), _pData(NULL)
	{
		allocate(clear);
	}

	template <class T> Matrix<T>::Matrix(const Matrix<T>& other) :
		_iRows(other._iRows), _iColumns(other._iColumns), _pData(NULL)
	{
		allocate(false);
		Util::copyArray(other._pData,_pData,_iRows*_iColumns);
	}

	template <class T> Matrix<T>::Matrix(int r, int c, T* contents, bool release) :
		_iRows(r), _iColumns(c), _pData(contents)
	{
		if (!release)
			{
				_pData = NULL;
				allocate(false);
				Util::copyArray(contents,_pData,_iRows*_iColumns);
			}
	}

	template <class T>
	template <class U> Matrix<T>::Matrix(const Matrix<U>& other) :
		_pData(NULL)
	{
		copy(other);
	}

	template <class T> Matrix<T>::Matrix(int r, int c, T firstElement, ...) :
		_iRows(r), _iColumns(c), _pData(NULL)
	{
		allocate(false);
		
		va_list argp;
		// initalize var ptr
		va_start(argp, firstElement);

		int size = r*c-1;
		T* data = _pData;
		*data = firstElement;
		data++;
		
		// repeat for each arg
		while (size--)
			*(data++) = va_arg(argp,T);

		// done with args
		va_end(argp);
	}
	
	template <class T> Matrix<T>::Matrix(const List<T>& vec, bool row) :
		_pData(NULL)
	{
		int len = vec.getLength();
		if (row)
			{
				_iRows = 1;
				_iColumns = len;
			}
		else
			{
				_iRows = len;
				_iColumns = 1;
			}

		allocate();
		Util::copyArray(vec.getData(),_pData,len);
	}
	
	template <class T> Matrix<T>::~Matrix()
	{
		delete[] _pData;
	}

	template <class T> void Matrix<T>::setSize(int r, int c)
	{
		if (_iRows != r || _iColumns != c)
			{
				_iRows = r;
				_iColumns = c;
				allocate();
			}
	}
	
	template <class T> List<T> Matrix<T>::getRow(int index) const throw (MatrixException&)
	{
		List<T> result(_iColumns);
		result.setLength(_iColumns);
		if (index < 0 || index >= _iRows)
			throw MatrixException("Matrix<T>::getRow(int): index exceeds matrix dimensions");
		for (int dindex=index*_iColumns,i=0;i<_iColumns;i++,dindex++)
			result[i] = _pData[dindex];
		return result;
	}

	template <class T> List<T> Matrix<T>::getColumn(int index) const throw (MatrixException&)
	{
		List<T> result(_iRows);
		result.setLength(_iRows);
		if (index < 0 || index >= _iColumns)
			throw MatrixException("Matrix<T>::getColumn(int): index exceeds matrix dimensions");
		for (int dindex=index,i=0;i<_iRows;i++,dindex+=_iColumns)
			result[i] = _pData[dindex];
		return result;
	}

	template <class T> void Matrix<T>::setRow(int index,List<T>& row) throw (MatrixException&)
	{
		if (index < 0 || index >= _iRows)
			throw MatrixException("Matrix<T>::setRow(int,List<T>&): index exceeds matrix dimensions");
		if(_iColumns != row.getLength())
			throw MatrixException("Matrix<T>::setRow(int,List<T>&): The length of rows differ.");		
		for (int dindex=index*_iColumns,i=0;i<_iColumns;i++,dindex++)
			_pData[dindex]=row[i];
	}
	template <class T> void Matrix<T>::setColumn(int index,List<T>& col) throw (MatrixException&)
	{
		if (index < 0 || index >= _iColumns)
			throw MatrixException("Matrix<T>::setColumn(int,List<T>&): index exceeds matrix dimensions");
		if(_iRows != col.getLength())
			throw MatrixException("Matrix<T>::setColumn(int,List<T>&): The length of columns differ.");
		for (int dindex=index,i=0;i<_iRows;i++,dindex+=_iColumns)
			_pData[dindex] = col[i];
	}
	
	template <class T> void Matrix<T>::removeRow(int index) throw (MatrixException&)
	{
		if(index < 0 || index >= _iRows)
			throw MatrixException("Matrix<T>::removeRow(int): index exceeds matrix dimensions");
		int lowLimit = index*_iColumns;
		int highLimit = lowLimit+_iColumns-1;
		T* tmp = new T[(_iRows-1)*_iColumns];
		for(int i=0,k=0;i<_iRows*_iColumns;i++)
			if(i<lowLimit || i>highLimit)
					{
						tmp[k]=_pData[i];
						k++;
					}
		delete[] _pData;
		_iRows--;
		_pData=tmp;
	}

	template <class T> void Matrix<T>::removeColumn(int index) throw (MatrixException&)
	{
		if(index < 0 || index >= _iColumns)
			throw MatrixException("Matrix<T>::removeColumn(int): index exceeds matrix dimensions");
		
		T* tmp = new T[_iRows*(_iColumns-1)];
		for(int i=0,k=0;i<_iRows*_iColumns;i++)
			{
				if(index != i)
					{
						tmp[k]=_pData[i];
						k++;
					}
				else index+=_iColumns;
			}
		delete[] _pData;
		_iColumns--;
		_pData=tmp;
	}

	template <class T> Matrix<T> Matrix<T>::getTranspose() const
	{
		Matrix<T> result(_iColumns,_iRows,false);
		if (_iColumns == 1 || _iRows == 1)
			Util::copyArray(_pData,result._pData,_iRows*_iColumns);
		else
			{
				for (int r=0,ri=0;r<_iColumns;r++)
					for (int c=0,mi=r;c<_iRows;c++,ri++,mi+=_iColumns)
						result._pData[ri] = _pData[mi];
			}
		return result;
	}
	
	template <class T> int Matrix<T>::pivot(int row) throw (MatrixException&)
	{
		int k=row;
		double data=-1.0;
		double temp=0.0;
	
		for (int i=row;i<_iRows;i++)
			if ((temp=absolute(_pData[i*_iColumns+row])) > data && temp != 0.0)
				{
					data = temp;
					k = i;
				}
		int tmpK=_iColumns*k;
		if (_pData[tmpK+row] == T(0))return -1;
		if (k != row)
			{ // swap the place of rows row and k
				int tmpRow=_iColumns*row;	
				T* tmp = new T[_iColumns];
				tmp=Util::copyArray(&_pData[tmpK],tmp,_iColumns);
				Util::copyArray(&_pData[tmpRow],&_pData[tmpK],_iColumns);
				Util::copyArray(tmp,&_pData[tmpRow],_iColumns);
				delete[] tmp;
				return k;
			}
		return 0;
	}

	template <class T> Matrix<T>& Matrix<T>::fliplr()
	{
		T tmp;
		T* p1,*p2;
		for (int i=0;i<_iRows;i++)
			{
				p1 = _pData + i*_iColumns;
				p2 = p1 + _iColumns-1;
				for (int j=_iColumns>>1;j--;p1++,p2--)
					{
						tmp = *p2;
						*p2 = *p1;
						*p1 = tmp;
					}
			}
		return *this;
	}
	
	template <class T> Matrix<T>& Matrix<T>::flipud()
	{
		T tmp;
		T *p1, *p2;
		for (int i=0;i<_iRows>>1;i++)
			{
				p1 = _pData + i*_iColumns;
				p2 = _pData + (_iRows-i-1)*_iColumns;
				for (int j=_iColumns;j--;p1++,p2++)
					{
						tmp = *p2;
						*p2 = *p1;
						*p1 = tmp;
					}
			}
		return *this;
	}

	template <class T> T Matrix<T>::sum(void) const
	{
		T sum = 0;
		for (int r=0,index=0;r<_iRows;r++)
			for (int c=0;c<_iColumns;c++,index++)
				sum += _pData[index];
		return sum;
	}

	template <class T> Matrix<T>& Matrix<T>::operator= (const Matrix<T>& other)
	{
		int size = other._iRows*other._iColumns;
		if (_iRows*_iColumns == size)
			{
				Util::copyArray(other._pData,_pData,size);

				_iRows = other._iRows;
				_iColumns = other._iColumns;
			}
		else
			{
				_iRows = other._iRows;
				_iColumns = other._iColumns;

				T* tmp = new T[size];
				Util::copyArray(other._pData,tmp,size);

				delete[] _pData;
				_pData = tmp;
			}
		return *this;
	}

	template <class T>
	template <class U> Matrix<T>& Matrix<T>::operator= (const Matrix<U>& other)
	{
		copy(other);
		return *this;
	}

	template <class T> Matrix<T>& Matrix<T>::operator= (T value)
	{
		for (int i=0;i<_iRows*_iColumns;i++)
			_pData[i] = value;
		return *this;
	}

	template <class T> Matrix<T>& Matrix<T>::operator= (const T* contents)
	{
		Util::copyArray(contents,_pData,_iRows*_iColumns);
		return *this;
	}

	template <class T>
	template <class U> void Matrix<T>::copy(const Matrix<U>& other)
	{
		_iRows = other.getRows();
		_iColumns = other.getColumns();

		const U* d = other.getData();
		T* tmp = new T[_iRows*_iColumns];
		for (int i=0;i<_iRows*_iColumns;i++)
			tmp[i] = T(d[i]);

		delete[] _pData;
		_pData = tmp;
	}

	template <class T> void Matrix<T>::allocate(bool clear)
	{
		//cerr << "Allocating ...\n";
		delete[] _pData;
		//cerr << "Allocating " << rows << "x" << columns << "=" << (rows*columns*sizeof(T)) << " bytes\n";
		_pData = new T[_iRows*_iColumns];
		//cerr << "Allocated memory\n";
		if (clear)
			for (int i=0;i<_iRows*_iColumns;i++)
				_pData[i] = 0;
		//cerr << "done.\n";
	}

	template <class T>
	template <class U> Matrix<T>::operator Matrix<U>()
	{
		Matrix<U> result(_iRows,_iColumns,false);
		U* resultData = result.getData();
		for (int i=0;i<_iRows*_iColumns;i++)
			resultData[i] = (U)_pData[i];
		return result;
	}

	template <class T> Matrix<T> Matrix<T>::operator() (int row, int column, int rws, int cls) const
	{
		Matrix<T> result(rws,cls,false);
		T* resultPtr = result._pData, *thisPtr = _pData + row*_iColumns + column;
		for (int i=rws; i--; thisPtr+=_iColumns, resultPtr+=cls)
			Util::copyArray(thisPtr, resultPtr, cls);
		return result;
	}


	template <class T> std::ostream& operator<< (std::ostream& sout, const Matrix<T>& m)
	{
		int rows=m.getRows();
		int cols=m.getColumns();
		
		sout << "<matrix rows='" << rows << "' columns='" << cols << "'>" << std::endl;
		const T* data = m.getData();
		for (int i=rows*cols;i--;data++)
			{
				Util::writeXMLItem(sout,*data);
				sout << std::endl;
			}
		sout << "</matrix>";
		return sout;
	}

#define getInt(variable,name) \
	child = element->getChildNode(std::string("matrix.") + name); \
	if (child && child->getNodeType() == Node::ATTRIBUTE_NODE) \
		variable = String::parseInt(((Attr*)child)->value);

	template <class T> std::istream& operator>> (std::istream& sin, Matrix<T>& m)
	{
		using namespace util::xml;
		
		XMLParser parser;
		sin >> std::ws;

		SmartPtr<Node> node(parser.getNextNode(sin));
		if (node.get() && node->getNodeType() == Node::ELEMENT_NODE)
			{
				Element* element = (Element*)node.get();
				if (element->tagType != Element::TAG_OPENING || element->tagName != "matrix")
					throw io::IOException("operator>> (istream&, Matrix<T>&): Expecting <matrix>.");
				const Node* child;
				m._iRows = m._iColumns = -1;
				getInt(m._iRows,"rows");
				getInt(m._iColumns,"columns");
				if (m._iRows < 0 || m._iColumns < 0)
					throw io::IOException("operator>> (istream&, Matrix<T>&): The size of a matrix must be specified.");
				sin >> std::ws;

				delete[] m._pData;
				m._pData = new T[m._iRows*m._iColumns];
				T* data = m._pData;
				for (int i=m._iRows*m._iColumns;i--;data++)
					{
						Util::readXMLItem(sin,*data);
						//cerr << "Read item " << i << endl;
						if (!sin)
							throw io::IOException("operator>> (istream&, List<T>&): Premature end of input after element " + String::toString(i) + ".");
					}

				sin >> std::ws;
				SmartPtr<Node> closing(parser.getNextNode(sin));
				if (!closing.get() || closing->getNodeType() != Node::ELEMENT_NODE ||
						!closing->getChildNode("matrix") ||
						((Element*)closing.get())->tagType != Element::TAG_CLOSING)
					throw io::IOException("operator>> (istream&, Matrix<T>&): Expecting </matrix>.");
			}
		else
			throw io::IOException("operator>> (istream&, List<T>&): Input stream does not contain a list.");

		return sin;
	}

#undef getInt

	template <class T> void Matrix<T>::operator+= (const Matrix<T>& other) throw (MatrixException&)
	{
		if (_iRows != other._iRows)
			throw MatrixException("Matrix dimensions do not match: different number of rows.");
		if (_iColumns != other._iColumns)
			throw MatrixException("Matrix dimensions do not match: different number of columns.");

		for (int i=0;i<_iRows*_iColumns;i++)
			_pData[i] += other._pData[i];
	}

	template <class T> void Matrix<T>::operator+= (T value)
	{
		for (int i=0;i<_iRows*_iColumns;i++)
			_pData[i] += value;
	}

	template <class T> void Matrix<T>::operator-= (const Matrix<T>& other) throw (MatrixException&)
	{
		if (_iRows != other._iRows)
			throw MatrixException("Matrix dimensions do not match: different number of rows.");
		if (_iColumns != other._iColumns)
			throw MatrixException("Matrix dimensions do not match: different number of columns.");

		for (int i=0;i<_iRows*_iColumns;i++)
			_pData[i] -= other._pData[i];
	}

	template <class T> void Matrix<T>::operator-= (T value)
	{
		for (int i=0;i<_iRows*_iColumns;i++)
			_pData[i] -= value;
	}

	template <class T> void Matrix<T>::operator*= (const Matrix<T>& other) throw (MatrixException&)
	{
		if (_iColumns != other._iRows)
			throw MatrixException("Matrix dimensions do not match: columns must be equal to other.rows.");
		for (int r=0, i1=0, index=0;r<_iRows;r++, i1+=_iColumns)
			for (int c=0;c<other._iColumns;c++,index++)
				for (int i=0, i2=c;i<_iColumns;i++, i2+=other._iColumns)
					_pData[index] += _pData[i1+i] * other._pData[i2];
	}		

	template <class T> void Matrix<T>::operator*= (T value)
	{
		for (int i=0;i<_iRows*_iColumns;i++)
			_pData[i] *= value;
	}

	template <class T> void Matrix<T>::operator/= (T value)
	{
		for (int i=0;i<_iRows*_iColumns;i++)
			_pData[i] /= value;
	}

	template <class T> Matrix<T> operator+ (const Matrix<T>& m1, const Matrix<T>& m2) throw (MatrixException&)
	{
		if (m1._iRows != m2._iRows)
			throw MatrixException("Matrix dimensions do not match: different number of rows.");
		if (m1._iColumns != m2._iColumns)
			throw MatrixException("Matrix dimensions do not match: different number of columns.");

		Matrix<T> result(m1);

		for (int i=0;i<m1._iRows*m1._iColumns;i++)
			result._pData[i] += m2._pData[i];
		return result;
	}

	template <class T> Matrix<T> operator+ (const Matrix<T>& mat, T value)
	{
		Matrix<T> result(mat);
		result += value;
		return result;
	}

	template <class T> Matrix<T> operator- (const Matrix<T>& m1, const Matrix<T>& m2) throw (MatrixException&)
	{
		if (m1._iRows != m2._iRows)
			throw MatrixException("Matrix dimensions do not match: different number of rows.");
		if (m1._iColumns != m2._iColumns)
			throw MatrixException("Matrix dimensions do not match: different number of columns.");

		Matrix<T> result(m1);

		for (int i=0;i<m1._iRows*m1._iColumns;i++)
			result._pData[i] -= m2._pData[i];
		return result;
	}

	template <class T> Matrix<T> operator- (const Matrix<T>& mat, T value)
	{
		Matrix<T> result(mat);
		result -= value;
		return result;
	}

	template <class T> Matrix<T> operator* (const Matrix<T>& m1, const Matrix<T>& m2) throw (MatrixException&)
	{
		Matrix<T> result(m1._iRows, m2._iColumns);
		if (m1._iColumns != m2._iRows)
			throw MatrixException("Matrix dimensions do not match: m1.columns must be equal to m2.rows.");
		for (int r=0, i1=0, index=0;r<m1._iRows;r++, i1+=m1._iColumns)
			for (int c=0;c<m2._iColumns;c++,index++)
				for (int i=0, i2=c;i<m1._iColumns;i++, i2+=m2._iColumns)
					result._pData[index] += m1._pData[i1+i] * m2._pData[i2];
		return result;
	}		

	template <class T> Matrix<T> operator* (const Matrix<T>& mat, T value)
	{
		Matrix<T> result(mat);
		result *= value;
		return result;
	}

	template <class T> Matrix<T> operator/ (const Matrix<T>& mat, T value)
	{
		Matrix<T> result(mat);
		result /= value;
		return result;
	}

	template <class T> bool operator== (const Matrix<T>& m1, const Matrix<T>& m2)
	{
		if (m1._iColumns != m2._iColumns ||
				m1._iRows != m2._iRows)
			return false;
		const T* data1 = m1.getData(), *data2 = m2.getData();
		for (int i=m1._iColumns*m2._iRows;i--;data1++, data2++)
			if (*data1 != *data2)
				return false;
		return true;
	}

	
	template <class T> double Matrix<T>::determinant(void) const throw (MatrixException&)
	{
		if(_iRows != _iColumns)
			throw MatrixException("Matrix<T>::determinant(void): Matrix must be square.");
		
		int i,j,k;
		double piv,det = 1.0;
		Matrix<double> tmp(*this);
		
		for (k=0;k<_iRows;k++)
			{
				int index = tmp.pivot(k);
				if (index == -1) return 0;
				else if (index != 0)det = -det;
				det *= tmp(k,k);
				for (i=k+1;i<_iRows;i++)
					{
						piv = tmp(i,k) / tmp(k,k);
						for (j=k+1;j<_iRows;j++)tmp(i,j) -= piv * tmp(k,j);
					}
			}

		return det;
	}

	template <class T> bool Matrix<T>::isSingular(void)
	{
		if (_iRows != _iColumns) return false;
		return (determinant() == 0);
	}
	
	template <class T> bool Matrix<T>::isDiagonal(void)
	{
		if (_iRows != _iColumns)return false;
		for (int i=0,index=0; i < _iRows; i++)
			for (int j=0; j < _iColumns; j++,index++)
				if (i != j && _pData[index] != T(0))return false;
		return true;
	}

	template <class T> bool Matrix<T>::isSymmetric(void)
	{
		if (_iRows != _iColumns)return false;
		for (int i=0,index=0; i < _iRows; i++)
      for (int j=0; j < _iColumns; j++,index++)
				if (_pData[index] != _pData[j*_iColumns+i])return false;
		return true;
	}
	template <class T> bool Matrix<T>::isSkewSymmetric(void)
	{
		if (_iRows != _iColumns)return false;
		for (int i=0,index=0; i < _iRows; i++)
      for (int j=0; j < _iColumns; j++,index++)
				if (_pData[index] != -_pData[j*_iColumns+i])return false;
   return true;
	}
	
	template <class T> double Matrix<T>::cofactor(int row,int col) const throw (MatrixException&)
	{	
		if(_iRows != _iColumns)
			throw MatrixException("Matrix<T>::cofactor(int,int): Matrix must be square.");
		if (row < 0 || row >= _iRows || col < 0 || col >= _iColumns)
			throw MatrixException("Matrix<T>::cofactor(int,int): index exceeds matrix dimensions.");
		
		Matrix<T> tmp(_iRows-1,_iColumns-1);
		
		for (int i=j=r=c=0;i<_iRows;i++)
			{
				if (i==row)continue;
				for (j=c=0;j<_iColumns;j++)
					{
						if (j==col)continue;
						tmp(r,c)=_pData[i*_iColumns+j];
						c++;
					}
				r++;
			}
		double cofactor = tmp.determinant();
		if ((row+col)%2 == 1)cofactor = -cofactor;
		
		return cofactor;
	}

	template <class T> void Matrix<T>::diagonal(void)
	{
    int size=minimum(_iRows,_iColumns);
		setSize(size,size);
    for (int i=0; i < _iRows*_iColumns; i+=(_iColumns+1))_pData[i]=1;
	}

	template <class T> Matrix<double> Matrix<T>::inverse(void) const throw (MatrixException&)
	{		
		if(_iRows != _iColumns)
			throw MatrixException("Matrix<T>::inverse(): Matrix must be square matrix.");
		
		Matrix<double> temp(*this);
		Matrix<double> result(_iRows,_iColumns);
		result.diagonal(); // make diagonal matrix

		// make the tmp row
		double* tmp = new double[_iColumns];
		double divider=0.0,value=0.0;
		for(int k=0,j=0,i=0;k<_iRows;k++)
			{
				int index = temp.pivot(k);
				if(index == -1)
					{
						delete[] tmp;
						throw MatrixException("Matrix<T>::inverse(): Inversion of a singular matrix");
					}
				if(index != 0)
					{
						int tmpIndex = _iColumns*index;
						int tmpK=_iColumns*k;	
						tmp=Util::copyArray(&result.getData()[tmpK],tmp,_iColumns);
						Util::copyArray(&result.getData()[tmpIndex],&result.getData()[tmpK],_iColumns);
						Util::copyArray(tmp,&result.getData()[tmpIndex],_iColumns);
					}
				divider = temp(k,k); // take the value from diagonal
				for(j=0;j<_iColumns; j++)
					{ // divide row k with divider
						result(k,j) /= divider;
						temp(k,j) /= divider;
					}
				for(i=0;i<_iRows;i++)
					if (i != k)
						{
							value = temp(i,k);
							for (j=0;j<_iColumns;j++)
								{
									temp(i,j) -= value * temp(k,j);
									result(i,j) -= value * result(k,j);
								}
						}
			}
		delete[] tmp;// remember to free the memory.
		return result;
	}

	template <class T> Matrix<T>& Matrix<T>::dotProduct(const Matrix& other) throw (MatrixException&)
	{
		if (other.getRows() != _iRows || other.getColumns() != _iColumns)
			throw MatrixException("Matrix::dotProduct(const Matrix&): Matix sizes differ.");

		const T* sData = other.getData();
		T* tData = _pData;
		for (int i=_iRows*_iColumns; i--; tData++, sData++)
			*tData *= *sData;
		return *this;
	}
}

#endif