matrix.h

矩阵的运算。建立了一个矩阵运算的模板类

//#include "abnormal.h"
#include <iostream.h>
#include<stdlib.h>
#include"exception.h"

template<class T>class Matrix;
template<class T>istream& operator>>(istream& in,Matrix<T> &m);
template<class T>ostream& operator<<(ostream& out,Matrix<T>& m);


template<class T>
class Matrix 
{ 
		int rows, cols; // 矩阵维数
		T *element; // 元素数组
		friend istream& operator>>(istream& in,Matrix<T> &m);
		friend ostream& operator<<(ostream& out,Matrix<T>& m);
	 
	public:

		Matrix(int r = 0, int c = 0)
		{// 类M a t r i x的构造函数
			// 验证r和c的合法性
			if (r < 0 || c < 0) throw BadInitializers();
			if ((!r || !c) && (r || c))
				throw BadInitializers();
			// 创建矩阵
			rows = r; cols = c;
			element = new T [r * c];
		}
		Matrix(const Matrix<T>& m)//复制构造函数
		{
			rows = m.rows;
			cols = m.cols;
			element = new T[rows * cols];
			for(int i =0; i<rows*cols;i++)
				element[i] = m.element[i];
		}	

		~Matrix() {delete [] element;}

		inline int GetRows() const {return rows;}

		inline int GetCols() const {return cols;}

		inline T& operator()(int i, int j) const
		{
			if (i < 1 || i > rows || j < 1 || j > cols) 
				throw OutOfBounds();
			return element[(i - 1) * cols + j - 1];
		}
		Matrix<T>& operator=(const Matrix<T>& m)
		{
			if(this != &m)
			 {
				rows = m.rows;
				cols = m.cols;
				delete []element;
				element = new T[rows * cols];
				for(int i =0; i<rows * cols;i++)
					element[i] = m.element[i];
			}
			return *this;
		}

		Matrix<T> operator+() const // 一元加法  
		{
			Matrix<T> w(rows,cols);
			for(int i=0;i<rows*cols;i++)
				w.element[i] = +element[i];
			return w;			  
		}
	       

		Matrix<T> operator+(const Matrix<T>& m) const
		{
			if(rows!=m.rows || cols!=m.cols)
				throw SizeMismatch();
		    // 返回w = -(*this)
            // 创建结果数组w
        	Matrix<T> w(rows,cols);
          	for (int i = 0; i < rows*cols; i++)
               w.element[i] = element[i] + m.element[i] ;
           return w;
		}
		Matrix<T> operator-() const // 一元减法  相当于在对一元素求相反值
		{
			Matrix<T> w(rows,cols);
			for(int i=0;i<rows*cols;i++)
				w.element[i] = -element[i];
			return w;
		}
		Matrix<T> operator-(const Matrix<T>& m) const
		{
			if(rows!=m.rows || cols!=m.cols)
				throw SizeMismatch();
		    // 返回w = -(*this)
            // 创建结果数组w
        	Matrix<T> w(rows,cols);
          	for (int i = 0; i < rows*cols; i++)
               w.element[i] = element[i] - m.element[i] ;
           return w;
		}
		Matrix<T> operator*(const Matrix<T>& m) const
		{
			// 矩阵乘法，返回w = (*this) * m.
			if (cols != m.rows) throw SizeMismatch();
			Matrix<T> w(rows, m.cols); // 结果矩阵
			// 为*this, m和w定义游标
			// 并设定初始位置为( 1 , 1 )
			int ct = 0, cm = 0, cw = 0;
			// 对所有的i和j计算w ( i , j )
			for (int i = 1; i <= rows; i++) 
			{
				// 计算出结果的第i 行
				for (int j = 1; j <= m.cols; j++) 
				{
					// 计算w ( i , j )的第一项
					T sum = element[ct] * m.element[cm];
					// 累加其余项
					for (int k = 2; k <= cols; k++) 
					{
						ct++; // 指向*this第i行的下一个元素
						cm += m.cols; // 指向m 的第j 列的下一个项
						sum += element[ct] * m.element[cm];
					}
					w.element[cw++] = sum; // 保存w ( i , j )
					// 重新调整至行首和下一列
					ct -= cols - 1;
					cm = j;
				}
				// 重新调整至下一行的行首和第一列
				ct += cols;
				cm = 0;
			}
			return w;
		}
		Matrix<T> operator* (const T& v)const
		{
			Matrix<T> w(rows,cols);
			for(int i =0; i<rows*cols;i++)
				w.element[i] =element[i]*v;
			return w;
		}
		Matrix<T> operator/ (const T& v)const
		{
			Matrix<T> w(rows,cols);
			for(int i =0; i<rows*cols;i++)
				w.element[i] =element[i]/v;
			return w;
		}
			

		Matrix<T>& operator+=(const Matrix<T>& w)
		{
			if(rows!=w.rows || cols!=w.cols)
				throw SizeMismatch();
			for(int i=0; i<rows*cols;i++)
			    element[i] +=w.element[i];
			return *this;
		} 
		Matrix<T>& operator-=(const Matrix<T>& w)
		{
			if(rows!=w.rows || cols!=w.cols)
				throw SizeMismatch();
			for(int i=0; i<rows*cols;i++)
			    element[i] -= w.element[i];
			return *this;
		}

        Matrix& operator*= (const T& t)
		{
			for(int i =0;i<rows*cols;i++)
				element[i]*=t;
			return *this;
		}
		Matrix& operator/= (const T& t)
		{
			for(int i =0;i<rows*cols;i++)
				element[i]/=t;
			return *this;
		}

		Matrix<T> TransPose()const
		{
			Matrix<T> w(cols,rows);
			for(int i =0;i< rows;i++)
			{
				for(int j=0; j<cols;j++)
				{
					w.element[j*rows + i] = element[i* cols + j];
				}
			}
			return w;
		}
		Matrix<T> Inverse() //矩阵求逆，注意：返回的矩阵为T型
		{
			if(rows != cols || rows == 0) {
				cout<<"The Inverse Matrix doesn't exist!\n";
				exit(1);
			}
			int n = rows;
			Matrix<T> m1(*this);
			Matrix<T> m2(n,n);
			
			int i,j;
			//初始化矩阵
			for(i = 1;i<=n;i++) {
				for(j = 1;j<=n;j++)
					m2(i,j) = 0;
			}
			for(i =1 ;i<= n; i++)
				m2(i,i) = 1;
			for( i = 1;i <= n; i++){		
				if(!Inverse1(i,m1,m2)) { //第i+n行减去第i行的  倍
					cout<<"The Inverse Matrix doesn't exist!\n";
				exit(1);
				}
			}
			  
			for(i = n-1; i >= 1; i--){
				if(!Inverse2(i,m1,m2)) { //第i-n行减去第i行的  倍
					cout<<"The Inverse Matrix doesn't exist!\n";
				exit(1);
				}
			}
		   

			return m2;			
		}
		double Det()  //求矩阵的行列式
		{
			if(rows != cols || rows == 0) {
				cout<<"The Inverse Matrix doesn't exist!\n";
				exit(1);
			}
			double result = 1.0;
			int i;
			Matrix<T> m(*this);			 
			for( i = 1;i <= rows; i++){
				if(!Det(i,m,result)) { //第i+n行减去第i行的  倍
					cout<<"The Det of the Matrix is zero!\n";
				exit(1);
				}
			}
			return result;
		}
	
			

	private:
		bool Inverse1(int r ,Matrix<T>& m1 ,Matrix<T>& m2)
		{//说明：j行减去i行，j>i ,j++
			//利用矩阵求逆公式
			int n = m1.rows;
			int i,j;
			T temp;
			if( m1(r,r) != 0 ){
				temp = m1(r,r);
				for(i = 1; i<=n;i++){//将需变换的第一个元素变为1
					m1(r,i) /= temp;
					m2(r,i) /= temp;
				}
			}
			else 
			{
				if(r==n) return false;
				for( i = r+1 ;i <= n ;i++){
					if( m1(i,r) != 0){//调换两行
						for(j = 1; j <= n; j++){
							temp = m1(r,j);
							m1(r,j) = m1(i,j);
							m1(i,j) = temp;

							temp = m2(r,j);
							m2(r,j) = m2(i,j);
							m2(i,j) = temp;
						}
					}
				}
				if(i > n)return false;
				temp = m1(r,r);
				for(i = 1; i<=n;i++){//将需变换的第一个元素变为1
					m1(r,i) /= temp;
					m2(r,i) /= temp;
				}
			}

			if(r == n) return true;

			for( i = r+1; i <= n; i++){
				if( m1(i,r) == 0) continue;
				temp = m1(i,r);
				for(j = 1; j <= n; j++){
					m1(i,j) -= temp * m1(r,j);
					m2(i,j) -= temp * m2(r,j);
				}
			}
			return true;
		}
		bool Inverse2(int r, Matrix<T>& m1 ,Matrix<T>& m2)
		{//说明：j行减去i行，j>i,j--
			//利用矩阵求逆公式
			int n = m1.rows;
			int i,j;
			T temp;
			for(i = r+1; i <= n; i++){
				if( m1(r , i ) == 0)continue;
				temp = m1(r,i);
				for(j = 1; j <= n; j++){
					m1(r,j) -= temp * m1(i,j);
					m2(r,j) -= temp * m2(i,j);
				}
			}
			return true;
		}
		bool Det(int r,Matrix& m,double& result)
		{
			int n = m.rows;
			int i,j;
			T temp;
			if( m(r,r) != 0 ){
				temp = m(r,r);
				result *= (float)m(r,r);
				for(i = 1; i<=n;i++) //将需变换的第一个元素变为1
					m(r,i) /= temp;
				 
			}
			else 
			{
				if(r==n) return false;
				for( i = r+1 ;i <= n ;i++){
					if( m(i,r) != 0){//调换两行
						for(j = r; j <= n; j++){
							temp = m(r,j);
							m(r,j) = m(i,j);
							m(i,j) = temp;
						}
						result *= -1;
						break;
					}
				}
				if(i>n) return false;
				temp = m(r,r);
				result *= m(r,r);
				for(i = r;i <= n;i++)
					m(r,i) /= temp;

			}
			for(i = r+1;i <= n; i++){
				if( m(i,r) == 0) break;
				temp = m(i,r);
				for( j = r; j<= n;j++)
					m(i,j) -= temp* m(r,j);
			}
			return true;
		}

			 
			

};

template<class T>
 istream& operator>>(istream& in,Matrix<T> &m)
		{ 
			int i,j;
			for( i = 0;i<m.rows;i++){
				cout<<"请输入第 "<<i+1<<" 行元素:"<<endl;
				for(j=0; j<m.cols;j++) 
					in>>m.element[i*m.cols + j];
			}
			return in;
				 
		}
template<class T>
ostream& operator<<(ostream& out,Matrix<T>& m)
{
	int i,j;
	for(i =0; i<m.rows;i++){
		 out<<"第 "<<i+1<<"行:"<<endl;
		 for(j=0; j<m.cols;j++) 
			 out<<m.element[i*m.cols + j]<<"\t";
			 
		 out<<endl;
	}
	return out;
	
}



/*/////////////////////////////////////////////////////////////////
                    TEST CODE
//////////////////////////////////////////////////////////////////

	Matrix<float> m(3,3);
	m(1,1) = 16 ;
	m(1,2) = 5 ;
	m(1,3) = 68 ;
	m(2,1) = 135 ;
	m(2,2) = 8 ;
	m(2,3) = 17 ;
	m(3,1) = 3 ;
	m(3,2) = 19 ;
	m(3,3) = 23 ;

//	cin>>m;
	cout<<m;
	Matrix<float> mInv = m.Inverse();
	cout<<mInv<<endl;
	cout<<m.Det()<<endl;
	cout<<mInv.Det()<<endl;

	Matrix<float> m2(2,2);
	cin>>m2;
	cout<<m2; 
	cout<<m2.Inverse();

*/
////////////////////////////////////////////////////////////////////