matrix.hpp

随机数类,产生随机数

// matrix.hpp
// declaration and implementation of template class Matrix 
// <Shugong Wang @ CAREERI/CAS> email: [e-mail]
// August 14, 2003
// Ver 1.0
// Modification history
// <Xin Li @ CAREERI/CAS> email: [e-mail]
// August, 2005


#ifndef __MATRIX_HPP__
#define __MATRIX_HPP__

// To program is just to understand!
#include <iostream>
#include <complex>
#include <fstream>
#include <cmath>
#include "conio.h"

using namespace std;

namespace ldas
{
	//以下这一部分定义了异常类，用来实现在数据同化计算中的异常处理
	template <class T> class Matrix;
	template <class T> class CMatrix;
	//Print Matrix to screen
	template <class T> ostream& operator <<(ostream& os,const Matrix<T>& matrix);
	template <class T> ostream& operator <<(ostream& os,const CMatrix<T>& cmatrix);
	//the add,minus,time,divid of a number and a matrix
	template <class T> Matrix<T> operator + (T a, Matrix<T>& m);
	template <class T> Matrix<T> operator - (T a, Matrix<T>& m);
	template <class T> Matrix<T> operator * (T a, Matrix<T>& m);
	template <class T> Matrix<T> operator / (T a, Matrix<T>& m);
	//mathematic operations of real matrix
	template <class T> Matrix<T> Sin(Matrix<T>& m);//对所有的元素取正弦，返回新的矩阵
	template <class T> Matrix<T> Cos(Matrix<T>& m);//对所有元素取余弦，返回新的矩阵
	template <class T> Matrix<T> Tan(Matrix<T>& m);//对所有元素取正切，返回新的矩阵
	template <class T> Matrix<T> Sinh(Matrix<T>& m);//对所有元素取双曲正弦，返回新的矩阵
	template <class T> Matrix<T> Cosh(Matrix<T>& m);//对所有元素取双曲余弦，返回新的矩阵
	template <class T> Matrix<T> Tanh(Matrix<T>& m);//对所有元素取双曲正切，返回新的矩阵
	template <class T> Matrix<T> Log(Matrix<T>& m);//对所有元素取自然对数，返回新的矩阵
	template <class T> Matrix<T> Log10(Matrix<T>& m);//对所有元素取10为底的对数，返回新的矩阵
	template <class T> Matrix<T> Exp(Matrix<T>& m);//对所有矩阵进行以e为底的指数运算，返回新的矩阵
	template <class T> Matrix<T> Sqrt(Matrix<T>& m);//对矩阵所有的元素取二次方根，返回新的矩阵
	template <class T> Matrix<T> Asin(Matrix<T>& m);//计算矩阵的反正弦函数，返回一个新的矩阵
	template <class T> Matrix<T> Acos(Matrix<T>& m);//计算矩阵的反余弦函数，返回一个新的矩阵
	template <class T> Matrix<T> Atan(Matrix<T>& m);//计算矩阵的反正切函数，返回一个新的矩阵
	//methematic operations of complex matrix
	template <class T> Matrix<T> Abs(CMatrix<T>& cm);//返回矩阵每个元素的模（modulus），The modulus of a complex number a + bi is sqrt(a^2 + b^2)
	template <class T> Matrix<T> Norm(CMatrix<T>& cm);//返回矩阵每一个元素的范数（norm）,The norm of a complex number a + bi is a^2 + b^2
	template <class T> CMatrix<T> Sin(CMatrix<T>& cm);//对所有的元素取正弦，返回新的矩阵
	template <class T> CMatrix<T> Cos(CMatrix<T>& cm);//对所有元素取余弦，返回新的矩阵
	template <class T> CMatrix<T> Tan(CMatrix<T>& cm);//对所有元素取正切，返回新的矩阵
	template <class T> CMatrix<T> Sinh(CMatrix<T>& cm);//对所有元素取双曲正弦，返回新的矩阵
	template <class T> CMatrix<T> Cosh(CMatrix<T>& cm);//对所有元素取双曲余弦，返回新的矩阵
	template <class T> CMatrix<T> Tanh(CMatrix<T>& cm);//对所有元素取双曲正切，返回新的矩阵
	template <class T> CMatrix<T> Log(CMatrix<T>& cm);//对所有元素取自然对数，返回新的矩阵
	template <class T> CMatrix<T> Log10(CMatrix<T>& cm);//对所有元素取10为底的对数，返回新的矩阵
	template <class T> CMatrix<T> Exp(CMatrix<T>& cm);//对所有矩阵进行以e为底的指数运算，返回新的矩阵
	template <class T> CMatrix<T> Sqrt(CMatrix<T>& cm);//对矩阵所有的元素取二次方根，返回新的矩阵	
	//template <class T> CMatrix<T> Asin(CMatrix<T>& cm);//计算矩阵的反正弦函数，返回一个新的矩阵
	//template <class T> CMatrix<T> Acos(CMatrix<T>& cm);//计算矩阵的反余弦函数，返回一个新的矩阵
	//template <class T> CMatrix<T> Atan(CMatrix<T>& cm);//计算矩阵的反正切函数，返回一个新的矩阵
	/////dot product and dot divide
	template <class T> Matrix<T> DotP(Matrix<T>& multiplicand, Matrix<T>& multiplicator);
	template <class T> Matrix<T> DotD(Matrix<T>& dividend, Matrix<T>& divisor);
	template <class T> CMatrix<T> DotP(CMatrix<T>& multiplicand, CMatrix<T>& multiplicator);
	template <class T> CMatrix<T> DotD(CMatrix<T>& dividend, CMatrix<T>& divisor);
	
	////////////////////////////////////////////////////
	class Exception
	{
	private:
		const char* m_Info;//存放异常信息
		const char* m_ClassName;//发生异常的类的名字
		const char* m_MethodName;//发生异常的方法名字
	public:
		void What(void){
			cout<<"Exception Occurred!"<<endl;
			cout<<"Class Name: "<<m_ClassName<<endl;
			cout<<"Method Name: "<<m_MethodName<<endl;
			cout<<"Exception Info: "<<m_Info<<endl; 
			cout<<"Press any key to terminate...";
			getch();
		}
		Exception(const char * ClassName, const char* MethodName, const char* Info){
			m_Info=Info;
			m_ClassName=ClassName;
			m_MethodName=MethodName;
		}
		~Exception(){
		}
	};
	//以上部分是异常类的定义以及实现的部分

	//以下部分是矩阵模板类的定义和实现
	template<class T> 
	class Matrix
	{
	private:
		T* m_buf;//存放数据的缓冲区
		T m_Tolerance;//数值比较时的容忍误差
		unsigned long m_row;//矩阵的行数
		unsigned long m_col;//矩阵的列数
	 	T RowMaxAbs(unsigned long& max_row, unsigned long cur_col,unsigned long start_row);//计算第k列的主元的位置和值
		T ValueAbs(T value);//计算绝对值
		T MaxAbs(unsigned long& row, unsigned long& col, unsigned long k);//计算第k行和第k列之后的主元的位置和值
		T DetAndRank(unsigned long &rank);//计算行列式的值和秩
		T Adjust(T a);//对数字小数点后五位进行四舍五入,取近似
		void DoAdjust();//对矩阵元素进行拉近
	public:
		Matrix();//构造一个空的矩阵，注意在使用之前一定要通过SetDim函数或者赋值操作，设定矩阵的维数
		Matrix(unsigned long row, unsigned long col);//使用行数和列数初始化矩阵
		Matrix(const Matrix<T>& matrix);//拷贝构造函数，生成的矩阵和matrix完全一致
		Matrix(T* const buf,unsigned long row, unsigned long col);//用一个数组初始化
		Matrix(const char* FileName);//通过存放在文件中的数据构造矩阵，数据文件有特定的格式
		~Matrix();
		unsigned long RowNum(void) const;//获得矩阵的行数
		unsigned long ColNum(void) const;//获得矩阵的列数
		bool IsSym(void);//判断矩阵是否为对称矩阵 is symmetry
		bool IsSquare(void);//判断矩阵是否为方阵 is square
		bool IsEqual(Matrix<T>& matrix);//判断两个矩阵是否相等
		bool IsI(void);//判断矩阵是否是单位矩阵
		bool IsSingular();//反断矩阵是否为奇异矩阵（退化矩阵）
		
		Matrix<T> ExtractColV(unsigned long col) const; //提取某一列向量
		Matrix<T> ExtractRowV(unsigned long row) const;//提取某一行向量
		void SetColV(const unsigned long col, const Matrix<T>& colElement);
		void SetRowV(const unsigned long row, const Matrix<T>& rowElement);
		
		T GetValue(unsigned long row, unsigned long col) const;//提取某一元素的值
		void SetValue(unsigned long row, unsigned long col,T value);//设置某一元素的值
		void SwapRows(unsigned long row1, unsigned long row2);//交换某两行
		void SwapCols(unsigned long col1, unsigned long col2);//交换某两列
		void GenI();//将当前矩阵设置成为单位矩阵
		void GenOnes();//将当前矩阵的所有元素都设置成为1
		void GenZeros();//将当前矩阵的所有元素的值都是之成为0
		Matrix<T>& Inv(void);//计算矩阵的逆矩阵,修改当前矩阵，返回当前矩阵
		Matrix<T>& Tran(void);//计算矩阵的转置矩阵，不产生新的矩阵
		Matrix<T> t(void);//计算矩阵的转置矩阵，产生新的矩阵
		Matrix<T>& Neg(void);//计算当前矩阵的负矩阵， 不产生新的矩阵，返回当前矩阵
		T Det();//计算矩阵的行列式的值。
		unsigned long Rank();//计算矩阵的秩
		T Trace();//如果矩阵是方阵，计算矩阵的迹
		Matrix<T> Cov();//计算当前矩阵的协方差矩阵，认为当前矩阵的每一个列向量为一个随机变量集，返回新的矩阵。
		T operator ()(unsigned long row, unsigned long col);//返回矩阵元素的值
		Matrix<T> operator ~();//计算逆矩阵，不产生新的矩阵，修改当前矩阵
		Matrix<T>& operator =(const Matrix<T>& rmatrix);
		// Matrix<T>& operator =(Matrix<T>& rmatrix);
		Matrix<T>& operator= (const T& scalar);

		
		Matrix<T> operator -(T a);//矩阵的每个元素都减去一个常数a
		Matrix<T> operator -(Matrix<T>& matrix);//两个矩阵相减，返回一个新的矩阵
		Matrix<T>& operator -=(T a);//当前矩阵的每一个元素都减去常数a，返回当前矩阵
		Matrix<T>& operator -=(Matrix<T>& matrix);//当前矩阵减去矩阵matrix,返回当前矩阵
		Matrix<T> operator -(void);//计算当前矩阵的负矩阵，产生新的矩阵
		
		Matrix<T> operator +(T a);//矩阵的每一个元素都加上常数a，返回一个新的矩阵
		Matrix<T> operator +(Matrix<T>& matrix);//两个矩阵相加，返回一个新的矩阵
		Matrix<T>& operator +=(T a);//矩阵的每一个元素都加上常数a，返回当前矩阵
		Matrix<T>& operator +=(Matrix<T>& matrix);//当前矩阵加上矩阵matrix，返回当前矩阵

		Matrix<T> operator *(T a);//当前矩阵乘以常数a，返回一个新的矩阵
		Matrix<T> operator *(Matrix<T>& matrix);//当前矩阵乘以矩阵matrix，返回一个新的矩阵
		Matrix<T>& operator *=(T a);//当前矩阵乘以常数a,返回当前矩阵
		Matrix<T>& operator *=(Matrix<T>& matrix);//当前矩阵乘以矩阵matrix,返回当前矩阵

		Matrix<T> operator /(T a);//当前矩阵的所有元素除以常数a,返回一个新的矩阵
		Matrix<T>& operator /=(T a);//当前矩阵的所有元素都除以常数a,返回当前矩阵
		
		Matrix<T>& operator ^=(double x);//当前矩阵的每一个元素都取y^x,返回当前矩阵
		Matrix<T> operator ^(double x);//意义同上，返回新的矩阵

		friend ostream& operator <<(ostream& os,const Matrix<T>& matrix);

		friend Matrix<T> operator + (T a, Matrix<T>& m);
		friend Matrix<T> operator - (T a, Matrix<T>& m);
		friend Matrix<T> operator * (T a, Matrix<T>& m);
		friend Matrix<T> operator / (T a, Matrix<T>& m);

		void Save(const char* FileName);
		void SetDim(unsigned long rowNum, unsigned long colNum);
	
	// Matrix decomposition
		void LU(Matrix<T>& L, Matrix<T>& U);//对矩阵进行三角分解，使用Doolittle算法
		Matrix<T> Chold();		// Cholesky decomposition

		
//mathematic operation
		Matrix<T>& Sin();//对所有的元素取正弦，修改当前矩阵，返回当前矩阵
		Matrix<T>& Cos();//对所有元素取余弦，修改当前矩阵，返回当前矩阵
		Matrix<T>& Tan();//对所有元素取正切，修改当前矩阵，返回当前矩阵
		Matrix<T>& Sinh();//对所有元素取双曲正弦，修改当前矩阵，返回当前矩阵
		Matrix<T>& Cosh();//对所有元素取双曲余弦，修改当前矩阵，返回当前矩阵
		Matrix<T>& Tanh();//对所有元素取双曲正切，修改当前矩阵，返回当前矩阵
		Matrix<T>& Log();//对所有元素取自然对数，修改当前矩阵，返回当前矩阵，如果某个矩阵元素小于等于零，抛出异常
		Matrix<T>& Log10();//对所有元素取10为底的对数，修改当前矩阵，返回当前矩阵，如果某个矩阵元素小于等于零，抛出异常
		Matrix<T>& Exp();//对所有矩阵进行以e为底的指数运算，修改当前矩阵，返回当前矩阵	
		Matrix<T>& Sqrt();//对所有元素取平方根，修改当前矩阵，返回当前矩阵，如果某个矩阵元素小于零，抛出异常
		Matrix<T>& Asin();//对当前矩阵所有元素取反正弦，修改当前矩阵，返回当前矩阵
		Matrix<T>& Acos();//对当前矩阵所有元素取反余弦，修改当前矩阵，返回当前矩阵
		Matrix<T>& Atan();//对当前矩阵所有元素取反正切，修改当前矩阵，返回当前矩阵
///// matrix dot procuct, dot divide operation
		Matrix<T>& DotP(Matrix<T>& multiplicator);//两个矩阵点乘，要求两个矩阵行数和列数相同，如果维数不同，抛出异常
		Matrix<T>& DotD(Matrix<T>& divisor);//两个矩阵点除，要求两个矩阵行数和列数相同，如果维数不同，抛出异常，如果被除数存在零元素，抛出异常
	};
	//以上式矩阵模板类的定义部分，一下是实现部分
	template <class T> Matrix<T>::Matrix()
	{
		m_col =0;
		m_row =0;
		m_Tolerance =(T)0.000001;
		m_buf=NULL;		
	}
	template <class T> Matrix<T>::Matrix(unsigned long row, unsigned long col)
	{
		m_row=row;
		m_col=col;
		m_Tolerance=(T)0.000001;
		unsigned long k;
		m_buf=new T[m_row*m_col];
		for(k=0;k<m_row*m_col;k++)
		{
			m_buf[k]=(T)0.0;
		}
	}

	template <class T> Matrix<T>::Matrix(const Matrix<T>& matrix)//拷贝构造函数
	{	
		unsigned long i,j;
		m_row=matrix.RowNum();
		m_col=matrix.ColNum();
		m_Tolerance=(T)0.000001;
		m_buf=new T[m_row*m_col];
		for(i=1;i<=m_row;i++)//将矩阵matrix的元素值复制到当前矩阵内
			for(j=1;j<=m_col;j++)
				this->SetValue(i,j,matrix.GetValue(i,j));  
	}

	template <class T> Matrix<T>::Matrix(T* const buf, unsigned long row, unsigned long col)
	{
		unsigned long k;
		m_row=row;
		m_col=col;
		m_Tolerance=(T)0.000001;
		m_buf=new T[row*col];
		for(k=0;k<m_row*m_col;k++)
			m_buf[k]=buf[k];
	}
	template <class T> Matrix<T>::Matrix(const char* FileName)
	{
		ifstream in;
		in.open(FileName) ; 
		if(!in)//如果文件不存在或者不能打开，抛出异常
			throw Exception("Matrix","Matrix(const char* FileName)","Can not open data file");
		in>>this->m_row;//读入行数
		in>>this->m_col;//读入列数
		m_Tolerance=(T)0.000001;
		m_buf=new T[m_row*m_col];
		T temp;
		unsigned long row,col;
		for(row=1;row<=m_row;row++)
			for(col=1;col<=m_col;col++)
			{
				in>>temp;
				this->SetValue(row,col,temp); 
			}	
	}

	template <class T> void Matrix<T>::SetDim(unsigned long rowNum,unsigned long colNum)
	{
		this->m_col=colNum;
		this->m_row=rowNum;
		this->m_buf =new T[rowNum*colNum];
		unsigned long k;
		for(k=0;k<m_row*m_col;k++)
			m_buf[k]=(T)0.0;
	}
	template <class T> Matrix<T>& Matrix<T>::Exp()
	{
		unsigned long row,col;
		for(row=1;row<=m_row;row++)
			for(col=1;col<=m_col;col++)
				this->SetValue(row,col,exp(this->GetValue(row,col)));
		return (*this);
	}
	template <class T> Matrix<T>& Matrix<T>::Sin()//对所有的元素取正弦，修改当前矩阵，返回当前矩阵
	{
		unsigned long row,col;
		for(row=1;row<=m_row;row++)
			for(col=1;col<=m_col;col++)
				this->SetValue(row,col,sin(this->GetValue(row,col)));
		return (*this);
	}
	template <class T> Matrix<T>& Matrix<T>::Cos()//对所有元素取余弦，修改当前矩阵，返回当前矩阵
	{
		unsigned long row,col;
		for(row=1;row<=m_row;row++)
			for(col=1;col<=m_col;col++)
				this->SetValue(row,col,cos(this->GetValue(row,col)));
		return (*this);
	}
	template <class T> Matrix<T>& Matrix<T>::Tan()//对所有元素取正切，修改当前矩阵，返回当前矩阵
	{
		unsigned long row,col;
		for(row=1;row<=m_row;row++)
			for(col=1;col<=m_col;col++)
				this->SetValue(row,col,tan(this->GetValue(row,col)));
		return (*this);
	}
	template <class T> Matrix<T>& Matrix<T>::Sinh()//对所有元素取双曲正弦，修改当前矩阵，返回当前矩阵
	{
		unsigned long row,col;
		for(row=1;row<=m_row;row++)
			for(col=1;col<=m_col;col++)
				this->SetValue(row,col,sinh(this->GetValue(row,col)));
		return (*this);
	}
	template <class T> Matrix<T>& Matrix<T>::Cosh()//对所有元素取双曲余弦，修改当前矩阵，返回当前矩阵
	{
		unsigned long row,col;
		for(row=1;row<=m_row;row++)
			for(col=1;col<=m_col;col++)
				this->SetValue(row,col,cosh(this->GetValue(row,col)));
		return (*this);
	}
	template <class T> Matrix<T>& Matrix<T>::Tanh()//对所有元素取双曲正切，修改当前矩阵，返回当前矩阵
	{
		unsigned long row,col;
		for(row=1;row<=m_row;row++)
			for(col=1;col<=m_col;col++)
				this->SetValue(row,col,sinh(this->GetValue(row,col)));
		return (*this);
	}
	template <class T> Matrix<T>& Matrix<T>::Log()//对所有元素取自然对数，修改当前矩阵，返回当前矩阵
	{
		unsigned long row,col;
		for(row=1;row<=m_row;row++)
			for(col=1;col<=m_col;col++)
			{
				if(this->GetValue(row,col)<=0)
					throw Exception("Matrix","Log","one or more element of matrix is less than or equal to zero");
			}
		for(row=1;row<=m_row;row++)
			for(col=1;col<=m_col;col++)
				this->SetValue(row,col,log(this->GetValue(row,col)));
		return (*this);
	}
	template <class T> Matrix<T>& Matrix<T>::Log10()//对所有元素取10为底的对数，修改当前矩阵，返回当前矩阵
	{
		unsigned long row,col;
		for(row=1;row<=m_row;row++)
			for(col=1;col<=m_col;col++)
			{
				if(this->GetValue(row,col)<=0)
					throw Exception("Matrix","Log10","one or more element of matrix is less than or equal to zero");
			}
		for(row=1;row<=m_row;row++)
			for(col=1;col<=m_col;col++)
				this->SetValue(row,col,log10(this->GetValue(row,col)));
		return (*this);
	}
	template <class T> Matrix<T>& Matrix<T>::Sqrt()
	{
		unsigned long row,col;
		for(row=1;row<=m_row;row++)
			for(col=1;col<=m_col;col++)
			{
				if(this->GetValue(row,col)<0)
					throw Exception("Matrix","Sqrt","one or more element of matrix is less than zero");
			}
		for(row=1;row<=m_row;row++)
			for(col=1;col<=m_col;col++)
				this->SetValue(row,col,sqrt(this->GetValue(row,col)));
		return (*this);
	}
	template <class T> Matrix<T>& Matrix<T>::Asin()
	{