matrix.inl

本程序是在CCS环境下运行的

// 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 _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(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);			//矩阵为正定阵
}

//全选主元法求矩阵行列式函数
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)) 
			break;	//return mRank;
		else				
			mRank++;		//子阵元素绝对值最大者不为0，矩阵秩加1

		if(k ==(ColRowMin - 1))		//已到最后一行(列)
			break;	//return mRank;
		
		if(iSign != k)				//主元不在当前行
		{
			for(size_t j = k; j < stCol; j ++)	//交换行元素
				swap(m(k, j), m(iSign, j));
		}

		if(jSign != k)				//主元不在当前列
		{
			for(size_t i = k; i < stRow; i ++)	//交换列元素
				swap(m(i, jSign), m(i, k));
		}

		for(i = k + 1; i < stRow; i ++)
		{
			const _Ty d(m(i, k) / m(k, k));		//消元因子
			for(size_t j = k + 1; j < stCol; j ++)	
				m(i, j) -= d * m(k, j);		//当前主元右下阵元素作变换
		}
	}
	return mRank;
}

//全选主元高斯-约当(Gauss-Jordan)法求矩阵逆
template <class _Ty>
int MatrixInversionGS(matrix<_Ty >& rhs)
{
	size_t stRank = rhs.GetColNum();	// 矩阵阶数
	if(stRank != rhs.GetRowNum())
		return int(-1);					//rhs不是方阵

	valarray<size_t> is(stRank);		//行交换信息
	valarray<size_t> js(stRank);		//列交换信息
	
	matrix<_Ty> m(rhs);					//生成一matrix对象

	for(size_t k = 0; k < stRank; k++)
	{	// 全选主元
		long double MaxValue(0);
		for(size_t i = k; i < stRank; i ++)
		{
			for(size_t j = k; j < stRank; j ++)
			{
				long double tmp(Abs(m(i, j)));	//求m(i,j)绝对值
				if(tmp > MaxValue)				//主元不在对角线上
				{
					MaxValue = tmp;
					is[k] = i;					//记录主元行数
					js[k] = j;					//记录主元列数
				}
			}
		}
		
		if(FloatEqual(MaxValue, 0)) 
			return int(0);						//主元为0，矩阵奇异
		
		if(is[k] != k)							//主元不在当前行
		{
			for(size_t j = 0; j < stRank; j ++)	//交换行元素
				swap(m(k, j), m(is[k], j));
		}

		if(js[k] != k)							//主元不在当前列
		{
			for(size_t i = 0; i < stRank; i ++)	//交换列元素
				swap(m(i, k), m(i, js[k]));
		}

		m(k, k) = 1.0 / m(k, k);				//主元倒数
		for(size_t j = 0; j < stRank; j ++)
			if(j != k)
				m(k, j) *= m(k, k);

		for(i = 0; i < stRank; i ++)
			if(i != k)
				for(size_t j = 0; j < stRank; j ++)	
					if(j != k)
						m(i, j) = m(i, j) - m(i, k) * m(k, j);

		for(i = 0; i < stRank; i ++)
			if(i != k)
				m(i, k) = -m(i, k) * m(k, k);
	}

	for(int r = stRank - 1; r >= 0; r--)
	{
		if(js[r] != r)
			for(size_t j = 0; j < stRank; j ++)
				swap(m(r, j), m(js[r], j));
		if(is[r] != r)
			for(size_t i = 0; i < stRank; i ++)
				swap(m(i, r), m(i, is[r]));
	}

	rhs = m;

	return int(1);
}

//用“变量循环重新编号法”求对称正定矩阵逆
//矩阵类型必须是浮点型
template <class _Ty>
int MatrixSymmetryRegularInversion(matrix<_Ty>& rhs)
{
	int iSym = MatrixSymmetryRegular(rhs, 1);	//判别是否对称正定
	if(iSym < 2)
		return (iSym);			//rhs不是对称正定阵

	size_t stRank = rhs.GetColNum();	// 矩阵阶数
	
	matrix<_Ty> msr(rhs);			//生成一matrix对象，用rhs初始化

	valarray<_Ty> b(stRank);
	
	for(size_t k=0; k<stRank; k++)
	{
		_Ty w= msr(0, 0);
		size_t m = stRank - k -1;
		for(size_t i = 1; i < stRank; i++)
		{
			_Ty g = msr(i, 0);
			b[i] = g / w;
			if(i <= m) b[i] = -b[i];
			for(size_t j = 1; j <= i; j ++)
				msr((i-1),(j-1)) = msr(i, j) + g * b[j];
		}
		msr(stRank-1, stRank-1) = 1.0 / w;
		for(i = 1; i < stRank; i ++)
			msr(stRank-1,(i-1)) =  b[i];
	}
	for(size_t i = 0; i < stRank-1; i ++)
		for(size_t j = i+1; j < stRank; j ++)
			msr(i,j) = msr(j, i);

	rhs = msr;

	return (iSym);
}


//特兰持(Trench)法求托伯利兹(Toeplitz)矩阵逆
//矩阵类型必须是浮点型
template <class _Ty>
int MatrixToeplitzInversionTrench(const valarray<_Ty>& t, const valarray<_Ty>& tuo, matrix<_Ty>& rhs)
{
	size_t stRank = rhs.GetColNum();	// 矩阵阶数
	if(stRank != rhs.GetRowNum())
		return int(-1);					//rhs不是方阵

	if(FloatEqual(t[0], 0))
		return int(0);	

	valarray<_Ty> c(stRank);
	valarray<_Ty> r(stRank);
	valarray<_Ty> p(stRank);

	_Ty a=t[0]; 
	c[0]=tuo[1]/t[0];
	r[0]=t[1]/t[0];

	matrix<_Ty> b(rhs);
	for(size_t k=0; k<=stRank-3; k++)
	{
		_Ty s=0.0;
		for(size_t j=1; j<=k+1; j++)
			s=s+c[k+1-j]*tuo[j];
		s=(s-tuo[k+2])/a;
		for(size_t i=0; i<=k; i++)
			p[i]=c[i]+s*r[k-i];
        c[k+1]=-s;
        s=0.0;
        for(j=1; j<=k+1; j++)
			s=s+r[k+1-j]*t[j];
        s=(s-t[k+2])/a;
        for(i=0; i<=k; i++)
        {
			r[i]=r[i]+s*c[k-i];
            c[k-i]=p[k-i];
        }
        r[k+1]=-s;
		a=0.0;
        for(j=1; j<=k+2; j++)
			a=a+t[j]*c[j-1];
        a=t[0]-a;
		if(FloatEqual(a, 0))
			return int(0);	
    }
    b(0,0)=1.0/a;
    for(size_t i=0; i<stRank-1; i++)
    {	
		k=i+1;
        b(0, k)=-r[i]/a;
		b(i+1,0)=-c[i]/a;
    }
    for(i=0; i<stRank-1; i++)
    for(size_t j=0; j<stRank-1; j++)
        b(i+1, j+1)=b(i,j)-c[i]*b(0,j+1)+c[stRank-j-2]*b(0,stRank-i-1);

	rhs = b;
    
    return int(1);
}

//实矩阵LU分解
//矩阵必须是浮点型
template <class _Ty>
int MatrixLU(const matrix<_Ty>& rhs, matrix<_Ty>& lhs, matrix<_Ty>& uhs)
{
	size_t stRank = rhs.GetColNum();	// 矩阵阶数
	if(stRank != rhs.GetRowNum())
		return int(-1);			//rhs不是方阵
	
	matrix<_Ty> m(rhs);
	for(size_t k=0; k<stRank-1; k++)
	{
		if(FloatEqual(m(k,k),0))
			return int(0);		//主元为0
		for(size_t i=k+1; i<stRank; i++)
			m(i,k) /= m(k,k);
		for(i=k+1; i<stRank; i++)
			for(size_t j=k+1; j<stRank; j++)
				m(i,j)=m(i,j)-m(i,k)*m(k,j);
	}
	//给上、下三角阵赋值
	for(size_t i=0; i<stRank; i++)
	{
		for(size_t j=0; j<i; j++)
		{
			lhs(i,j)=m(i,j);
			uhs(i,j)=0.0;
		}
		lhs(i,i)=1.0;
		uhs(i,i)=m(i,i);
		for(j=i+1; j<stRank; j++)
		{
			lhs(i,j)=0.0;
			uhs(i,j)=m(i,j);
		}
	}
	
	return (1);		//分解成功
}

//用豪斯荷尔德(Householder)变换对一般m*n阶的实矩阵进行QR分解
template <class _Ty>
int MatrixQR(matrix<_Ty>& rhs, matrix<_Ty>& rhq)
{
	size_t stRow = rhs.GetRowNum();	// 矩阵行数
	size_t stCol = rhs.GetColNum();	// 矩阵列数
	
	if(stRow < stCol)
		return (0);					//行不能小于列

	for(size_t i=0; i<stRow; i++)
    for(size_t j=0; j<stRow; j++)
    { 
		rhq(i,j)=0.0;
        if(i==j) rhq(i,j)=1.0;
    }

	size_t nn=stCol;
    
	if(stRow == stCol) nn=stRow-1;

	for(size_t k=0; k<nn; k++)
    {
		_Ty u=0.0;

        for(size_t i = k; i < stRow; i++)
        { 
			_Ty w = Abs(rhs(i,k));
            if(w > u) u = w;
        }
        _Ty alpha=0.0;
        for(i = k; i < stRow; i++)