matrix.hpp

随机数类,产生随机数

	}
	template <class T> void Matrix<T>::GenZeros()//产生所有元素都为零的矩阵
	{
		unsigned long row, col;
		for(row=1; row<= this->m_row; row++ )
			for(col=1; col<= this->m_col; col++)
			{
				SetValue(row,col,0.0);
			}
	}
	template <class T> T Matrix<T>::RowMaxAbs(unsigned long& max_row,unsigned long cur_col, unsigned long start_row)//用在计算逆矩阵的函数中
	{//私有函数
		unsigned long i;
		T temp;
		temp=GetValue(start_row,cur_col);//  *this(start_row,cur_col);//初始化变量
		max_row=start_row;
		for(i=start_row;i<=m_row;i++)
		{
			if(ValueAbs(temp)<=ValueAbs(GetValue(i,cur_col)))
			{
				temp=GetValue(i,cur_col);
				max_row=i;
			}
		}
		return temp;
	}
	template <class T> T Matrix<T>::MaxAbs(unsigned long &row,unsigned long &col,unsigned long k)//计算的第k行和第k列之后的主元的位置和值，私有函数
	{//私有函数
		T temp=0.0;
		unsigned long i,j;
		for(i=k;i<=m_row;i++)
			for(j=k;j<=m_col;j++)
			{
				if(temp<ValueAbs(GetValue(i,j)))//如果有比零大的元素，临时变量就等于这个元素
				{
					row=i;
					col=j;
					temp=ValueAbs(GetValue(i,j));
				}
			}
		return temp;
	}
	template <class T> T Matrix<T>::DetAndRank(unsigned long &rank)
	{
		if(this->m_col!=this->m_row||this->m_col==0)//如果这个矩阵不是方阵或者矩阵为空，抛出异常
			 throw Exception("Matrix","DetAndRank","The matrix is not square or its dimension is zero.");
		else{
			unsigned long i,j,k,is,js;//其中is和js用来存放选中主元的位置
			T f,detv,q,d;
			Matrix<T> temp(*this);
			rank=0;//矩阵的秩初值设置为零
			f=1.0;
			detv=1.0;
			for(k=1;k<m_row;k++)
			{
				q=MaxAbs(is,js,k);//q是k行k列之后的主元的数值，is和js是选中主元的位置。
                if(q<=m_Tolerance)//如果整个矩阵最大主元是零，行列式的值
					return 0.0;
				rank++;//
				if(is!=k)//如果主元不在当前行，交换行
				{
					f=-f;
					temp.SwapRows(k,is);
				}
				if(js!=k)//如果主元不在当前列，交换列
				{
					f=-f;
					temp.SwapCols(k,js);
				}
				q=temp.GetValue(k,k);
				detv*=q;
				for(i=k+1;i<=m_row;i++)
				{
					d=temp.GetValue(i,k)/q;
					for(j=k+1;j<=m_row;j++)
					{
						temp.SetValue(i,j,(temp.GetValue(i,j)-d*temp.GetValue(k,j)));
					}
				}
			}//k循环结束了，呵呵
			q=temp.GetValue(m_row,m_col);
			if(ValueAbs(q)>m_Tolerance)
			{
				rank++;
				return Adjust(f*detv*q);
			}
			else{
				return 0.0;
			}
		}
	}
	template <class T> T Matrix<T>::Det ()//计算矩阵的行列式的值
	{
		T temp;
		unsigned long rank;
		temp=DetAndRank(rank);
		return Adjust(temp);
	}
	template <class T> unsigned long Matrix<T>::Rank()//计算矩阵的秩
	{
		unsigned long rank;
		this->DetAndRank(rank);
		return rank;
	}

	template <class T>
		void Matrix<T>::LU(Matrix<T>& L,Matrix<T>& U)
	{//注意，输入参数L,U的大小应该和当前矩阵相同
		//unsigned long k,i,j,jj;
		T det;
		/*T temp;*/
		Matrix<T> matrix(*this);
		if(!this->IsSquare())//如果当前矩阵不是方阵，抛出异常，程序结束
			throw Exception("Matrix","LU","the matrix is not square.");
		det=matrix.Det();//计算当前矩阵的行列式的数值，如果行列式的数值为零，那么退出程序
		if(ValueAbs(det)<this->m_Tolerance)//如果行列式的值等于零，抛出异常
			throw Exception("Matrix","LU","the matrix is singular.");
		L=*this;
		U=*this;
		L.GenZeros();
		U.GenZeros(); 
		unsigned long i,r,k;
		T temp;
		for(i=1;i<=this->m_row;i++)
		{
			U.SetValue(1,i,this->GetValue(1,i));
			//L.SetValue(i,1,this->GetValue(i,1)/U.GetValue(1,1));
		}
		L.SetValue(1,1,1); 
		for(i=2;i<=this->m_row;i++)
		{
			//U.SetValue(1,i,this->GetValue(1,i));
			L.SetValue(i,1,this->GetValue(i,1)/U.GetValue(1,1));
			L.SetValue(i,i,1); 
		}
		for(r=2;r<=this->m_col;r++)
		{
			for(i=r;i<=this->m_col;i++)
			{
				temp=0;
				for(k=1;k<=r-1;k++)
					temp+=L(r,k)*U(k,i);
				U.SetValue(r,i,(this->GetValue(r,i)-temp));
			}
			for(i=r+1;i<=m_col;i++)
			{
				if(r!=m_row)
				{
					temp=0;
					for(k=1;k<=(r-1);k++)
						temp+=L(i,k)*U(k,r);
					//cout<<U(r,r)<<endl;
					L.SetValue(i,r,(this->GetValue(i,r)-temp)/U(r,r)); 
				}
			}
		}
	}

	template <class T>
	Matrix<T> Matrix<T>::Chold()
	{
		if ( !this->IsSym() )
			throw Exception("Matrix","Chold","the matrix is not symmetric matrix");

		Matrix<T> matrix(*this);
		int n = this->RowNum();
		int i, j;
		T s;
	
		for (i=1; i<=n; i++)
		{
			for (j=i; j<=n; j++)
			{	
				s = this->GetValue(i,j);
				for (int m=i-1; m>=1; m--)
					s -= matrix.GetValue(i,m)*matrix.GetValue(j,m);
				if (i==j)
				{
					if (s<m_Tolerance)
						throw Exception("Matrix","Chold","not positive definite");
					else
						matrix.SetValue(i, i, sqrt(s));
				}
				else
					matrix.SetValue(j, i, s/matrix.GetValue(i,i));
			}
		}
		
		for (j=2; j<=n; j++)
			for (i=1; i<j; i++)
				matrix.SetValue(i, j, T(0)); 
						
		return matrix;
	}

	template <class T> Matrix<T>& Matrix<T>::Inv()//应用高斯-约当列主元消去法，计算当前矩阵的逆矩阵，返回一个新的矩阵
	{
		
		if(!this->IsSquare())//如果当前矩阵不是方阵，则抛出异常，不能进行求逆运算
			throw Exception("Matrix","Inv","The current matrix is not square, can't calculate inverse.");
		//cout<<"妈的"<<endl<<*this<<endl;
		if(ValueAbs(this->Det())<=this->m_Tolerance)//如果矩阵是奇异矩阵，抛出异常，不能进行求逆运算
			throw Exception("Matrix","Inv","the current matrix is singular, can't calculate inverse");
		else{
			unsigned long i,j,k,dim;
			T det;//记录矩阵的行列式的值
			dim=m_row+1;
			unsigned long* const Ip=new unsigned long[dim];//用来存放主元位置
			//Ip=new[dim];//初始化，可以作为数组使用
			
			Matrix<T> matrix(*this);//临时矩阵
			Matrix<T> m(*this);
			//matrix=*this;//将当前矩阵设置到临时矩阵中
			//m=*this;
			m.GenZeros(); 
			det=1.0;//(1)
			for(k=1;k<=m_col;k++)
			{
				//(2)按列选主元
				T c0,h;
				unsigned long max_row;
				c0=matrix.RowMaxAbs(max_row,k,k);
			//	cout<<c0<<" "<<max_row<<endl;
				Ip[k]=max_row;//记录
				if(c0<m_Tolerance)//(3)奇异矩阵，抛出异常
					Exception("Matrix","Inv","the current matrix is singular, can't calculate inverse");
				if(max_row!=k)//(4)如果最大行不等于K,交换行，否则，下一步
				{
					matrix.SwapRows(max_row,k);
					det=-det;
				}
				det=det*c0;//(5)
				//(6)
				//matrix. (k,k)=1/c0;
				matrix.SetValue(k,k,1/c0);
				h=matrix(k,k);
				for(i=1;i<=m_row;i++)
				{
					if(i!=k)
					{
						m.SetValue(i,k,(-matrix(i,k)*h));
						matrix.SetValue(i,k,m(i,k));
					}
				}
				//end(6)
				//(7)消元计算
				for(i=1;i<=m_row;i++)
					for(j=1;j<=m_col;j++)
					{
						if((i!=k)&&(j!=k))//不是主元
						{
							matrix.SetValue(i,j,(matrix(i,j)+m(i,k)*matrix(k,j)));
						}
					}
				//end(7)
				//(8)计算主行
				for(j=1;j<=m_col;j++)
				{
					if(j!=k)
					{
						matrix.SetValue(k,j,(matrix(k,j)*h));
					}
				}
				//(8)end
			}//k循环结束
			unsigned long t;
			for(k=(m_col-1);k>=1;k--)
			{
				t=Ip[k];
				if(k!=t)
				{
					matrix.SwapCols(t,k);
				}
			}
			matrix.DoAdjust(); 
			(*this)=matrix;
			return *this;
		}
	}//Inv()结束
	template <class T> bool Matrix<T>::IsSingular()
	{
		T det=this->Det() ;
		if(this->ValueAbs(det)<=this->m_Tolerance)
			return true;
		else 
			return false;
	}
	template <class T> Matrix<T> Matrix<T>::Cov()
	{
		Matrix<T> cov(m_col,m_col);//矩阵自动初始化为零矩阵
		T* mean= new T[m_col];
		unsigned long row, col,k;
		for(col=1;col<=m_col;col++)//计算每一个随即变量的平均值（期望）
		{
			mean[col-1]=(T)0.0;
			for(row=1;row<=m_row;row++)
				mean[col-1]+=this->GetValue(row,col); 
			mean[col-1]/=m_row;
		}
		T temp;
		for(row=1;row<=m_col;row++)
			for(col=row;col<=m_col;col++)
			{
				temp=0;
				for(k=1;k<=m_row;k++)
				{
					temp+=(this->GetValue(k,row)-mean[row-1])*(this->GetValue(k,col)-mean[col-1]);
				}
				cout<<endl;
				temp/=(m_row-1);
				cov.SetValue(row,col,temp);
			}
		for(row=1;row<=cov.RowNum();row++)
			for(col=1;col<row;col++)
			{
				cov.SetValue(row,col,cov.GetValue(col,row));
			}
		delete[] mean;
		return cov;
	}
	template <class T> Matrix<T> Matrix<T>::operator ~()//对当前矩阵进行逆运算，运算结果保存在当前矩阵内。
	{
		Matrix<T> temp(*this);
		/**this=temp.Inv();*/
		return (temp.Inv());
	}
	template <class T> Matrix<T>& Matrix<T>::operator -=(T a)
	{
		return (*this+=(-a));
	}
	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)); 
			}
		return *this;
	}
	template <class T> Matrix<T> Matrix<T>::operator +(T a)//当前矩阵的所有元素都加上常数a，返回一个新的矩阵
	{
		Matrix<T> matrix(this->m_row,this->m_col);
		matrix=*this;
		matrix+=a;
		return matrix;
	}
	template <class T> Matrix<T>& Matrix<T>::operator +=(Matrix<T>& matrix)//当前矩阵加上矩阵matrix,返回当前矩阵
	{
		unsigned long row,col;
		if((m_row!=matrix.RowNum())||(m_col!=matrix.ColNum()))//如果两个矩阵的规模不同，不能进行相加运算，抛出异常
			throw Exception("Matrix","Operator +=","the two matrice are not in same size, can not perform addtion operation.");
		for(row=1;row<=m_row;row++)
			for(col=1;col<=m_col;col++)
			{
				this->SetValue(row,col,(matrix(row,col)+(*this)(row,col))); 
			}
		return *this;
	}
	template <class T> Matrix<T> Matrix<T>::operator +(Matrix<T>& matrix)//当前矩阵加上matrix,返回一个新的矩阵
	{
		if(m_row!=matrix.RowNum()||m_col!=matrix.ColNum())
			throw Exception("Matrix","Operator +","the two matrice are not in same size, can not perform addtion operation.");
		Matrix<T> temp(this->m_row,this->m_col);
		temp=*this;
		temp+=matrix;
		return temp;
	}
	template <class T> Matrix<T> Matrix<T>::operator -(void)//当前矩阵的所有的元素都取复制负值，返回当前矩阵
	{
		Matrix<T> matrix(*this);
		matrix.Neg();
		return matrix;
	}
	template <class T> Matrix<T>& Matrix<T>::operator -=(Matrix<T>& matrix)//当前矩阵减去矩阵matrix，返回当前矩阵
	{
		(*this)+=(-matrix);
		return *this;
	}
	template <class T> Matrix<T> Matrix<T>::operator -(T a)//当前矩阵的所有元素都减去常数a，返回一个新的矩阵
	{