matrix.hpp

简单的C++计算库.其中实现了矩阵运算,大数运算.矩阵处理中使用lazy calculate,并且封装常见算法,可以将普通matlab程序用SDL改写.

#ifndef		_SDL_MATHS_MATRIX_MATRIX_HPP_
#define		_SDL_MATHS_MATRIX_MATRIX_HPP_

#define		MatrixData	DataBlock

#include	"MatrixIterator.hpp"
#include	"MatrixRef.hpp"
#include	"MatrixPolicy.hpp"
#include	"MatrixTrait.hpp"

SDL_MATHS_MATRIX_BEGIN

#ifdef		SDL_MATHS_NUMBERS_NP	
using	namespace	SDL_MATHS_NUMBERS_NP;
#endif

/*
*
*	BasicMatrix
*
*/

template<
		typename _NumberType = Float,
		typename _MatrixPolicy = FullMatrixPolicy<_NumberType>
		>
class	BasicMatrix
	:	protected _MatrixPolicy,
		protected _MatrixPolicy::MatrixData
{
public:

	typedef	typename NonConstType<BasicMatrix>::Type 		SelfType;
	typedef	_NumberType										NumberType;
	typedef	_MatrixPolicy									MatrixPolicy;
	typedef	typename MatrixPolicy::MatrixData 				MatrixData;
	typedef	typename MatrixPolicy::
					StoredElementIterator					SEIterator;
	typedef	typename MatrixPolicy::
					ConstStoredElementIterator				CSEIterator;

public:

	BasicMatrix(Int32 _Row, Int32 _Col,
			typename MatrixPolicy::MatrixDataType * pOtherData,
			UInt32 size)
		:	MatrixPolicy(_Row, _Col),
			MatrixData(pOtherData, size,
						size == 0 ? MatrixData::NOT_COPY_DATA : MatrixData::COPY_DATA,
						MatrixPolicy::AllocateSize(_Row, _Col))
	{
	}

	BasicMatrix(Int32 _Row, Int32 _Col = 1)
		:	MatrixPolicy(_Row, _Col),
			MatrixData(MatrixPolicy::AllocateSize(_Row, _Col)) 
	{
	}

	BasicMatrix(Int32 _Row, Int32 _Col, NumberType _InitValue)
		:	MatrixPolicy(_Row, _Col),
			MatrixData(MatrixPolicy::AllocateSize(_Row, _Col)) 
	{
		for (SEIterator	iter = SEBegin();
				iter.good(); ++iter)
			*iter = _InitValue;
	}

	BasicMatrixExpr<SelfType> expr()
	{
		return BasicMatrixExpr<SelfType>(*this);
	}

/*
*
*	Zero element
*
*/

	void	Zero()
	{
		for (SEIterator	iter = SEBegin();
				iter.good(); ++iter)
			*iter = 0;
	}

	void	Diag(NumberType _Value = 1)
	{
		Zero();

		for (int i = INDEX_OFFSET; i < row() + INDEX_OFFSET; ++i)
		{
			operator () (i, i) = _Value;
		}
	}

/*
*
*	Matrix Size
*
*/

	Int32	row() const
	{
		return MatrixPolicy::Row;
	}

	Int32	col() const
	{
		return MatrixPolicy::Col;
	}

/*
*
*	Element Visit
*
*/

	ConstRowRef<SelfType>	operator [] (Int32 _row) const
	{
		return ConstRowRef<SelfType>(*this, _row);
	}

	RowRef<SelfType>	operator [] (Int32 _row)
	{
		return RowRef<SelfType>(*this, _row);
	}

	NumberType	operator () (Int32 _row, Int32 _col) const
	{
		return MatrixPolicy::at(MatrixData::pData, _row, _col);
	}

	NumberType	operator () (Int32 _row) const
	{
		if (col() == 1)
		{
			return MatrixPolicy::at(MatrixData::pData, _row, INDEX_OFFSET);
		}
		else
		{
			return MatrixPolicy::at(MatrixData::pData, INDEX_OFFSET, _row);
		}
	}

	NumberType&	operator () (Int32 _row, Int32 _col)
	{
		return MatrixPolicy::at(MatrixData::pData, _row, _col);
	}

	NumberType&	operator () (Int32 _row)
	{
		if (col() == 1)
		{
			return MatrixPolicy::at(MatrixData::pData, _row, INDEX_OFFSET);
		}
		else
		{
			return MatrixPolicy::at(MatrixData::pData, INDEX_OFFSET, _row);
		}
	}

/*
*
*	Iterators
*
*/

	SEIterator	SEBegin()
	{
		return SEIterator(MatrixData::pData, row(), MatrixPolicy::SecondPara());
	}

	CSEIterator	SEBegin() const
	{
		return CSEIterator(MatrixData::pData, row(), MatrixPolicy::SecondPara());
	}

	CSEIterator	CSEBegin() const
	{
		return CSEIterator(MatrixData::pData, row(), MatrixPolicy::SecondPara());
	}

/*
*
*	SubMatrix
*
*/

	ConstSubMatrixRef<SelfType>	operator () (Int32 _RowBegin, Int32 _RowEnd,
											Int32 _ColBegin, Int32 _ColEnd) const
	{
		return ConstSubMatrixRef<SelfType>(*this, _RowBegin, _RowEnd, _ColBegin, _ColEnd);
	}

	SubMatrixRef<SelfType>	operator () (Int32 _RowBegin, Int32 _RowEnd,
										Int32 _ColBegin, Int32 _ColEnd)
	{
		return SubMatrixRef<SelfType>(*this, _RowBegin, _RowEnd, _ColBegin, _ColEnd);
	}

/*
*
*	Swap
*
*/

	void	Swap(SelfType& other)
	{
		MatrixPolicy::Swap(other);
		MatrixData::Swap(other);
	}

/*
*
*	operator =
*
*/

	struct	MatrixAssigner
	{
		MatrixAssigner(const SEIterator& _iter) :
			iter(_iter)
		{
		}

		MatrixAssigner& operator , (NumberType Value)
		{
			if (iter.good())
			{
				*iter++ = Value;
			}
			return *this;
		}

		MatrixAssigner& operator = (NumberType Value)
		{
			if (iter.good())
			{
				*iter++ = Value;
			}
			return *this;
		}

	private:

		SEIterator	iter;

	};

	MatrixAssigner	assigner()
	{
		return	MatrixAssigner(SEBegin());
	}

	template<typename _OtherType>
	SelfType&	operator = (const _OtherType& other)
	{
		RT_CONDITION(row() == other.row() && col() == other.col());
		
		for (SEIterator	iter = SEBegin(); iter.good(); ++iter)
		{
			*iter = other(iter.RowIdx(), iter.ColIdx());
		}

		return *this;
	}

	SelfType&	operator = (const SelfType& other)
	{
		if (&other != this)
		{
			RT_CONDITION(row() == other.row() && col() == other.col());
			for (SEIterator	iter = SEBegin(); iter.good(); ++iter)
			{
				*iter = other(iter.RowIdx(), iter.ColIdx());
			}
		}

		return *this;
	}

	template<typename _OtherType>
	SelfType&	SafeAssign(const _OtherType& other)
	{
		RT_CONDITION(row() == other.row() && col() == other.col());

		SelfType	temp(row(), col());

		for (SEIterator	iter = temp.SEBegin(); iter.good(); ++iter)
		{
			*iter = other(iter.RowIdx(), iter.ColIdx());
		}

		Swap(temp);

		return *this;
	}

	SelfType&	SafeAssign(const SelfType& other)
	{
		if (&other != this)
		{
			RT_CONDITION(row() == other.row() && col() == other.col());

			SelfType	temp(row(), col());

			for (SEIterator	iter = temp.SEBegin(); iter.good(); ++iter)
			{
				*iter = other(iter.RowIdx(), iter.ColIdx());
			}

			Swap(temp);
		}

		return *this;
	}

	template<typename _OtherType>
	SelfType&	ChangeTo(const _OtherType& other)
	{
		SelfType	new_mat(other.row(), other.col());

		new_mat = other;
		Swap(new_mat);

		return *this;
	}

	SelfType&	ChangeTo(const SelfType& other)
	{
		if (&other != this)
		{
			SelfType	new_mat(other.row(), other.col());

			new_mat = other;
			Swap(new_mat);
		}

		return *this;
	}

/*
*
*	Operators
*
*/

	template<typename _OtherType>
	Bool	operator == (const _OtherType& other) const
	{
		if (row() != other.row() || col() != other.col())
		{
			return false;
		}

		for (int i = INDEX_OFFSET; i < row()+INDEX_OFFSET; ++i)
		for (int j = INDEX_OFFSET; j < col()+INDEX_OFFSET; ++j)
		{
			if (operator () (i, j) != (NumberType)other(i, j))
			{
				return false;
			}
		}

		return true;
	}

	Bool	operator == (const SelfType& other) const
	{
		if (&other == this)
		{
			return true;
		}

		if (row() != other.row() || col() != other.col())
		{
			return false;
		}

		if (MatrixData::pData == other.MatrixData::pData)
		{
			return true;
		}

		CSEIterator		self_iter = CSEBegin();
		CSEIterator		other_iter = other.CSEBegin();

		for (; self_iter.good() && other_iter.good(); ++self_iter, ++other_iter)
		{
			if (*self_iter != *other_iter)
			{
				return false;
			}
		}

		return true;
	}

	template<typename _OtherType>
	Bool	operator != (const _OtherType& other) const
	{
		return ! operator == (other);
	}

	Bool	operator != (const SelfType& other) const
	{
		return ! operator == (other);;
	}

	SelfType&	operator ^= (Int32 Value)
	{
		RT_CONDITION(row() == col());

		SelfType	temp(row(), col());
		temp = *this ^ Value;
		Swap(temp);

		return *this;
	}

};

template<typename _NumberType>
class	BasicMatrix<_NumberType, FullMatrixPolicy<_NumberType> > 
	:	protected FullMatrixPolicy<_NumberType>,
		protected FullMatrixPolicy<_NumberType>::MatrixData
{
public:
	
	typedef	typename NonConstType<BasicMatrix>::Type 		SelfType;
	typedef	_NumberType										NumberType;
	typedef	FullMatrixPolicy<_NumberType>					MatrixPolicy;
	typedef	typename MatrixPolicy::MatrixData 				MatrixData;
	typedef	typename MatrixPolicy::
					StoredElementIterator					SEIterator;
	typedef	typename MatrixPolicy::
					ConstStoredElementIterator				CSEIterator;

public:

	BasicMatrix(Int32 _Row, Int32 _Col,
			typename MatrixPolicy::MatrixDataType * pOtherData,
			UInt32 size)
		:	MatrixPolicy(_Row, _Col),
			MatrixData(pOtherData, size, size, MatrixPolicy::AllocateSize(_Row, _Col))
	{
	}

	BasicMatrix(Int32 _Row, Int32 _Col = 1)
		:	MatrixPolicy(_Row, _Col),
			MatrixData(MatrixPolicy::AllocateSize(_Row, _Col)) 
	{
	}

	BasicMatrix(Int32 _Row, Int32 _Col, NumberType _InitValue)
		:	MatrixPolicy(_Row, _Col),
			MatrixData(MatrixPolicy::AllocateSize(_Row, _Col)) 
	{
		for (int i = 0; i < row(); ++i)
		for (int j = 0; j < col(); ++j)
		{
			MatrixData::pData[MatrixPolicy::Col * i + j] = _InitValue;
		}
	}

	BasicMatrixExpr<SelfType> expr()
	{
		return BasicMatrixExpr<SelfType>(*this);
	}

/*
*
*	Zero element
*
*/

	void	Zero()
	{
		for (int i = 0; i < row(); ++i)
		for (int j = 0; j < col(); ++j)
		{
			MatrixData::pData[MatrixPolicy::Col * i + j] = 0;
		}
	}

	void	Diag(NumberType _Value = 1)
	{
		Zero();