matrix.h

虽然网上有一些现成的类

			return;
		}

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

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

		// Self scalar multiplication by P.Li(Y) [1/21/2005]
		void operator*=(const T& value)
		{
			transform(_data.begin(), _data.end(), _data.begin(), bind2nd(std::multiplies<T>(), value));
			return;
		}

		// Matrix multiplication by P.Li(Y) [1/17/2005]
		Matrix<T> operator*(const Matrix<T>& rhs) const throw(MatrixException)
		{
			if (_col != rhs.RowNum()) {
				throw MatrixException("EXECPTION: Dimensions don't agree");
			}
			Matrix<T> *pThis = const_cast<Matrix<T>*> (this);
			Matrix<T> *pRhs = const_cast<Matrix<T>*> (&rhs);
			unsigned row = _row;
			unsigned col = rhs.ColNum();
			T tmp;
			Matrix<T> a(row, col);
			for (unsigned i = 0; i < row; i++) {
				for (unsigned j = 0; j < col; j++) {
					tmp = static_cast<T>(0);
					for (unsigned k = 0; k < _col; k++) {
						tmp += pThis->At(i, k) * pRhs->At(k, j);
					}
					a.At(i, j) = tmp;
				}
			}
			return a;
		}

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

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


		//////////////////////////////////////////////////////////////////////////
		// What category? by P.Li(Y) [1/17/2005]

		// Trace by P.Li(Y) [1/20/2005]
		T Tr() throw(MatrixException)
		{
			if (_traceCalculated) {
				return _trace;
			}
			if (_row != _col) {
				throw MatrixException("EXECPTION: Not a square matrix");
			}
			_trace = 0;
			for (unsigned i = 0; i < _row; i++) {
				_trace += At(i, i);
			}
			_traceCalculated = true;
			return _trace;
		}

		// Determinant by P.Li(Y) [1/17/2005]
		T Det() throw(MatrixException)
		{
			if (_detCalculated) {
				return _det;
			}
			if (_row != _col) {
				throw MatrixException("EXECPTION: Not a square matrix");
			}
			Matrix<double> a(_row, _col);
			transform(_data.begin(), _data.end(), a._data.begin(), type_cast<T, double>());
			double det = 1.0;
			for (unsigned i = 0; i < _row; i++) {
				unsigned j;
				for (j = i; j < _col && fabs(a.At(i, j)) < _epsilon; j++);
				if (j == _col) {
					return static_cast<T>(0);
				}
				if (j != i) {
					det *= -1;
					a.SwapTwoColumns(i, j);
				}
				det *= a.At(i, i);
				double ratio;
				for (j = i + 1; j < _row; j++) {
					ratio = a.At(j, i) / a.At(i, i);
					a.At(j, i) = 0.0;
					for (unsigned k = i + 1; k < _col; k++) {
						a.At(j, k) -= a.At(i, k) * ratio;
					}
				}
			}
			_detCalculated = true;
			const char* dataType = typeid(T).name();
			// 不能四舍五入真麻烦 by P.Li(Y) [1/25/2005]
			if (strcmp(dataType, "double") != 0 && strcmp(dataType, "float") != 0) {
				if (det > 0) {
					det += 0.5;
				}
				if (det < 0) {
					det -= 0.5;
				}
			}
			return _det = static_cast<T>(det); // CAUTION: maybe not safe, by P.Li(Y) [1/24/2005]
		}

		// Cofactor of the element (i, j) by P.Li(Y) [1/24/2005]
		T Cofactor(unsigned row, unsigned col) throw(MatrixException)
		{
			if (_row != _col) {
				throw MatrixException("EXECPTION: Not a square matrix");
			}
			if (row >= _row || col >= _col) {
				throw MatrixException("EXECPTION: Index out of boundary");
			}
			Matrix<T> a(_row - 1, _col - 1);
			int idx = 0;
			for (unsigned i = 0; i < _row; i++) {
				if (i == row) {
					continue;
				}
				for (unsigned j = 0; j < _col; j++) {
					if (j == col) {
						continue;
					}
					a._data[idx++] = At(i, j);
				}
			}
			return (row + col) % 2 == 0 ? a.Det() : -1 * a.Det();
		}

		Matrix<T> RowCofactor(unsigned row) throw(MatrixException)
		{
			Matrix<T> c(1, _col);
			for (unsigned i = 0; i < _col; ++i) {
				c[0][i] = Cofactor(row, i);
			}
			return c;
		}

		// Matrix inverse, only the matrix of type float and double can be inversed,
		// otherwise an exception is thrown, by P.Li(Y) [1/20/2005]
		Matrix<T> Inv() throw(MatrixException)
		{
			const char* dataType = typeid(T).name();
			if (strcmp(dataType, "float") != 0 && strcmp(dataType, "double") != 0) {
				throw MatrixException("EXECPTION: Only the matrix of type float and double can be inversed");
			}
			if (_row != _col) {
				throw MatrixException("EXECPTION: Not a square matrix");
			}
			if (Det() == 0) {
				throw MatrixException("EXCEPTION: Singular matrix cannot be inversed");
			}

			vector<T> dataBackup = _data;
			Matrix<T> inverse = Eye(_row);
			
			T pivot;
			T pivotCandidate;
			T ratio;
			unsigned pivotRow;
			for (unsigned col = 0; col < _col; col++) {
				pivot = At(col, col);
				pivotRow = col;
				for (unsigned row = col + 1; row < _row; row++) {
					pivotCandidate = At(row, col);
					if (std::fabs(pivotCandidate) > fabs(pivot)) {
						pivot = pivotCandidate;
						pivotRow = row;
					}
				}
				// 要不要写一个初等行变换的函数？ by P.Li(Y) [1/20/2005]
				SwapTwoRows(col, pivotRow);
				inverse.SwapTwoRows(col, pivotRow);
				transform(
					_data.begin() + col * _col, 
					_data.begin() + (col + 1) * _col, 
					_data.begin() + col * _col,
					bind2nd(std::divides<T>(), pivot)
					);
				transform(
					inverse._data.begin() + col * _col, 
					inverse._data.begin() + (col + 1) * _col, 
					inverse._data.begin() + col * _col,
					bind2nd(std::divides<T>(), pivot)
					);
				for (unsigned row = 0; row < _row; row++) {
					if (row == col) {
						continue;
					}
					ratio = At(row, col);
					for (unsigned j = 0; j < _col; j++) {
						At(row, j) -= ratio * At(col, j);
						inverse.At(row, j) -= ratio * inverse.At(col, j);
					}
				}
			}
			_data = dataBackup;
			return inverse;
		}

		// In-place inverse, by P.Li(Y) [1/21/2005]
		void InvInPlace() throw(MatrixException)
		{
			try {
				*this = Inv();
			} catch (MatrixException&) {
				throw;
			}
			return;
		}


		Matrix<T> Transpose()
		{
			Matrix<T> a(_col, _row);
			for (unsigned i = 0; i < _col; i++) {
				for (unsigned j = 0; j < _row; j++ ) {
					a.At(i, j) = At(j, i);
				}
			}
			return a;
		}

		// In-place transpose by P.Li(Y) [1/21/2005]
		void TransposeInPlace()
		{
			*this = Transpose();
		}



		//////////////////////////////////////////////////////////////////////////
		// Change content, by P.Li(Y) [1/18/2005]

		// Fill the matrix with value by P.Li(Y) [1/21/2005]
		void Fill(T value)
		{
			fill(_data.begin(), _data.end(), value);
		}

		void SwapTwoRows(unsigned row1, unsigned row2) throw(MatrixException)
		{
			if (row1 >= _row || row2 >= _row) {
				throw MatrixException("EXECPTION: Index out of boundary");
			}
			swap_ranges(_data.begin() + (row1 * _col), _data.begin() + (row1 * _col + _col), _data.begin() + (row2 * _col));
		}

		void SwapTwoColumns(unsigned col1, unsigned col2) throw(MatrixException)
		{
			if (col1 >= _col || col2 >= _col) {
				throw MatrixException("EXECPTION: Index out of boundary");
			}
			T temp;
			for (unsigned i = 0; i < _row; i++) {
				temp = At(i, col1);
				At(i, col1) = At(i, col2);
				At(i, col2) = temp;
			}
		}


	private:
		static bool AcceptType()
		{
			const char* dataType = typeid(T).name();
			for (int i = 0; i < _supportedTypeNum; i++) {
				if (strcmp(dataType, _supportedTypes[i]) == 0)
					return true;
			}
			return false;
		}

	public:
		vector<T> _data;
	private:
		static const char* _supportedTypes[];
		static const int _supportedTypeNum;
		static const double _epsilon;
		unsigned _row;
		unsigned _col;
		T _trace;
		bool _traceCalculated;
		T _det;
		bool _detCalculated;
	};


	template<typename T>
		const char* Matrix<T>::_supportedTypes[] = {
			"short", "int", "long", "float", "double"
		};

	template<typename T>
		const int Matrix<T>::_supportedTypeNum = 5;

	template<typename T>
		const double Matrix<T>::_epsilon = 1e-300;

} // namespace pli

#endif // !defined(AFX_MATRIX_H__B4444097_816B_46D9_83D2_C70317AF03B3__INCLUDED_)