matrix.inl

里面有非常丰富的数值分析函数。编写很规范

// Matrix.inl		矩阵模板类函数(方法)定义
// Ver 1.0.0.0
// 版权所有(C) 何渝, 2002
// 最后修改: 2002.5.31

#ifndef _MATRIX_INL
#define _MATRIX_INL

//矩阵乘法函数
template <class _Tyout, class _Tylhs, class _Tyrhs>	//最后结果在mOut中
matrix<_Tyout>& MatrixMultiply(matrix<_Tyout>& mOut, const matrix<_Tylhs>& lhs, const matrix<_Tyrhs>& rhs)
{	//断定左边矩阵的列数与右边矩阵的行数相等
	Assert(lhs.GetColNum() == rhs.GetRowNum());
	//生成矩阵新对象，用lhs的行作为新阵的行数，用rhs的列数作为新阵的列数
	matrix<_Tyout> mTmp(lhs.GetRowNum(), rhs.GetColNum());

	for(size_t i = 0; i < mTmp.GetRowNum(); i ++)
	{
		for(size_t j = 0; j < mTmp.GetColNum(); j ++)
		{
			mTmp(i, j) = _Tyout(0);		//赋初值0
			for(size_t k = 0; k < lhs.GetColNum(); k ++)
			{
				mTmp(i, j) += lhs(i, k) * rhs(k, j);
			}
		}
	}

	mOut = mTmp;		//将最后结果转放入mOut矩阵中
	return mOut;		//返回结果矩阵mOut
}

//矩阵乘以向量
template <class _Tyout,class _Tylhs,class _Tyrhs>
valarray<_Tyout>& MatrixMultiplyVector(valarray<_Tyout>& mOut,matrix<_Tylhs>& lhs,valarray<_Tyrhs>& rhs)
{
	//断定左边矩阵的列数与右边向量的行数相等
	Assert(lhs.GetColNum() == rhs.size());

	valarray<_Tyout> mTmp(rhs.size());

	for(size_t i = 0; i < lhs.GetRowNum(); i ++)
	{
		mTmp[i] = _Tyout(0);		//赋初值0
		for(size_t j = 0; j <lhs.GetColNum(); j ++)
		{
			mTmp[i] += lhs(i, j) * rhs[j];

		}
	}

	mOut = mTmp;		//将最后结果转放入mOut向量中
	return mOut;		//返回结果向量mOut
}

//输出矩阵函数		按一行一行进行输出
template <class _Ty>
void MatrixLinePrint(const matrix<_Ty>& mOut)
{	
	size_t sR, sC;
	sR=mOut.GetRowNum();
	sC=mOut.GetColNum();

	for(size_t stR=0; stR<mOut.GetRowNum(); stR++)
	{
		for(size_t stC=0; stC<mOut.GetColNum(); stC++)
		{
			cout.width(15);				//元素对齐，让每个元素占15列
			cout << mOut(stR, stC) << ' ';
		}
		cout << endl;
	}
}

//输出矩阵函数		按指定行进行输出
template <class _Ty>
void MatrixLinePrint(const matrix<_Ty>& mOut, size_t LineNo)
{	
	size_t sR, sC;

	sR=mOut.GetRowNum();
	sC=mOut.GetColNum();
	
	for(size_t stC=0; stC<mOut.GetColNum(); stC++)
	{
		cout.width(15);				//元素对齐，让每个元素占15列
		cout << mOut(LineNo, stC) << ' ';
	}
	cout << endl;
}

//矩阵转置	==	原阵在mIn，转置后的矩阵在mOut  ==
template <class _Ty>
void MatrixTranspose(const matrix<_Ty>& mIn,matrix<_Ty>& mOut)
{
	size_t sR, sC;
	sR = mIn.GetRowNum();			//取原矩阵行数
	sC = mIn.GetColNum();			//取原矩阵列数

	matrix<_Ty> mTemp(sC, sR);		//生成一新阵，行数与列数与原阵互换
	
	for(size_t stC=0; stC<sC; stC++)
		for(size_t stR=0; stR<sR; stR++)
			mTemp(stC, stR) = mIn(stR, stC);	//对新阵赋值
					
	mOut= mTemp;//返回新的转置阵
}

//判断矩阵对称
template <class _Ty>	
bool MatrixSymmetry(const matrix<_Ty>& rhs)
{
	bool bSy = true;
	size_t stRow = rhs.GetRowNum();	//取矩阵行数

	if(rhs.GetColNum() == stRow)	// 必须是方阵
	{
		for(size_t i = 1; i < stRow; i ++)			//判断是否对称
			for(size_t j = 0; j < i; j ++)
				if(FloatNotEqual((long double)rhs(i, j), (long double)rhs(j, i)))
				{
					bSy =  false;
					goto RET;
				}
	}
	else
		bSy = false;
	RET: return bSy;	//矩阵对称
}

//判断矩阵对称正定	
template <class _Ty>	
int MatrixSymmetryRegular(const matrix<_Ty>& rhs, int sym)
{							
	long double ldDet;
	size_t i, j, k;

	size_t sC = rhs.GetColNum();	//矩阵列数
	size_t sR = rhs.GetRowNum();	//矩阵行数

	size_t stRank = sR;				// 矩阵阶数
	if(stRank != rhs.GetRowNum())
		return int(-1);				// 不是方阵

	if(sym > 0)
		if(MatrixSymmetry(rhs)==false)
			return int(-2);			//rhs不是对称阵

	cout << " K = 1 \t Determinant = " << rhs(0, 0) <<endl;
	
	for(k = 0; k < stRank; k ++)	//若要判别半正定，负定，这句要修改
	{
		if(FloatEqual(rhs(k, k), 0)||rhs(k, k) < 0)
			return int(-3);	//对角元不大于0，矩阵不是正定阵
	}

	for(k = 2; k <= sR; k++)
	{
		matrix<long double> m(k, k);	//生成一matrix对象

		for(i=0; i<k; i++)
			for(j=0; j<k; j++)
				m(i, j) = (long double)rhs(i, j);	//初始化

		ldDet = MatrixDeterminant(m);				// 顺序主子式的值

		cout << " K = " << k << "\t Determinant = " << ldDet << endl; 
		
		if(FloatEqual(ldDet,0) || ldDet < 0.0)
			return (0);					//不是正定阵
	}
	if(sym == 1) return int(2);			//矩阵为正定对称阵
	else		 return int(1);			//矩阵为正定阵
}

//求向量范数
//t 是向量，NormType是范数类型
//NormType=1时，表示1-范数或列范数；NormType=2时，表示2-范数或Euclid范数，为默认值；NormType=-1时，表示∞-范数；
template <class _Ty>
double VectorNorm(const valarray<_Ty>& t,int NormType=2)
{
	double Norm;
	size_t n;

	//初始化
	Norm=0;
	n=t.size();

	switch(NormType)
	{
		case 1://1-范数
			{
				for(int i=0;i<n;i++)
					Norm+=fabs(t[i]);
				
				break;
			}
		case 2://2-范数
			{
				for(int i=0;i<n;i++)
					Norm+=fabs(t[i])*fabs(t[i]);				
			
				Norm=sqrt(Norm);				
				break;			
			}
		case -1://∞-范数
			{
				for(int i=0;i<n;i++)
				{
					if(fabs(t[i])>Norm)
						Norm=fabs(t[i]);
				}				
				break;			
			}	
	}
	return Norm;
}

//求矩阵范数
//rhs 是矩阵，NormType是范数类型
//NormType=0时，表示Frobenius范数，为默认值；
//NormType=1时，表示1-范数；NormType=2时，表示2-范数；NormType=-1时，表示∞-范数；
template <class _Ty>
double MatrixNorm(const matrix<_Ty>& rhs,int NormType=0)
{
	double	Norm;
	size_t	m_stRow;		//矩阵行数变量
	size_t	m_stCol;		//矩阵列数变量

	//初始化
	Norm	=0;
	m_stRow	=rhs.GetRowNum();
	m_stCol	=rhs.GetColNum();

	switch(NormType)
	{
		case 0:
			{
				double NormTemp=0;
				for(int i=0;i<m_stRow;i++)
				{
					for(int j=0;j<m_stCol;j++)
					{
						NormTemp+=rhs(i,j)*rhs(i,j);
					}
				}				
				Norm=sqrt(NormTemp);
				break;			
			}
		case 1://1-范数
			{
				double NormTemp;

				for(int j=0;j<m_stCol;j++)
				{
					NormTemp=0;

					for(int i=0;i<m_stRow;i++)
						NormTemp+=fabs(rhs(i,j));
					
					if(NormTemp>Norm)
						Norm=NormTemp;				
				}
				break;			
			}

		case 2://2-范数
			{
				matrix<double> mOut(m_stCol,m_stCol);
				matrix<double> rhsTranspose(m_stCol,m_stRow);
				MatrixTranspose(rhs,rhsTranspose);
				MatrixMultiply(mOut, rhsTranspose,rhs);//mOut=rhs'*rhs;
		
				//-----------QR法求特征值
				HessenbergTransform(mOut);	
				valarray<complex<double> > uv(m_stCol);
				EigenvalueVectorHessenbergQR(mOut,uv,1e-8,1e6);
				
				complex<double> EigMax=uv[0];
				
				for(int i=1;i<m_stCol;i++)
				{
					if(uv[i]>EigMax)
						EigMax=uv[i];				
				}
				Norm=abs(EigMax);
				
				break;			
			}

		case -1://∞-范数
			{
				double NormTemp;

				for(int i=0;i<m_stRow;i++)
				{
					NormTemp=0;

					for(int j=0;j<m_stCol;j++)
						NormTemp+=fabs(rhs(i,j));
					
					if(NormTemp>Norm)
						Norm=NormTemp;				
				}
				break;			
			}	
	}

	return Norm;
}


//全选主元法求矩阵行列式函数
template <class _Ty>
long double MatrixDeterminant(const matrix<_Ty>& rhs)		
{	
	long double  MaxValue, tmp;
	size_t stRank = rhs.GetColNum();	// 矩阵阶数
	if(stRank != rhs.GetRowNum())
		return long double(0);			//rhs不是方阵

	matrix<long double> m(stRank, stRank);	//生成一matrix对象

	for(size_t i=0; i<stRank; i++)
		for(size_t j=0; j<stRank; j++)
			m(i, j) = (long double)rhs(i, j);	//初始化

	size_t iSign, jSign;				// 主元的位置标志

	long double Det(1);					// 行列式的值

	int	nSgn = 1;						// 符号

	for(size_t k = 0; k < stRank - 1; k ++)	// 全选主元
	{	
		MaxValue = 0.0;
		for(i = k; i < stRank; i ++)
		{
			for(size_t j = k; j < stRank; j ++)
			{
				tmp = Abs(m(i, j));		//求m(i,j)绝对值
				if(tmp > MaxValue)
				{
					MaxValue = tmp;
					iSign = i;			//记下主元位置
					jSign = j;
				}
			}
		}
		
		if(FloatEqual(MaxValue, 0))	//绝对值最大元素为0，行列式也为0
			return long double(0);

		if(iSign != k)			//绝对值最大元素不在当前行
		{
			nSgn = -nSgn;		//变换行列式符号
			for(size_t j = k; j < stRank; j ++)		//交换行
				swap(m(k, j), m(iSign, j));
		}

		if(jSign != k)			//绝对值最大元素不在当前列
		{
			nSgn = -nSgn;		//变换行列式符号
			for(size_t i = k; i < stRank; i ++)		//交换列
				swap(m(i, jSign), m(i, k));
		}

		Det *= m(k, k);					//对角元相乘
		for(i = k + 1; i < stRank; i ++)
		{
			long double d(m(i, k) / m(k, k));		//消元因子
			for(size_t j = k + 1; j < stRank; j ++)	//将主元下方元素消为0
				m(i, j) -= d * m(k, j);	//当前主元行下行其余元素作变换
		}
	}

	return Det * nSgn * m(stRank - 1, stRank - 1);	//返回行列式数值
}

//全选主元高斯(Gauss)消去法求一般矩阵的秩
template <class _Ty>						//返回值为秩数
size_t MatrixRank(const matrix<_Ty>& rhs)
{
	size_t iSign, jSign;					//主元的位置标志
	size_t mRank = 0;						//矩阵秩数
	size_t stRow = rhs.GetRowNum();			//取矩阵行数
	size_t stCol = rhs.GetColNum();			//取矩阵列数
	size_t ColRowMin = Min(stRow, stCol);	//取行或列最小值

	matrix<_Ty> m(rhs);				//生成一matrix对象，用rhs初始化

	for(size_t k = 0; k < ColRowMin; k ++)
	{	// 全选主元
		long double MaxValue(0);
		for(size_t i = k; i < stRow; i ++)
		{
			for(size_t j = k; j < stCol; j ++)
			{
				long double tmp(Abs(m(i, j)));	//求m(i,j)绝对值
				if(tmp > MaxValue)
				{
					MaxValue = tmp;
					iSign = i;					//记下主元位置
					jSign = j;
				}
			}
		}
		
		//子阵元素绝对值最大者为0，	注意浮点数与0相等的定义，见comm.h
		if(FloatEqual(MaxValue, 0))