matrix.h

矩阵计算的c++代码

}
CMatrix& CMatrix::operator-=(__MatrixType__ right)
{	
	if(IsEmpty())
	{// if the matrix is null, return a null matrix
		return (*this);
	}
	for(int i=0; i<m_nRows; ++i)
	{
		for(int j=0; j<m_nCols; ++j)
		{	*(*(m_pData+i)+j) -= right;	}
	}
	return (*this);
}
CMatrix& CMatrix::operator-=(const CMatrix& right)
{	
	assert(right.m_nRows==m_nRows && right.m_nCols==m_nCols);

	if(IsEmpty())
	{// if the matrix is null, return a null matrix
		return (*this);
	}
	for(int i=0; i<m_nRows; ++i)
	{
		for(int j=0; j<m_nCols; ++j)
		{	*(*(m_pData+i)+j) -= *(*(right.m_pData+i)+j);	}
	}
	return (*this);
}


// operator *=
CMatrix& CMatrix::operator*=(__MatrixType__ right)
{
	if(IsEmpty())
	{// if the matrix is null, return a null matrix
		return (*this);
	}

	for(int i=0; i<m_nRows; ++i)
	{
		for(int j=0; j<m_nCols; ++j)
		{	*(*(m_pData+i)+j) *= right;	}
	}
	return (*this);
}

// get Hilbert matrix, a famous bad matrix
CMatrix& CMatrix::Hilbert(int nDim)
{
	SetSize(nDim, nDim);

	for(int i=0; i<m_nRows; ++i)
	{
		for(int j=0; j<m_nCols; ++j)
		{	*(*(m_pData+i)+j) = 1.0/(i+j+1);	}
	}
	return (*this);
}

//eye, unit matrix
CMatrix& CMatrix::Eye(void)
{
	if(IsEmpty())
	{// if the matrix is null, return a null matrix
		return (*this);
	}

	assert(m_nRows == m_nCols);

	Zeros();	// clear zero
	for(int i=0; i<m_nRows; ++i)
	{	*(*(m_pData+i)+i) = 1.0;	}
	return (*this);
}
CMatrix& CMatrix::Eye(int nDim)
{
	Zeros(nDim, nDim);	// clear zero
	for(int i=0; i<m_nRows; ++i)
	{	*(*(m_pData+i)+i) = 1.0;	}
	return (*this);
}
//ones, all matrix elements is one
CMatrix& CMatrix::Ones(void)
{
	if(IsEmpty())
	{// if the matrix is null, return a null matrix
		return (*this);
	}

	for(int i=0; i<m_nRows; ++i)
	{
		for(int j=0; j<m_nCols; ++j)
		{	*(*(m_pData+i)+j) = 1.0;	}
	}

	return (*this);
}
CMatrix CMatrix::Ones(int nRows, int nCols)
{
	SetSize(nRows, nCols);
	return Ones();
}
//zeros, zero matrix
CMatrix& CMatrix::Zeros(void)
{
	if(IsEmpty())
	{// if the matrix is null, return a null matrix
		return (*this);
	}

	for(int i=0; i<m_nRows; ++i)
	{
		for(int j=0; j<m_nCols; ++j)
		{	*(*(m_pData+i)+j) = 0.0;	}
	}

	return (*this);
}
CMatrix CMatrix::Zeros(int nRows, int nCols)
{
	SetSize(nRows, nCols);
	return Zeros();
}

//triangle lower
CMatrix CMatrix::TriL(bool bDiagonal)const
{	// bDiagonal==true 表示取对角线上的元素，否则对角线上的元素为0
	if(IsEmpty())
	{	return (*this);	}
	assert(m_nRows == m_nCols);

	int i, j, kk;
	CMatrix _result;
	_result.Zeros(m_nRows, m_nCols);

	if(bDiagonal)	{	kk = 1;	}
	else			{	kk = 0;	}
	for(i=0; i<m_nRows; ++i)
	{
		for(j=0; j<i+kk; ++j)
		{	*(*(_result.m_pData+i)+j) = *(*(m_pData+i)+j);	}
	}

	return _result;
}
//triangle upper
CMatrix CMatrix::TriU(bool bDiagonal)const
{	// bDiagonal==true 表示取对角线上的元素，否则对角线上的元素为0
	if(IsEmpty())
	{	return (*this);	}
	assert(m_nRows == m_nCols);

	int i, j, kk;
	CMatrix _result;
	_result.Zeros(m_nRows, m_nCols);

	if(bDiagonal)	{	kk = 1;	}
	else			{	kk = 0;	}

	for(i=0; i<m_nRows; ++i)
	{
		for(j=m_nCols-1; j>i-kk; --j)
		{	*(*(_result.m_pData+i)+j) = *(*(m_pData+i)+j);	}
	}
	return _result;
}

int CMatrix::LuDcmp(CMatrix &lower, CMatrix &upper, vector<int>& index) const	// make LU decomposition
{	// 列主元LU分解
	assert(m_nRows == m_nCols);
	if(m_nRows == 0)	return -3;

	int i, imax, j, k, det;
	__MatrixType__ big, dum, sum, temp;
	CVector vv(m_nRows);
	CMatrix AA(*this);	// AA的作用：由于在lu分解过程中改变了原来矩阵的内容，这里用一个变量代替使得原来的矩阵不改变

	// 进行矩阵的行初等变换时，记录行列式值的正负号
	det = 1;	// 初始值设为 1
	upper.Zeros();	// clear zero
	lower.Zeros();

	for(i=0; i<m_nRows; ++i)
	{	// 选出矩阵中每行绝对值最大的元素，也就是缩放信息
		big = 0.0;
		for (j=0; j<m_nCols; ++j)
		{
			temp = fabs( *(*(AA.m_pData+i)+j) );
			if ( temp > big )	big = temp;
		}
		if (big == 0.0)
		{	cout<<"Singular Matrix in Routine LUDCMP"<<endl;	return(-4);	}
		vv[i] = 1.0/big;	// 记录其倒数
	}

	for (j=0; j<m_nCols; ++j)
	{	// 开始分解
		for (i=0; i<j; ++i)
		{	// 计算临时值sum，并冲掉AA(i,j)
			sum = *(*(AA.m_pData+i)+j);
			for(k=0; k<i; ++k)
			{	sum -= *(*(AA.m_pData+i)+k)   *    *(*(AA.m_pData+k)+j);	}
			*(*(AA.m_pData+i)+j) = sum;
		}
		big = 0.0;
		for (i=j; i<m_nRows; ++i)
		{	// 计算临时值sum，并冲掉AA(i,j)
			sum = *(*(AA.m_pData+i)+j);
			for (k=0; k<j; ++k)
			{	sum -= *(*(AA.m_pData+i)+k)   *    *(*(AA.m_pData+k)+j);	}
			*(*(AA.m_pData+i)+j) = sum;
			if ( (dum=vv[i]*fabs(sum)) >= big)	// 选主元，记录行号
			{	big = dum;	imax = i;	}
		}
		if(j != imax)	// 是否需要换行？
		{	// 换行、记录行号
			for (k=0; k<m_nCols; ++k)
			{
				dum = *(*(AA.m_pData+imax)+k);
				*(*(AA.m_pData+imax)+k) = *(*(AA.m_pData+j)+k);
				*(*(AA.m_pData+j)+k) = dum;
			}
			det = -det;	// 行列式符号改变
			vv[imax] = vv[j];
		}
		index[j] = imax;	// 记录变换矩阵

		// 如果出现这种情况，说明这个矩阵是奇异矩阵。
		// 但为了使得计算能够进行下去，采取了这种变通方式。
		if(*(*(AA.m_pData+j)+j) == 0.0)	{	*(*(AA.m_pData+j)+j) = ERROR20;	}
		// 也可以采取如下方式
	//	if(*(*(AA.m_pData+j)+j) == 0.0)	{	cout<<"Singular Matrix in Routine LUDCMP"<<endl;	return (-4);	}
		
		if(j != m_nRows-1)
		{
			dum = 1.0 / *(*(AA.m_pData+j)+j);
			for(i=j+1; i<m_nRows; ++i)
			{	*(*(AA.m_pData+i)+j) *= dum;	}
		}
	}
	upper = AA.TriU(true);
	lower = AA.TriL(false);
	for(i=0; i<m_nRows; ++i)
	{	lower(i,i) = 1.0;	}

	return det;
}

int CMatrix::LuBksb(CVector & xx, const CVector &bb)const	// solve a linear equations
{	// LU BackSubstitution
	assert(m_nRows==m_nCols && m_nRows==bb.size() && xx.size()>=bb.size());
	if(m_nRows == 0)	return -3;

	int det, i, j, ip;
	__MatrixType__ tSum;
	vector<int> index(m_nRows);
	CMatrix upper(m_nRows, m_nCols), lower(m_nRows, m_nCols);

	xx = bb;	// 初始化xx为bb, xx=bb，这样做的目的是不改变bb中的值
	// 进行 LU 分解
	det = LuDcmp(lower, upper, index);

	// sovle L*y = P*b
	ip = index[0];
	if(ip != 0)
	{	// 给出变换后的yy的第一个值
		tSum = xx[ip];
		xx[ip] = xx[0];
		xx[0] = tSum;
	}
	for(i=1; i<m_nRows; ++i)
	{
		ip = index[i];
		tSum = xx[ip];
		xx[ip] = xx[i];
		for(j=0; j<i; ++j)
		{	tSum -= lower(i,j) * xx[j];	}
		xx[i] = tSum;
	}

	// sovle U*x = y
	for(i=m_nRows-1; i>=0; --i)
	{
		tSum = xx[i];
		for(j=m_nRows-1; j>i; --j)
		{	tSum -= upper(i,j) * xx[j];	}
		xx[i] = tSum / upper(i,i);
	}
	return det;
}
int CMatrix::LuBksb(DVector & xx, const DVector &bb)const	// solve a linear equations
{	// LU BackSubstitution
	int det;
	CVector _xx(xx), _bb(bb);
	det = LuBksb(_xx, _bb);
	for(int i=0; i<_xx.size(); ++i)
	{	xx[i] = _xx[i];	}
	return det;
}

CMatrix CMatrix::Inverse(void)const	// get inverse matrix
{
	if(IsEmpty())
	{	return (*this);	}
	assert(m_nRows==m_nCols);

	int i, j, det;
	vector<int> index(m_nRows);
	CMatrix upperInverse(m_nRows, m_nCols), upper(m_nRows, m_nCols), lowerInverse(m_nRows, m_nCols),
		lower(m_nRows, m_nCols), trans(m_nRows, m_nCols), II(m_nRows, m_nCols);
	CVector xUpper(m_nRows), xLower(m_nRows), col(m_nRows);

	det = LuDcmp(lower, upper, index);	// 进行 LU 分解
	upperInverse.Zeros();
	lowerInverse.Zeros();
	trans.Eye();	// 矩阵单位化

	for(j=0; j<m_nCols; ++j)
	{
		for(i=0; i<m_nRows; ++i)
		{	col[i] = 0.0;	}
		col[j] = 1.0;
		lower.LuBksb(xLower, col);
		upper.LuBksb(xUpper, col);
		for(i=0; i<m_nRows; ++i)
		{
			lowerInverse(i,j) = xLower[i];
			upperInverse(i,j) = xUpper[i];
		}
	}
	for(i=0; i<m_nRows; ++i)
	{
		if(i != index[i])
		{
			II.Eye();	// 矩阵单位化
			II(i,i) = 0.0;	II(index[i],index[i]) = 0.0;
			II(i,index[i]) = 1.0;	II(index[i],i) = 1.0;
			trans = II * trans;
		}
	}
	return (upperInverse * lowerInverse * trans);	
}

// return Determinant
__MatrixType__ CMatrix::Det(void)const
{
	if(IsEmpty())
	{	return 0.0;	}

	int i, det;
	__MatrixType__ _result=1.0;
	vector<int> index(m_nRows);
	CMatrix upper(m_nRows, m_nCols), lower(m_nRows, m_nCols);

	det = LuDcmp(lower, upper, index);

	for(i=0; i<m_nRows; ++i)
	{	_result *= upper(i,i);	}

	return det*_result;
}

// the condition number of the matrix
__MatrixType__ CMatrix::Cond(void)const
{
	if(IsEmpty())
	{	return 0;	}

	return Norm() * (Inverse().Norm());
}

// solve a linear equations using sor iterative method
int CMatrix::Sor(CVector & xx, const CVector &bb, double omega, bool bInitlized, int nMaxIter, double dError)const
{	// SorSolver() 使用逐次超松弛迭代方法求解线性方程组。sor means successive over-relaxation
	// omega(ω)是松弛因子，0<ω<2 是迭代收敛的必要条件。
	// 当1<ω<2时，是超松弛迭代，当0<ω<1时，是欠松弛迭代，当ω=1时，是Gauss-Seidel迭代。ω初始值为1
	// nMaxIter 是迭代的最大次数，nMaxIter的值不宜太小。初始值为100
	// bInitlized 表示是否使用已经初始化的初始向量xx来迭代

	assert(omega>0.0 && omega<2.0);	// 0<ω<2 是迭代收敛的必要条件
	assert(nMaxIter >= 50);	// 迭代的最大次数不宜太小
	assert(m_nRows==m_nCols && m_nRows==bb.size() && xx.size()>=bb.size());	//确定xx不比bb小

	if(m_nRows == 0)	return (-3);	// 如果方程数目为 0，则返回0

	const int DIM = bb.size();	// 取得未知数的个数
	int nTime, i, j;
	double norm=1.0;
	CVector xx1(DIM);
	CMatrix AA(*this);
	
	if(bInitlized)	{	xx1 = xx;	}	// bInitlized==true 用xx来迭代
	else	{	xx1 = xx = 1.0;	}	// 否则，初始化xx1,xx为1.0

	for(nTime=0; nTime<nMaxIter; ++nTime)
	{	// 
		if( norm >= dError && nTime>=nMaxIter-1 )	// 如果不收敛，就退出。
		{	cout<<"The iteration does not converge."<<endl;	return (-4);}
		for(i=0; i<DIM; ++i)
		{	//	迭代开始
			xx1[i] = bb[i];
			for(j=0; j<i; ++j)	// 减去(n+1)部分
			{	xx1[i] -= AA(i,j)*xx1[j];	}
			for(j=i; j<DIM; ++j)	// 减去(n)部分
			{	xx1[i] -= AA(i,j)*xx[j];	}
			xx1[i] = xx[i] + omega*xx1[i]/AA(i,i);	// 乘上松弛因子
		}
		norm = (xx-xx1).Norm() / xx.Norm();
		xx = xx1;
		if( norm < dError )  break;
	}
	return ( nTime );
}

int CMatrix::Sor(DVector & xx, const DVector &bb, double omega, bool bInitlized, int nMaxIter, double dError)const	// solve a linear equations using sor method
{
	int nTime;
	CVector _xx(xx), _bb(bb);
	nTime = Sor(_xx, _bb, omega, bInitlized, nMaxIter, dError);
	for(int i=0; i<_xx.size(); ++i)
	{	xx[i] = _xx[i];	}
	return nTime;
}

	// solve a linear equations using Jacobi iterative method
int CMatrix::ItJacobi(CVector & xx, const CVector &bb, bool bInitlized, int nMaxIter, double dError)const
{	// ItJacobi() 使用Jacobi迭代方法求解线性方程组。
	// nMaxIter 是迭代的最大次数，nMaxIter的值不宜太小。初始值为100
	// bInitlized 表示是否使用已经初始化的初始向量xx来迭代

	assert(nMaxIter >= 50);	// 迭代的最大次数不宜太小
	assert(m_nRows==m_nCols && m_nRows==bb.size() && xx.size()>=bb.size());	//确定xx不比bb小

	if(m_nRows == 0)	return (-3);	// 如果方程数目为 0，则返回0

	const int DIM = bb.size();	// 取得未知数的个数
	int nTime, i, j;
	double norm=1.0;
	CVector xx1(DIM);
	CMatrix AA(*this);
	
	if(bInitlized)	{	xx1 = xx;	}	// bInitlized==true 用xx来迭代
	else	{	xx1 = xx = 1.0;	}	// 否则，初始化xx1,xx为1.0

	for(nTime=0; nTime<nMaxIter; ++nTime)
	{	// 
		if( norm >= dError && nTime>=nMaxIter-1 )	// 如果不收敛，就退出。
		{	cout<<"The iteration does not converge."<<endl;	return (-4);}
		for(i=0; i<DIM; ++i)
		{	// 减去两部分，不包括i
			xx1[i] = bb[i];
			for(j=0; j<i; ++j)	// 小于i的部分
			{	xx1[i] -= AA(i,j)*xx[j];	}
			for(j=i+1; j<DIM; ++j)	// 小于i的部分
			{	xx1[i] -= AA(i,j)*xx[j];	}
			xx1[i] /= AA(i,i);
		}
		norm = (xx-xx1).Norm() / xx.Norm();
		xx = xx1;
		if( norm < dError )  break;
	}
	return ( nTime );
}

// define some additional friend operations of matrix
#include "xxMatrix"

// namespace end
__FMXL_END_NAMESPACE

#endif //__MATRIX_H