matrix.h

虽然网上有一些现成的类

/********************************************************************
created:	2005-01-17
filename: 	Matrix.h
author:		Peng Li (Younger) @ CASIA
email:		[e-mail]

purpose:	To provide a template Matrix class with basic matrix 
			operations
caution:	Exceptions may be thrown out from some of the functions
*********************************************************************/

#if !defined(AFX_MATRIX_H__B4444097_816B_46D9_83D2_C70317AF03B3__INCLUDED_)
#define AFX_MATRIX_H__B4444097_816B_46D9_83D2_C70317AF03B3__INCLUDED_

#if _MSC_VER > 1000
#pragma once
#endif // _MSC_VER > 1000

#pragma warning( disable : 4290 )

#include <iostream>
#include <iomanip>
#include <vector>
#include <string>
#include <algorithm>
#include <functional>
#include <cmath>
using std::string;
using std::vector;

namespace pli {

	class MatrixException
	{
	public:
		MatrixException(const char* info)
		{
			cause = info;
		}
		virtual ~MatrixException()
		{
		}

		string GetCause()
		{
			return cause;
		}

		void Print(std::ostream& out)
		{
			out << cause << std::endl;
		}
	private:
		string cause;
	};

	// Function obeject class by P.Li(Y) [1/18/2005]
	template<typename T1, typename T2>
	struct type_cast : public std::unary_function<T1, T2> {
		T2 operator()(T1 x) const {
			return static_cast<T2>(x);
		}
	};

	// Forward declaration by P.Li(Y) [1/17/2005]
	template<typename T> class Matrix;

	// Print contents of the Matrix by P.Li(Y) [1/17/2005]
	template<typename T>
		std::ostream& operator<<(std::ostream& out, const Matrix<T>& m)
	{
		int nIdx = 0;
		for (unsigned i = 0; i < m._row; i++) {
			for (unsigned j = 0; j < m._col; j++) {
				out << std::setw(14) << std::right << m._data[nIdx++];
			}
			out << endl;
		}
		return out;
	}

	// Scalar plus by P.Li(Y) [1/17/2005]
	template<typename T>
		Matrix<T> operator+ (const T& value, const Matrix<T>& m)
	{
		return m + value;
	}

	// Scalar multiplication by P.Li(Y) [1/17/2005]
	template<typename T>
		Matrix<T> operator* (const T& value, const Matrix<T>& m)
	{
		return m * value;
	}



	// class definition by P.Li(Y) [1/19/2005]
	template<typename T> 
	class Matrix
	{
	public:
		//////////////////////////////////////////////////////////////////////////
		// Static fuctions by P.Li(Y) [1/18/2005]

		// Create an identity matrix with size of row by row, by P.Li(Y) [1/18/2005]
		static Matrix<T> Eye(unsigned row) throw(MatrixException)
		{
			if (!AcceptType()) {
				throw MatrixException("EXECPTION: Not supported type");
			}
			Matrix<T> a(row, row);
			for (unsigned i = 0; i < row; i++) {
				a.At(i, i) = 1;
			}
			return a;
		}

		//////////////////////////////////////////////////////////////////////////
		// Constructions and Destruction  by P.Li(Y) [1/17/2005]
		Matrix() throw(MatrixException)
		{
			Matrix(0, 0);
		}
		Matrix(unsigned row, unsigned col, T iniValue = static_cast<T>(0)) throw(MatrixException)
		{
			if (!AcceptType()) {
				throw MatrixException("EXECPTION: Not supported type");
			}
			Init(row, col);
			fill(_data.begin(), _data.end(), iniValue);
		}

		template<typename TT>
			Matrix(const Matrix<TT>& ref)
		{
			if (!AcceptType()) {
				throw MatrixException("EXECPTION: Not supported type");
			}
			Init(ref.RowNum(), ref.ColNum());
			transform(ref._data.begin(), ref._data.end(), _data.begin(), type_cast<TT, T>());
		}

		template<typename TT>
			void operator=(const Matrix<TT>& rhs)
		{
			if (!AcceptType()) {
				throw MatrixException("EXECPTION: Not supported type");
			}
			Init(rhs.RowNum(), rhs.ColNum());
			transform(rhs._data.begin(), rhs._data.end(), _data.begin(), type_cast<TT, T>());
			return;
		}


		void Init(unsigned row, unsigned col)
		{
			_row = row;
			_col = col;
			_data.resize(_row * _col);
			_traceCalculated = false;
			_detCalculated = false;
		}

		virtual ~Matrix()
		{
		}

		//////////////////////////////////////////////////////////////////////////
		// Accessments, by P.Li(Y) [1/17/2005]
		unsigned RowNum() const
		{
			return _row;
		}
		void SetRowNum(unsigned nRowNum)
		{
			_row = nRowNum;
			_data.resize(_row * _col);
			fill(_data.begin(), _data.end(), static_cast<T>(0));
		}
		unsigned ColNum() const
		{
			return _col;
		}
		void SetColNum(unsigned nColNum)
		{
			_col = nColNum;
			_data.resize(_row * _col);
			fill(_data.begin(), _data.end(), static_cast<T>(0));
		}
		void SetSize(unsigned nRowNum, unsigned nColNum)
		{
			Init(nRowNum, nColNum);
			fill(_data.begin(), _data.end(), static_cast<T>(0));
		}

		T& At(unsigned row, unsigned col) throw(MatrixException)
		{
			if (row >= _row || col >= _col) {
				throw MatrixException("EXECPTION: Index out of boundary");
			}
			_detCalculated = false;
			return _data[row * _col+ col];
		}

		Matrix<T> Row(unsigned row) throw(MatrixException)
		{
			if (row >= _row) {
				throw MatrixException("EXECPTION: Index out of boundary");
			}
			Matrix<T> rowVector(1, _col);
			copy(_data.begin() + row * _col, _data.begin() + (row + 1) * _col ,rowVector._data.begin());
			return rowVector;
		}

		Matrix<T> Col(unsigned col)
		{
			if (col >= _col) {
				throw MatrixException("EXECPTION: Index out of boundary");
			}
			Matrix<T> colVector(_row, 1);
			for (unsigned row = 0; row < _row; row++) {
				colVector._data[row] = At(row, col);
			}
			return colVector;
		}

		//////////////////////////////////////////////////////////////////////////
		// Operators

		// Print contents by P.Li(Y) [1/17/2005]
		template<typename T>
			friend std::ostream& operator<< (std::ostream& os, const Matrix<T>& m);

		T* operator[](unsigned row) // For convenient use like a[i][j], by P.Li(Y) [Mar 9, 2005]
		{
			return &_data[row * _col];
		}

		const T* operator[](unsigned row) const
		{
			return &_data[row * _col];
		}

		bool operator==(const Matrix<T>& rhs)
		{
			return (_data == rhs._data);
		}

		Matrix<T> operator-()
		{
			Matrix<T> a(*this);
			transform(a._data.begin(), a._data.end(), a._data.begin(), negate<T>());
			return a;
		}

		// Scalar plus by P.Li(Y) [1/17/2005]
		Matrix<T> operator+(const T& value) const
		{
			Matrix<T> a(*this);
			a += value;
			return a;
		}

		// In-place scalar plus, by P.Li(Y) [1/21/2005]
		void operator+=(const T& value)
		{
			transform(_data.begin(), _data.end(), _data.begin(), bind2nd(std::plus<T>(), value));
			return;
		}

		// Scalar subtraction by P.Li(Y) [1/17/2005]
		Matrix<T> operator-(const T& value) const
		{
			Matrix<T> a(*this);
			a -= value;
			return a;
		}

		// In-place scalar subtraction, by P.Li(Y) [1/21/2005]
		void operator-=(const T& value)
		{
			transform(_data.begin(), _data.end(), _data.begin(), bind2nd(std::minus<T>(), value));
			return;
		}

		// Matrix plus by P.Li(Y) [1/17/2005]
		Matrix<T> operator+(const Matrix<T>& rhs) const throw(MatrixException)
		{
			Matrix<T> a(*this);
			try {
				a += rhs;
			} catch (MatrixException& ) {
				throw;
			}
			return a;
		}

		// In-place matrix plus by P.Li(Y) [1/17/2005]
		void operator+=(const Matrix<T>& rhs) throw(MatrixException)
		{
			if (_row != rhs.RowNum() || _col != rhs.ColNum()) {
				throw MatrixException("EXECPTION: Dimensions don't agree");
			}
			transform(_data.begin(), _data.end(), rhs._data.begin(), _data.begin(), std::plus<T>());
			return;
		}

		// Matrix substraction by P.Li(Y) [1/17/2005]
		Matrix<T> operator-(const Matrix<T>& rhs) const throw(MatrixException)
		{
			Matrix<T> a(*this);
			try {
				a -= rhs;
			} catch (MatrixException&) {
				throw;
			}
			return a;
		}

		// In-place matrix substraction by P.Li(Y) [1/17/2005]
		void operator-=(const Matrix<T>& rhs) throw(MatrixException)
		{
			if (_row != rhs.RowNum() || _col != rhs.ColNum()) {
				throw MatrixException("EXECPTION: Dimensions don't agree");
			}
			transform(_data.begin(), _data.end(), rhs._data.begin(), _data.begin(), std::minus<T>());