matrix.hpp

随机数类,产生随机数

		Matrix<T> matrix(*this);
		return(matrix+(-a));
	}
	template <class T> Matrix<T> Matrix<T>::operator -(Matrix<T>& matrix)//当前矩阵减去矩阵matriz,返还一个新的矩阵
	{
		Matrix<T> temp(*this);
		if(m_row!=matrix.RowNum()||m_col!=matrix.ColNum())//如果两个矩阵的大小不一样，抛出异常，函数结束
			throw Exception("Matrix","Operator-","the two matices are not in same size.");
		temp-=matrix;
		temp.DoAdjust(); 
		//cout<<"here"<<endl;
		return temp;
	}
	template <class T> Matrix<T>& Matrix<T>::operator *=(T a)//矩阵的每一个元素乘以常数a，返回当前矩阵
	{
		unsigned long row, col;
		for(row=1;row<=m_row;row++)
			for(col=1;col<=m_col; col++)
			{
				this->SetValue(row,col,((*this)(row,col)*a)); 
			}
		DoAdjust();//误差调整
		return *this;
	}
	template <class T> Matrix<T> Matrix<T>::operator *(T a)//矩阵数乘常数a，返回一个新的矩阵
	{
		Matrix<T> matrix(*this);
		matrix*=a;
		return matrix;
	}
	template <class T> Matrix<T>& Matrix<T>::operator *=(Matrix<T>& matrix)//当前矩阵乘以矩阵matrix,返回当前矩阵
	{
		if(this->m_col!=matrix.RowNum())//如果矩阵的规模不匹配，抛出异常，程序结束
			throw Exception("Matrix","operator *=","the size of two matrices is mismatch.");
		unsigned k,j,i;
		T t;
		Matrix<T> temp(*this);//临时矩阵，当前矩阵的数值复制到这个矩阵之中
		m_col=matrix.ColNum();//当前矩阵的列数修改成乘数矩阵的列数
		delete[] this->m_buf;
		this->m_buf=new T[this->m_row*this->m_col];
		this->GenZeros();//当前矩阵的所有元素的值设置为零 
		//n=m_row;
		for(k=1;k<=m_row;k++)
			for(j=1;j<=m_col;j++)
			{
				t=0.0;
				for(i=1;i<=matrix.RowNum();i++)
					t+=temp(k,i)*matrix(i,j);
				this->SetValue(k,j,t); 
			}
		this->DoAdjust();  
		return *this;
	}
	template <class T> Matrix<T> Matrix<T>::operator *(Matrix<T>& matrix)//当前矩阵乘以矩阵matrix,返回一个新的矩阵
	{
		if(this->m_col!=matrix.RowNum())//如果矩阵的规模不匹配，抛出异常，程序结束
			throw Exception("Matrix","operator *","the size of two matrices is mismatch.");
		Matrix<T> temp(*this);
		temp*=matrix;
		return temp;
	}
	template <class T> Matrix<T>& Matrix<T>::operator /=(T a)
	{
		if(ValueAbs(a)<=this->m_Tolerance)//如果a为零，抛出异常，终止程序
			throw Exception("Matrix","operator /=","Divided by zero");
		T c=1/a;
		(*this)*=c;
		return *this;
	}
	template <class T> Matrix<T> Matrix<T>::operator /(T a)
	{
		if(ValueAbs(a)<=this->m_Tolerance)//如果a为零，抛出异常，终止程序
			throw Exception("Matrix","operator /","Divided by zero");
		Matrix<T> matrix(*this);
		matrix/=a;
		return matrix;
	}
	template <class T> Matrix<T>& Matrix<T>::operator ^=(double x)
	{
		unsigned long row, col;
		T temp1,temp2;
		for(row=1;row<=m_row;row++)
			for(col=1;col<=m_col;col++)
			{
				temp1=(*this)(row,col);
				temp2=(T)pow((double)temp1,x);
				this->SetValue(row,col,temp2); 
			}
		return *this;
	}
	template <class T> Matrix<T> Matrix<T>::operator ^(double x)
	{
		Matrix<T> matrix(*this);
		matrix^=x;
		return matrix;
	}
	template <class T> T Matrix<T>::Trace()
	{
		if(!this->IsSquare())//如果矩阵不是方阵，抛出异常，结束
			throw Exception("Matrix","Trace","the matrix is not square, can not calculate trace");
		T temp;
		unsigned long row;
		temp=(T)0.0;
		for(row=1;row<=m_row;row++)
			temp+=this->GetValue(row,row);
		return temp;
	}
	template <class T> ostream& operator <<(ostream& os,const Matrix<T>& matrix)
	{
		unsigned long row, col;
		for(row=1;row<=matrix.m_row;row++)
		{
			for(col=1;col<=matrix.m_col;col++)
			{
				//os.setf(ios::left);
				os.width(15); 
				os<<matrix.m_buf[(row-1)*matrix.m_col+col-1]<<"  ";
			}
			os<<endl;
		}
		return os;
	}
	template <class T> Matrix<T> operator +(T a, Matrix<T>& m)
	{
		Matrix<T> temp(m);
		return (temp+=a);
	}
	template <class T> Matrix<T> operator -(T a, Matrix<T>& m)
	{
		Matrix<T> temp(m);
		unsigned long row,col;
		for(row=1;row<=m.RowNum();row++)
			for(col=1;col<=m.ColNum();col++)
				temp.SetValue(row,col,(a-m(row,col)));
		return temp;
	}
	template <class T> Matrix<T> operator *(T a, Matrix<T>& m)
	{
		Matrix<T> temp(m);
		return (temp*=a);
	}
	template <class T> Matrix<T> operator /(T a, Matrix<T>& m)
	{
		unsigned long row, col;
		for(row=1;row<=m.RowNum();row++)
			for(col=1;col<=m.ColNum();col++)
			{
				if(m.ValueAbs(m(row,col))<=m.m_Tolerance)
					throw Exception("Template operator and Matrix","/","divied by zero");
			}
		Matrix<T> temp(m);
		for(row=1;row<=m.RowNum();row++)
			for(col=1;col<=m.ColNum();col++)
			{
				temp.SetValue(row,col,a/temp(row,col)); 	
			}
		return temp;
		
	}
	//function template
	template <class T> Matrix<T> Sin(Matrix<T>& m)//对所有的元素取正弦，返回新的矩阵
	{
		Matrix<T> temp(m);
		return temp.Sin();
	}
	template <class T> Matrix<T> Cos(Matrix<T>& m)//对所有元素取余弦，返回新的矩阵
	{
		Matrix<T> temp(m);
		return m.Cos();
	}
	template <class T> Matrix<T> Tan(Matrix<T>& m)//对所有元素取正切，返回新的矩阵
	{
		Matrix<T> temp(m);
		return temp.Tan();
	}
	template <class T> Matrix<T> Sinh(Matrix<T>& m)//对所有元素取双曲正弦，返回新的矩阵
	{
		Matrix<T> temp(m);
		return temp.Sinh();
	}
	template <class T> Matrix<T> Cosh(Matrix<T>& m)//对所有元素取双曲余弦，返回新的矩阵
	{
		Matrix<T> temp(m);
		return temp.Cosh();
	}
	template <class T> Matrix<T> Tanh(Matrix<T>& m)//对所有元素取双曲正切，返回新的矩阵
	{
		Matrix<T> temp(m);
		return temp.Tanh();
	}
	template <class T> Matrix<T> Log(Matrix<T>& m)//对所有元素取自然对数，返回新的矩阵
	{
		Matrix<T> temp(m);
		return temp.Log();
	}
	template <class T> Matrix<T> Log10(Matrix<T>& m)//对所有元素取10为底的对数，返回新的矩阵
	{
		Matrix<T> temp(m);
		return temp.Log10();
	}
	template <class T> Matrix<T> Exp(Matrix<T>& m)//对所有矩阵进行以e为底的指数运算，返回新的矩阵
	{
		Matrix<T> temp(m);
		return temp.Exp();
	}
	template <class T> Matrix<T> Sqrt(Matrix<T>& m)//对矩阵所有的元素取二次方根，返回新的矩阵
	{
		Matrix<T> temp(m);
		return temp.Sqrt();
	}
	template <class T> Matrix<T>& Matrix<T>::DotP(Matrix<T>& multiplicator)
	{
		if((this->m_col-multiplicator.ColNum())*(this->m_row-multiplicator.RowNum())!=0)
			throw Exception("Matrix","DotP","the dimensions of two matrices are not equal");
		unsigned long row, col;
		for(row=1;row<=m_row;row++)
			for(col=1;col<=m_col;col++)
			{
				this->SetValue(row,col,this->GetValue(row,col)*multiplicator.GetValue(row,col));
			}
		return *this;
	}
	template <class T> Matrix<T>& Matrix<T>::DotD(Matrix<T>& divisor)
	{
		if((this->m_col-divisor.ColNum())*(this->m_row-divisor.RowNum())!=0)
			throw Exception("Matrix","DotD","the dimensions of two matrices are not equal");
		unsigned long row, col;
		for(row=1;row<=m_row;row++)
			for(col=1;col<=m_col;col++)
			{
				if(this->ValueAbs(divisor.GetValue(row,col))<=this->m_Tolerance)
					throw Exception("Matrix","DotD","Divided by zero");
				else
					this->SetValue(row,col,this->GetValue(row,col)/divisor.GetValue(row,col));
			}
		return *this;
	}
	template <class T> Matrix<T> DotP(Matrix<T>& multiplicand, Matrix<T>& multiplicator)
	{
		Matrix<T> temp(multiplicand);
		return temp.DotP(multiplicator);
	}
	template <class T> Matrix<T> DotD(Matrix<T>& dividend,Matrix<T>& divisor)
	{
		Matrix<T> temp(dividend);
		return temp.DotD(divisor);
	}
	//以上是矩阵模板类的定义和实现
	///////////////////////////////////
	//以下是复数矩阵模板类的定义部分
	template <class T> 
	class CMatrix
	{
	private:
		unsigned long m_row;//存放矩阵的行数
		unsigned long m_col;//存放矩阵的列数
		T m_Tolerance;//存放矩阵的容许误差，小于等于这个数认为为零
		T* m_RealBuf;//矩阵实部的内存缓冲区
		T* m_ImagBuf;//矩阵虚部的内存缓冲区
		T ValueAbs(T value) const;//计算元素的绝对值
		T MaxAbs(unsigned long& row, unsigned long& rol, unsigned long k);//计算第k行和第k列之后的主元和位置
		T Adjust(T a);//
		void DoAdjust();//应用拉近技术，调正数值
	public:
		CMatrix();
		CMatrix(unsigned long row, unsigned long col);//定义一个row行，col列，实部和虚部都为零的复数矩阵
		CMatrix(CMatrix& cmatrix);//拷贝构造函数
		CMatrix(const char *FileName);
		CMatrix(T* RealBuf,T* ImagBuf, unsigned long row, unsigned long col);//通过两个一维数组和复数矩阵的行数和列数，构造一个复数矩阵
		CMatrix(unsigned long row, unsigned long col, Matrix<T> real,Matrix<T> imag);//通过实部矩阵和虚部矩阵构造复数矩阵
		~CMatrix();//析构函数，释放内存

		unsigned long RowNum();//返回矩阵的行数
		unsigned long ColNum();//返回矩阵的列数

		void SetValue(unsigned long row, unsigned long col, T real, T image);//根据下标给矩阵的元素赋值
		void SetValue(unsigned long row, unsigned long col, complex<T>& c);//根据下标，给矩阵的元素赋值
		complex<T> GetValue(unsigned long row, unsigned long col);//根据下标返回矩阵元素的值，返回类型：复数

		CMatrix<T>& Neg();//当前矩阵的实部和虚部都乘以－1，修改当前矩阵,返回当前矩阵
		CMatrix<T>& Tran();//对当前矩阵进行转置操作，修改当前矩阵,返回当前矩阵
		CMatrix<T> t();//对当前矩阵进行转置，产生一个新的矩阵，不改变原来矩阵，返回新的矩阵
		CMatrix<T>& Conj();//对当前矩阵取共轭操作，返回当前矩阵
		CMatrix<T> conj_tran();//对当前矩阵共轭，再取转置，返回新的矩阵		

		Matrix<T> Real();//返回当前矩阵的实部
		Matrix<T> Imag();//返回当前矩阵的虚部
		Matrix<T> Arg();//返回当前矩阵所有元素的幅角（argument）,For a complex number a + bi, the argument is equal to arctan(b/a). 
//复数矩阵的数学运算
		Matrix<T> Abs();//返回矩阵每个元素的模（modulus），The modulus of a complex number a + bi is sqrt(a^2 + b^2)
		Matrix<T> Norm();//返回矩阵每一个元素的范数（norm）,The norm of a complex number a + bi is a^2 + b^2
		CMatrix<T>& Sin();//对所有的元素取正弦，修改当前矩阵，返回当前矩阵
		CMatrix<T>& Cos();//对所有元素取余弦，修改当前矩阵，返回当前矩阵
		CMatrix<T>& Tan();//对所有元素取正切，修改当前矩阵，返回当前矩阵
		CMatrix<T>& Sinh();//对所有元素取双曲正弦，修改当前矩阵，返回当前矩阵
		CMatrix<T>& Cosh();//对所有元素取双曲余弦，修改当前矩阵，返回当前矩阵
		CMatrix<T>& Tanh();//对所有元素取双曲正切，修改当前矩阵，返回当前矩阵
		CMatrix<T>& Log();//对所有元素取自然对数，修改当前矩阵，返回当前矩阵
		CMatrix<T>& Log10();//对所有元素取10为底的对数，修改当前矩阵，返回当前矩阵
		CMatrix<T>& Exp();//对所有矩阵进行以e为底的指数运算，修改当前矩阵，返回当前矩阵
		CMatrix<T>& Sqrt();//对矩阵所有的元素取二次方根，修改当前矩阵，返回当前矩阵
//复数矩阵的数学运算结束
//复数矩阵的操作符重载
		complex<T> operator ()(unsigned long row, unsigned long col);//根据下标，返回一个复数
		CMatrix<T>& operator =(CMatrix<T>& rcmatrix);//当前矩阵等于right,返回当前矩阵
		CMatrix<T>& operator +=(CMatrix<T>& rcmatrix);//当前矩阵加上right,返回当前矩阵
		CMatrix<T> operator +(CMatrix<T>& rcmatrix);//当前矩阵加上right,返回一个新的矩阵
		CMatrix<T> operator -();//矩阵求负，产生新的矩阵，返回新的矩阵
		CMatrix<T>& operator -=(CMatrix<T>& rcmatrix);//当前矩阵减去rcmatrix,修改当前矩阵，返回当前矩阵
		CMatrix<T> operator -(CMatrix<T>& rcmatrix);//当前矩阵减去rcmatrix,产生新的矩阵，返回新的矩阵
		CMatrix<T>& operator *=(T k);//矩阵数乘k,修改当前矩阵，返回当前矩阵
		CMatrix<T> operator *(T k);//矩阵数乘k，返回新的矩阵
		CMatrix<T>& operator /=(T k);//矩阵的每一个元素都除以k,修改当前矩阵，返回当前矩阵
		CMatrix<T> operator /(T k);//矩阵的每一个元素都除以k,返回一个新的矩阵
		CMatrix<T>& operator *=(CMatrix<T>& rcmatrix);//当前矩阵乘以矩阵rcmatrix,修改当前矩阵，返回当前矩阵
		CMatrix<T> operator *(CMatrix<T>& rcmatrix);//当前矩阵乘以矩阵rcmatrix，产生新的矩阵，返回新的矩阵
		CMatrix<T>& operator ^=(T x);//矩阵的每一个元素都进行参数为y^x的幂运算，修改当前矩阵，返回当前矩阵
		CMatrix<T> operator ^(T x);//矩阵的每一个元素都进行按照y^x格式进行幂运算，返回一个新的矩阵
		friend ostream& operator <<(ostream& os,const CMatrix<T>& cmatrix);//屏幕输出当前矩阵
//以上匙复数模板类的操作符重载
		void SwapRows(unsigned long row1, unsigned long row2);//当前矩阵交换两行
		void SwapCols(unsigned long col1, unsigned long col2);//当前矩阵交换两列
		CMatrix<T>& Inv();//矩阵求逆，修改当前矩阵，返回当前矩阵
		CMatrix<T> operator ~();//矩阵求逆，返回新的矩阵
		void GenI();//将当前矩阵设置成为但为矩阵，主元实部为1，虚部为零，其它实部虚部都为零
		bool IsSquare();//当前矩阵是否为方阵
		bool IsSym();//当前矩阵是否为对称矩阵
		void SaveData(const char *FileName);
		void Save(const char *FileName);
		void SetDim(unsigned long rowNum, unsigned long colNum);
////the dot product and dot divide operation of complex matrix
		CMatrix<T>& DotP(CMatrix<T>& multiplicator);
		CMatrix<T>& DotD(CMatrix<T>& divisor);
	};
	//以下是复数矩阵模板类的实现部分
	template <class T> CMatrix<T>::CMatrix(T* RealBuf, T* ImagBuf,unsigned long row, unsigned long col)
	{
		m_row=row;
		m_col=col;
		m_Tolerance=(T)(1.0e-20);
		m_RealBuf= new T[m_row*m_col];
		m_ImagBuf= new T[m_row*m_col];
		unsigned long i,j;

		for(i=1;i<=m_row;i++)
			for(j=1;j<=m_col;j++)
			{
				SetValue(i,j,RealBuf[(i-1)*m_col+j-1],ImagBuf[(i-1)*m_col+j-1]);
			}
	}
	template <class T> CMatrix<T>::CMatrix()
	{
		m_row=0;
		m_col=0;
		m_Tolerance=(T)(0.000001);
		m_RealBuf= NULL;
		m_ImagBuf= NULL;
	}
	template <class T> CMatrix<T>::CMatrix(unsigned long row, unsigned long col)
	{
		m_row=row;
		m_col=col;
		m_Tolerance=(T)(0.000001);
		m_RealBuf= new T[m_row*m_col];
		m_ImagBuf= new T[m_row*m_col];
		unsigned long k;
		for(k=0;k<m_row*m_col;k++)
		{
			m_RealBuf[k]=(T)0.0;
			m_ImagBuf[k]=(T)0.0;
		}
	}
	template <class T> CMatrix<T>::CMatrix(CMatrix<T>& cmatrix)
	{
		m_row=cmatrix.RowNum();
		m_col=cmatrix.ColNum();
		m_Tolerance=(T)(0.000001);
		m_RealBuf= new T[m_row*m_col];
		m_ImagBuf= new T[m_row*m_col];
		unsigned long row, col;
		for(row=1;row<=m_row; row++)