matrix.cpp

各种矩阵运算的实现

// CMATRIX.cpp: implementation of the CMATRIX class.
//
//////////////////////////////////////////////////////////////////////
#include "stdafx.h"
#include "MATRIX.h"
#include "math.h"

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
//#define new DEBUG_NEW
#endif
//////////////////////////////////////////////////////////////////////
int couter=0;

inline Radian ToRad(Angle angle)
{ return (angle/180)*WKPI; }

inline Angle ToAng(Radian rad)
{ return (rad/WKPI)*180; }

double ctan(Radian rad)
{ return 1/tan(rad); }

double csc(Radian rad)
{ return 1/sin(rad); }

Radian actan(double r)
{ return atan(1/r); }

double Cos(Radian r)
{
	if(((int)ToAng(r)/90)==ToAng(r)/90&&((int)ToAng(r)/90)%2!=0)
		return 0;
	else
		return cos(r);
}

double round(double num,unsigned n=0)//四舍五入
{
	double dex=num-(long)num;
	dex*=pow(10,n);
	double here=dex-(long)dex;
	if(here>=0.5)
	{
		dex=((long)dex+1)/pow(10,n);
		return (long)num+dex;
	}
	else
	{
		dex=(long)dex/pow(10,n);
		return (long)num+dex;
	}
}
//========================================
struct _EXPARA
{
	unsigned ROW;
	unsigned COL;
	unsigned r1,c1;
	unsigned r2,c2;
} EXPARA;

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CMATRIX::CMATRIX(unsigned r/*=1*/,unsigned c/*=1*/)
{
	if(Init()==false) return ;
//	ASSERT(Init());
		//无法分配堆内存以构建对象
	if(row<=0)	row=1;
		//指定的矩阵的行为零！将强制行为1！
	if(col<=0)	col=1;
		//指定的矩阵的列为零！将强制行为1！
	row=r; col=c;
/*
	couter++;
	cout<<"\nCON "<<couter<<endl;
//*/
}

CMATRIX::~CMATRIX()
{
	Destroy();

}

bool CMATRIX::Clear()
{
	Destroy();
	return Init();
}

bool CMATRIX::Init()
{
	if((head=fence=tail=new CElement)==NULL) return false;
	row=col=1; len=0; 
	return true;
}

void CMATRIX::Destroy()
{
	while(tail!=NULL)
	{
		fence=tail;
		tail=tail->prev;
		delete fence;
	}
/*
	couter--;
	cout<<"\nDE  "<<couter<<endl;
//*/
}

bool CMATRIX::Append(const double &_item_)
{
	if(tail==head&&len==0)
	{ fence->Set(_item_); len=1; }

	else
	{
		if(len==row*col) return true;
		else if((tail=tail->next=new CElement(_item_,tail,NULL))==NULL)
			return false;//无法分配堆内存以继续！
		else len++;
	}
	return true;
}

inline unsigned CMATRIX::RowColToLen(unsigned r, unsigned c)
{ return ((r-1)*col+c); }

inline double CMATRIX::GetCell(unsigned l)
{
	if(l>=1&&l<=len)
	{
		fence=head;
		if(l==1)
		{ return fence->Cell();}
		else
		{
			for(unsigned i=1;i<l;i++)
			{ fence=fence->next; }
			return fence->Cell();
		}
	}
	else return (0.0);//越界读取矩阵元素！
}

bool CMATRIX::SetCell(unsigned r, unsigned c, double it)
{
	if(r<1||r>row||c<1||c>col) return false;
		//越界设置矩阵元素！
	unsigned _l=RowColToLen(r,c);
	return SetCell(_l,it);
}

inline bool CMATRIX::SetCell(unsigned l, double cel)
{
	if(l>=1&&l<=len)
	{
		fence=head;
		if(l==1) { fence->Set(cel); return true; }
		for(unsigned i=1;i<l;i++)
		{ fence=fence->next; }
		fence->Set(cel);
		return true;
	}
	else return false;//越界设置矩阵元素！
}

inline void CMATRIX::SetRow(unsigned r)
{ row=r; }

inline void CMATRIX::SetCol(unsigned c)
{ col=c; }

inline unsigned CMATRIX::Row()
{ return row; }

inline unsigned CMATRIX::Col()
{ return col; }

inline unsigned CMATRIX::Len()
{ return len; }

bool CMATRIX::operator ==(CMATRIX &_m_)
{
	if((head==_m_.head&&tail==_m_.tail)||this==(&_m_)) return true;
	else
	{
		if(len==_m_.Len()&&row==_m_.Row()&&col==_m_.Col())
		{
			CElement* p;
			unsigned i=0;
			for(fence=head,p=_m_.head;fence,p;fence=fence->next,p=p->next)
			{ 
				if(fence->Cell()==p->Cell()) i++;
			}
			if(i==len) return true;
			else return false;
		}
		else return false;
	}
}

void CMATRIX::operator =(CMATRIX &_m_)
{
	if((*this)==_m_) return ;
	Clear();
	row=_m_.Row(); col=_m_.Col();
	for(unsigned i=1;i<=_m_.Len();i++)
	{ Append(_m_.GetCell(i)); }
}

double CMATRIX::operator ()(unsigned r,unsigned c/*=1*/)
{ return GetCell(RowColToLen(r,c)); }

inline double CMATRIX::Cell(unsigned int r, unsigned int c)
{
	return (*this)(r,c);
}

CMATRIX& operator +(CMATRIX &_m1_, CMATRIX &_m2_)
{
	TMPPMAT(_m1_.Row(),_m2_.Col());
	if(_m1_.Len()==0||_m2_.Len()==0)
	{
		//矩阵为空！加法无效！
		PMATSET(1,1);
		BACK;
	}
	if(_m1_.Row()!=_m2_.Row()||_m1_.Col()!=_m2_.Col()||_m1_.Len()!=_m2_.Len())
	{
		//行或列数不等，无法相加！
		BACK;
	}
	for(unsigned i=1;i<=_m2_.Len();i++)
	{
		if(!(pmat->Append(_m1_.GetCell(i)+_m2_.GetCell(i))))
			BACK;//无法分配堆内存以继续加法！
	}
	BACK;
}

CMATRIX& operator -(CMATRIX &_m1_, CMATRIX &_m2_)
{
	TMPPMAT(_m1_.Row(),_m2_.Col()); 
	if(_m1_.Len()==0||_m2_.Len()==0)
	{
		//矩阵为空！减法无效！
		PMATSET(1,1);
		BACK;
	}
	if(_m1_.Row()!=_m2_.Row()||_m1_.Col()!=_m2_.Col()||_m1_.Len()!=_m2_.Len())
	{
		//行或列数不等，无法相减！
		BACK;
	}
	for(unsigned i=1;i<=_m2_.Len();i++)
	{ 
		if(!(pmat->Append(_m1_.GetCell(i)-_m2_.GetCell(i))))
			BACK;//无法分配堆内存以继续减法！
	}
	BACK;
}

CMATRIX& operator *(CMATRIX &_m1_, CMATRIX &_m2_)
{
	TMPPMAT(_m1_.Row(),_m2_.Col());
	if(_m1_.Len()==0||_m2_.Len()==0)
	{
		//矩阵为空！乘法无效！
		PMATSET(1,1);
		BACK;
	}
	if(_m1_.Col()!=_m2_.Row())
	{
		//行列不匹配，无法相乘！
		BACK;
	}
	double tmpdat=0;
	for(unsigned i=1;i<=_m1_.Row();i++)
	{
		for(unsigned j=1;j<=_m2_.Col();j++)
		{
			for(unsigned k=1;k<=_m2_.Row();k++)
			{ tmpdat=tmpdat+_m1_(i,k)*_m2_(k,j); }
			if(!(pmat->Append(tmpdat)))
				BACK;//无法分配堆内存以继续两矩阵相乘！
			tmpdat=0;
		}
	}
	BACK;
}

CMATRIX& operator *(double k, CMATRIX &_m_)
{
	TMPPMAT(_m_.Row(),_m_.Col());
	if(_m_.Len()==0)
	{
		//矩阵为空！数乘无效！
		PMATSET(1,1);
		BACK;
	}
	for(unsigned i=1;i<=_m_.Len();i++)
	{
		if(!(pmat->Append(k*_m_.GetCell(i))))
		{
			//无法分配堆内存以继续数乘矩阵！
			BACK;
		}
	}
	BACK;
}

CMATRIX& CMATRIX::T()
{
	TMPPMAT(col,row);
	if(len==0)
	{
		//矩阵为空！无法转置！
		PMATSET(1,1);
		BACK;
	}
	for(unsigned i=1;i<=col;i++)
	{
		for(unsigned j=1;j<=row;j++)
		{
			if(!(pmat->Append((*this)(j,i))))
				BACK;//无法分配堆内存以继续矩阵转置！
		}
	}
	BACK;
}

CMATRIX& CMATRIX::Remmat(unsigned r, unsigned c)
{
	if(len==0)
	{
		//矩阵为空！求取余子式矩阵无效！
		return *this;
	}
	if(r>row||c>col||r<=0||c<=0)
	{
		//越界求取余子式矩阵！
		return *this;
	}
	TMPPMAT(row-1,col-1);
	if((row==1&&col==1)||pmat->Row()==0||pmat->Col()==0)
	{
		if(!(pmat->Append(0)))
		{
			//无法分配堆内存以继续求矩阵余子式！
			BACK;
		}
		PMATSET(row,col);
		BACK;
	}
	double tmpdat=-1;
	for(unsigned i=1;i<=row;i++)
	{
		if(i==r) { continue; }
		for(unsigned j=1;j<=col;j++)
		{
			if(j==c) { continue; }
			tmpdat=(*this)(i,j);
			if(!(pmat->Append(tmpdat)))
			{
				//无法分配堆内存以继续求矩阵余子式！
				BACK;
			}
		}
	}
	BACK;
}

CMATRIX& adjoint(CMATRIX &_m_)
{
	TMPPMAT(_m_.Row(),_m_.Col());
	if(_m_.Len()==0)
	{
		//矩阵为空！无法求其对应的伴随矩阵！
		PMATSET(1,1);
		BACK;
	}
	int modulus=-1;
	for(unsigned i=1;i<=_m_.Row();i++)
	{
		for(unsigned j=1;j<=_m_.Col();j++)
		{
			modulus=modulus*(-1);
			if(!(pmat->Append(modulus*det(_m_.Remmat(i,j)))))
			{
				//无法分配堆内存以继续求其对应的伴随矩阵！
				BACK;
			}
		}
	}
//	*pmat=pmat->T();
	pmat=&pmat->T();
	BACK;
}

double det(CMATRIX &_m_)
{
	if(_m_.Len()==0)
	{
		//矩阵为空！其对应的行列式值为零！
		return (0);
	}
	int cor=-1;
	double tmpdat=0;
	if(_m_.Row()!=_m_.Col())
	{
		//行列不同，无法求其对应的行列式值！
		return tmpdat;
	}
	if(_m_.Len()==1) { return _m_(1,1); } 
	if((_m_.Row()==1&&(_m_.Col()==1)||_m_.Len()==1)) { return _m_(1,1); }
	if((_m_.Row()==1&&(_m_.Col()==1)||_m_.Len()==1)) { return _m_(1,1); }
	for(unsigned j=1;j<=_m_.Col();j++)
	{
		cor=cor*(-1);
		tmpdat=tmpdat+_m_(1,j)*det(_m_.Remmat(1,j))*cor;
	}
	return tmpdat;
}

CMATRIX& CMATRIX::GetLine(unsigned lin)
{
	TMPPMAT(1,col);
	if(len==0)
	{
		//矩阵为空！无法提取其中的某一行！
		PMATSET(1,1);
		BACK;
	}
	if(lin<=0||lin>row)
	{
		//越界获取行向量！
		PMATSET(1,1);
		BACK;
	}
	for(unsigned j=1;j<=col;j++)
	{
		if(!(pmat->Append((*this)(lin,j))))
		{
			//无法分配堆内存以继续提取其中的某一行！
			BACK;
		}
	}
	BACK;
}

CMATRIX& CMATRIX::GetQue(unsigned que)
{
	TMPPMAT(row,1);
	if(len==0)
	{
		//矩阵为空！无法提取其中的某一列！
		PMATSET(1,1);
		BACK;
	}
	if(que<=0||que>col)
	{
		//越界获取列向量！
		PMATSET(1,1);
		BACK;
	}
	for(unsigned i=1;i<=row;i++)
	{
		if(!(pmat->Append((*this)(i,que))))
		{
			//无法分配堆内存以继续提取其中的某一列！
			BACK;
		}
	}
	BACK;
}

CMATRIX& CMATRIX::KidMat(unsigned l1, unsigned l2, unsigned q1, unsigned q2)
{
	TMPPMAT(1,1);
	if(len==0)
	{
		//矩阵为空！无法获取子矩阵！
		BACK;
	}
	if(l1<=0||l2<=0||q1<=0||q2<=0||q1>col||q2>col||l1>row||l2>row)
	{
		//越界获取子矩阵！
		BACK;
	}
	unsigned L1=(l1<=l2)?l1:l2; unsigned L2=(l1>l2)?l1:l2;
	unsigned Q1=(q1<=q2)?q1:q2; unsigned Q2=(q1>q2)?q1:q2;
	PMATSET(L2-L1+1,Q2-Q1+1); 
	for(unsigned i=L1;i<=L2;i++)
	{
		for(unsigned j=Q1;j<=Q2;j++)
		{
			if(!(pmat->Append((*this)(i,j))))
			{
				//无法分配堆内存以继续获取子矩阵！
				BACK;
			}
		}
	}
	BACK;
}

CMATRIX& inverse(CMATRIX &_m_)
{
	TMPPMAT(_m_.Row(),_m_.Col());
	if(_m_.Len()==0)
	{
		//矩阵为空！无法求其对应的逆矩阵！
		PMATSET(1,1);cout<<"\ndet=0"<<endl;
		BACK;
	}
	if(det(_m_)==0)
	{
		//奇异矩阵无法求解逆矩阵！
		PMATSET(1,1);
		BACK;
	}
//	*pmat=(1/det(_m_))*adjoint(_m_);
	pmat=&((1/det(_m_))*adjoint(_m_));
	BACK;
}

CMATRIX& CMATRIX::Identity()
{
	TMPPMAT(row,col);
	for(unsigned i=1;i<=pmat->Row();i++)
	{
		for(unsigned j=1;j<=pmat->Col();j++)
		{
			if(i==j)
			{
				if(!(pmat->Append(1)))
				{
					//无法分配堆内存以继续生成标准形的矩阵！
					BACK;
				}

			}
			else
			{
				if(!(pmat->Append(0)))
				{
					//无法分配堆内存以继续生成标准形的矩阵！
					BACK;
				}
			}
		}
	}
	BACK;
}

CMATRIX& CMATRIX::IDENTITY(unsigned dim)
{
	TMPPMAT(dim,dim);
	if(dim==0)
	{
		//指定的方阵维数为零！
		PMATSET(1,1);
		BACK;
	}
	for(unsigned i=1;i<=pmat->Row();i++)
	{
		for(unsigned j=1;j<=pmat->Col();j++)
		{	
			if(i==j)
			{
				if(!(pmat->Append(1)))
				{
					//无法分配堆内存以继续生成单位矩阵！
					BACK;
				}
			}
			else
			{
				if(!(pmat->Append(0)))
				{
					//无法分配堆内存以继续生成单位矩阵！
					BACK;
				}
			}
		}
	}
	BACK;
}

CMATRIX& CMATRIX::Basis(unsigned lin,unsigned que)
{
	TMPPMAT(row,col);
	if(len==0)
	{
		//矩阵为空！无法求其基矩阵！