matrix.hpp

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

		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
	{
		INDEX_ADJUST(_row);
		INDEX_ADJUST(_col);

#ifdef	INDEX_CHECK
		RT_CONDITION(_row >= 0 && _row < MatrixPolicy::Row && _col >= 0 && _col < MatrixPolicy::Col);
#endif

		return MatrixData::pData[MatrixPolicy::Col * _row + _col];
	}

	NumberType	operator () (Int32 _row) const
	{
		INDEX_ADJUST(_row);

		if (col() == 1)
		{
			RT_CONDITION(_row >= 0 && _row < MatrixPolicy::Row);

			return MatrixData::pData[MatrixPolicy::Col * _row];
		}
		else
		{
			RT_CONDITION(_row >= 0 && _row < MatrixPolicy::Col);

			return MatrixData::pData[_row];
		}
	}

	NumberType&	operator () (Int32 _row, Int32 _col)
	{
		INDEX_ADJUST(_row);
		INDEX_ADJUST(_col);

#ifdef	INDEX_CHECK
		RT_CONDITION(_row >= 0 && _row < MatrixPolicy::Row && _col >= 0 && _col < MatrixPolicy::Col);
#endif

		return MatrixData::pData[MatrixPolicy::Col * _row + _col];
	}

	NumberType&	operator () (Int32 _row)
	{
		INDEX_ADJUST(_row);

		if (col() == 1)
		{
			RT_CONDITION(_row >= 0 && _row < MatrixPolicy::Row);

			return MatrixData::pData[MatrixPolicy::Col * _row];
		}
		else
		{
			RT_CONDITION(_row >= 0 && _row < MatrixPolicy::Col);

			return MatrixData::pData[_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);
	}

	void	SwapRow(Int32 _first, Int32 _second)
	{
		INDEX_ADJUST(_first);
		INDEX_ADJUST(_second);

#ifdef	INDEX_CHECK
		RT_CONDITION(_first>= 0 && _first < MatrixPolicy::Row && _second >= 0 && _second < MatrixPolicy::Row);
#endif
		
		if (_first == _second) return;

		Int32	first_offset = MatrixPolicy::Col * _first;
		Int32	second_offset = MatrixPolicy::Col * _second;

		for (int i = 0; i < col(); ++i)
		{
			std::swap(MatrixData::pData[first_offset + i], MatrixData::pData[second_offset + i]);
		}
	}

	void	SwapCol(Int32 _first, Int32 _second)
	{
		INDEX_ADJUST(_first);
		INDEX_ADJUST(_second);

#ifdef	INDEX_CHECK
		RT_CONDITION(_first>= 0 && _first < MatrixPolicy::Col && _second >= 0 && _second < MatrixPolicy::Col);
#endif

		if (_first == _second) return;

		Int32	last = row() * col();

		for (int i = 0; i < last; i += MatrixPolicy::Col)
		{
			std::swap(MatrixData::pData[i + _first], MatrixData::pData[i + _second]);
		}
	}

/*
*
*	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 (int i = INDEX_OFFSET; i < row()+INDEX_OFFSET; ++i)
		for (int j = INDEX_OFFSET; j < col()+INDEX_OFFSET; ++j)
		{
			operator () (i, j) = (NumberType)other(i, j);
		}

		return *this;
	}

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

			for (int i = INDEX_OFFSET; i < row()+INDEX_OFFSET; ++i)
			for (int j = INDEX_OFFSET; j < col()+INDEX_OFFSET; ++j)
			{
				operator () (i, j) = (NumberType)other(i, j);
			}
/*
			memcpy(MatrixData::pData, other.MatrixData::pData,
				MatrixPolicy::Row * MatrixPolicy::Col * sizeof(NumberType));
*/		}

		return *this;
	}

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

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

		Swap(temp);

		return *this;
	}

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

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

			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;
		}

		/*
				Candidate method : use memcmp
		*/

		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);;
	}

	template<typename _OtherType>
	SelfType&	operator += (const _OtherType& other)
	{
		*this = *this + other;
		return *this;
	}

	template<typename _OtherType>
	SelfType&	operator -= (const _OtherType& other)
	{
		*this = *this - other;
		return *this;
	}

	template<typename _OtherType>
	void	Assign(SelfType& self, const _OtherType& other, Int2Type<1>)
	{
		SelfType	temp(self.row(), self.col());

		temp = *this * other;

		self.Swap(temp);
	}

	template<typename _OtherType>
	void Assign(SelfType& self, const _OtherType& other, Int2Type<0>)
	{
		self = self * other;
	}

	template<typename _OtherType>
	SelfType&	operator *= (const _OtherType& other)
	{
		Assign(*this, other, Int2Type<IsMatConcept<_OtherType>::Result >());
		return *this;
	}
	
	template<typename _OtherType>
	SelfType&	operator /= (const _OtherType& other)
	{
		*this = *this / other;
		return *this;
	}

	SelfType&	operator /= (const SelfType& other)
	{
		if (&other == this)
		{
			Diag();
		}
		else
		{
			*this = *this / other;
		}
		return *this;
	}

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

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

		return *this;
	}

};

template<typename _NumberType = Float>
struct	matrix
{
	typedef	_NumberType	NumberType;
	typedef	BasicMatrix<NumberType, FullMatrixPolicy<NumberType> >		FullMatrix;
	typedef	BasicMatrix<NumberType, DiagMatrixPolicy<NumberType> >		DiagMatrix;
	typedef	BasicMatrix<NumberType, KDiagMatrixPolicy<NumberType> >		KDiagMatrix;
	typedef	BasicMatrix<NumberType, UpperMatrixPolicy<NumberType> >		UpperMatrix;
	typedef	BasicMatrix<NumberType, LowerMatrixPolicy<NumberType> >		LowerMatrix;
	typedef	BasicMatrix<NumberType, SymMatrixPolicy<NumberType> >		SymMatrix;
	typedef	BasicMatrix<NumberType, SparseMatrixPolicy<NumberType> >	SparseMatrix;
};

#define		FullMatrix matrix<Float>::FullMatrix
#define		DiagMatrix matrix<Float>::DiagMatrix
#define		KDiagMatrix matrix<Float>::KDiagMatrix
#define		UpperMatrix matrix<Float>::UpperMatrix
#define		LowerMatrix matrix<Float>::LowerMatrix
#define		SymMatrix matrix<Float>::SymMatrix
#define		SparseMatrix matrix<Float>::SparseMatrix

SDL_MATHS_MATRIX_END

#endif

#include	"MatrixExpr.hpp"
#include	"MatrixOpe.hpp"
#include	"MatrixFun.hpp"
#include	"MatrixAlgo.hpp"