arraymatrix.cpp

实现常用的矩阵运算

// ArrayMatrix.cpp: implementation of the CArrayMatrix class.
//
//////////////////////////////////////////////////////////////////////

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

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
IMPLEMENT_SERIAL(CArrayMatrix, CObject, 1)

// ArrayMatrix.cpp : implementation file
//

#include "stdafx.h"
#include "ArrayMatrix.h"
#include <cmath>


// CArrayMatrix

CArrayMatrix::CArrayMatrix()
{
	m_Name=_T("");
	m_wCon=0;
	m_wRow=0;
}

CArrayMatrix::CArrayMatrix(UINT m_row,UINT m_col)
{
		m_Name=_T("");
		m_wRow=m_row;
		m_wCon=m_col;
		m_Array.SetSize(m_wRow*m_wCon);
		for(int i=0;i<int(m_wRow*m_wCon);i++)
		{
			m_Array[i]=0;
		}
	
}
CArrayMatrix::CArrayMatrix(const CArrayMatrix & m)
{
	m_Name=m.m_Name;
	m_wRow=m.m_wRow;
	m_wCon=m.m_wCon;
	m_Array.SetSize(m_wRow*m_wCon);
	for(int i=0;i<int(m_wRow*m_wCon);i++)
	{
		m_Array[i]=m.m_Array[i];
	}
}

CArrayMatrix::CArrayMatrix(CString m_name,UINT m_row,UINT m_col, double**data)
{
	m_Name=m_name;
	m_wRow=m_row;
	m_wCon=m_col;
	m_Array.SetSize(m_row*m_col);
	for (UINT i=0;i<m_row;i++)
		for(UINT j=0;j<m_col;j++)
			m_Array[i*m_col+j]=data[i][j];

}
CArrayMatrix::CArrayMatrix(CString m_name,UINT m_row,UINT m_col, float**data)
{
	m_Name=m_name;
	m_wRow=m_row;
	m_wCon=m_col;
	m_Array.SetSize(m_row*m_col);
	for (UINT i=0;i<m_row;i++)
		for(UINT j=0;j<m_col;j++)
			m_Array[i*m_col+j]=data[i][j];

}
CArrayMatrix::CArrayMatrix(UINT m_row,UINT m_col,double**data)
{
	m_Name="";
	m_wRow=m_row;
	m_wCon=m_col;
	m_Array.SetSize(m_row*m_col);
	for (UINT i=0;i<m_row;i++)
		for(UINT j=0;j<m_col;j++)
			m_Array[i*m_col+j]=data[i][j];

}
CArrayMatrix::CArrayMatrix(UINT m_row,UINT m_col,float**data)
{
	m_Name="";
	m_wRow=m_row;
	m_wCon=m_col;
	m_Array.SetSize(m_row*m_col);
	for (UINT i=0;i<m_row;i++)
		for(UINT j=0;j<m_col;j++)
			m_Array[i*m_col+j]=data[i][j];

}
CArrayMatrix::CArrayMatrix(UINT m_row,UINT m_col,double *data)
{
	m_Name="";
	m_wRow=m_row;
	m_wCon=m_col;
	m_Array.SetSize(m_row*m_col);
	for (UINT i=0;i<m_row;i++)
		for(UINT j=0;j<m_col;j++)
			m_Array[i*m_col+j]=data[i*m_col+j];

}
CArrayMatrix::CArrayMatrix(UINT m_row,UINT m_col,float *data)
{
	m_Name="";
	m_wRow=m_row;
	m_wCon=m_col;
	m_Array.SetSize(m_row*m_col);
	for (UINT i=0;i<m_row;i++)
		for(UINT j=0;j<m_col;j++)
			m_Array[i*m_col+j]=data[i*m_col+j];

}
BOOL CArrayMatrix::Initialize(UINT m_row,UINT m_col)
{
	
	m_Name="";
	m_wRow=m_row;
	m_wCon=m_col;
	m_Array.SetSize(m_wCon*m_wRow);
	for(UINT i=0;i<m_wRow*m_wCon;i++)
		m_Array[i]=0;
	return TRUE;

}
void CArrayMatrix::SetData(UINT m_row,UINT m_col,double data)
{
	int temple=this->GetCol();
	m_Array[m_row*temple+m_col]=data;
}
CArrayMatrix::~CArrayMatrix()
{
}




// CArrayMatrix member functions
BOOL CArrayMatrix::CanAddSub(const CArrayMatrix &m)const
{
	if((int)m_wCon==(int)m.m_wCon&&(int)m.m_wRow==(int)m_wRow) return TRUE;
	return FALSE;
}

BOOL CArrayMatrix::CanContrary()const
//矩阵是方阵
{
	if(int(m_wRow)!=int(m_wCon)) return FALSE;
	if(GetRange()==0.0) return FALSE;
	return TRUE;
}
BOOL CArrayMatrix::CanMutiply(const CArrayMatrix &m)const
{
	if(m_wCon==0||m.m_wCon==0) return FALSE;
	if(m_wCon==1&&m_wRow==1) return TRUE;
	if(m.m_wCon==1&&m.m_wRow==1) return TRUE;
	if(m_wCon==m.m_wRow) return TRUE;
	return FALSE;
}

CArrayMatrix CArrayMatrix::GetAccompany()const
{
	CArrayMatrix m;
	m.m_wCon=m_wCon;
	m.m_wRow=m_wRow;
	m.m_Array.SetSize(m_Array.GetSize());
	for(int i=0;i<int(m_wRow);i++)
	{
		for(int j=0;j<int(m_wCon);j++)
		{
			CArrayMatrix tp;
			tp.m_wCon=m_wCon-1;
			tp.m_wRow=m_wRow-1;
			tp.m_Array.SetSize(tp.m_wCon*tp.m_wRow);
			int num=0;
			for(int g=0;g<int(m_wRow);g++)
			{
				if(g==i) continue;
				for(int h=0;h<int(m_wCon);h++)
				{
					if(h==j) continue;
					else
					{
						tp.m_Array[num]=m_Array[g*m_wCon+h];
						num++;
					}
				}
			}
			int tip;
			double fuhao;
			tip=i+j+2;
			fuhao=1.0;
			if(tip%2==1) fuhao=-1.0;
			m.m_Array[i*m_wCon+j]=fuhao*tp.GetRange();
		}
	}
	m=m.T();
	return m;
}
CString CArrayMatrix::GetName()
{
	return m_Name;
}
UINT CArrayMatrix::GetRow()
{
	return m_wRow;
}
UINT CArrayMatrix::GetCol()
{
	return m_wCon;
}
double CArrayMatrix::GetData(UINT row,UINT col)
{
	return m_Array[row*m_wCon+col];
}

double CArrayMatrix::GetRange() const
//求行列式值的全选主元高斯消去法
{
	if(m_wRow==0||m_wCon==0) return 0.0;
	if(m_wRow==1&&m_wCon==1) return m_Array[0];
	if(m_wRow==2&&m_wCon==2) return m_Array[0]*m_Array[3]-m_Array[1]*m_Array[2];
	double f,det;
	CArray<double,double&>temple;
	temple.SetSize(m_wCon*m_wRow);
	for(int i=0;i<int(m_wRow*m_wCon);i++)
		temple[i]=m_Array[i];
	
	int i,j,k,is,js,l,u,v;
	double q,d;

	// 初值
	f=1.0; 
	det=1.0;

	// 消元
	for (k=0; k<=int(m_wRow-2); k++)
	{ 
		q=0.0;
	    for (i=k; i<=int(m_wCon)-1; i++)
		{
			for (j=k; j<=INT(m_wCon)-1; j++)
			{ 
				l=i*m_wCon+j; 
				d=fabs(temple[l]);
				if (d>q) 
				{ 
					q=d; 
	                is=i; 
	                js=j;
				}
			}
		}
		if (q == 0.0)
		{ 
			det=0.0; 
			return(det);
		}
		if (is!=k)
		{ 
			f=-f;
			for (j=k; j<=int(m_wCon)-1; j++)
			{ 
				u=k*m_wCon+j; 
				v=is*m_wCon+j;
				d=temple[u]; 
				temple[u]=temple[v]; 
				temple[v]=d;
			}
		}
		if (js!=k)
		{ 
			f=-f;
			for (i=k; i<=int(m_wCon)-1; i++)
			{
				u=i*m_wCon+js; 
	            v=i*m_wCon+k;
	            d=temple[u]; 
	            temple[u]=temple[v]; 
	            temple[v]=d;
			}
		}
		l=k*m_wCon+k;
	    det=det*temple[l];
	    for (i=k+1; i<=int(m_wCon)-1; i++)
		{ 
			d=temple[i*m_wCon+k]/temple[l];
			for (j=k+1; j<=int(m_wCon)-1; j++)
			{ 
				u=i*m_wCon+j;
				temple[u]=temple[u]-d*temple[k*m_wCon+j];
			}
		}
	}

	// 求值
	det=f*det*temple[m_wCon*m_wCon-1];
	

	return(det);
}
BOOL CArrayMatrix::SetName(CString m)
{
	m.TrimLeft();
	m.TrimRight();
	if(lstrlen(m)==0) return FALSE;
	m_Name=m;
	return TRUE;
}
CArrayMatrix CArrayMatrix::T()const
//求转置矩阵
{
	CArrayMatrix temp;
	if(m_wCon==0||m_wRow==0) return temp;
	if(m_wCon==1&&m_wRow==1) return *this;
	temp.m_Array.SetSize(m_Array.GetSize());
	temp.m_wCon=m_wRow;
	temp.m_wRow=m_wCon;
	for(int i=0;i<int(m_wRow);i++)
	{
		for(int j=0;j<int(m_wCon);j++)
		{
			temp.m_Array[j*temp.m_wCon+i]=m_Array[i*m_wCon+j];
		}
	}
	return temp;
}
CArrayMatrix CArrayMatrix::operator +(const CArrayMatrix & m)
{
	CArrayMatrix temp;
	temp.m_wCon=m_wCon;
	temp.m_wRow=m_wRow;
	temp.m_Array.SetSize(m_wCon*m_wRow);
	for(int i=0;i<int(m_wCon*m_wRow);i++)
	{
		temp.m_Array[i]=m_Array[i]+m.m_Array[i];
	}
	return temp;
}
CArrayMatrix & CArrayMatrix::operator =(const CArrayMatrix & m)
{
	m_wRow=m.m_wRow;
	m_wCon=m.m_wCon;
	m_Array.SetSize(m_wRow*m_wCon);
	for(int i=0;i<int(m_wRow*m_wCon);i++)
	{
		m_Array[i]=m.m_Array[i];
	}
	return *this;
}
CArrayMatrix & CArrayMatrix::operator =(double m)
{
	m_wRow=1;
	m_wCon=1;
	m_Array.SetSize(1);
	m_Array.SetAt(0,m);
	return *this;
}
CArrayMatrix CArrayMatrix::operator -(const CArrayMatrix & m)
{
	CArrayMatrix temp;
	temp.m_wCon=m_wCon;
	temp.m_wRow=m_wRow;
	temp.m_Array.SetSize(m_wCon*m_wRow);
	for(int i=0;i<int(m_wCon*m_wRow);i++)
	{
		temp.m_Array[i]=m_Array[i]-m.m_Array[i];
	}
	return temp;	
}
CArrayMatrix CArrayMatrix::operator *(const CArrayMatrix &m)
{
	/*if(m_wCon==1)
	{
		CArrayMatrix temp;
		temp.m_Array.SetSize(m.m_Array.GetSize());
		temp.m_wCon=m.m_wCon;
		temp.m_wRow=m.m_wRow;
		for(int i=0;i<m.m_Array.GetSize();i++)
		{
			temp.m_Array[i]=m_Array[0]*m.m_Array[i];
		}
		return temp;
	}
	if(m.m_wCon==1)
	{
		CArrayMatrix temp;
		temp.m_Array.SetSize(m_Array.GetSize());
		temp.m_wCon=m_wCon;
		temp.m_wRow=m_wRow;
		for(int i=0;i<m_Array.GetSize();i++)
		{
			temp.m_Array[i]=m_Array[i]*m.m_Array[0];
		}
		return temp;
	}*/
	CArrayMatrix temp;
	temp.m_wRow=m_wRow;
	temp.m_wCon=m.m_wCon;
	temp.m_Array.SetSize(temp.m_wCon*temp.m_wRow);
	for(int i=0;i<int(temp.m_wRow);i++)
	{
		for(int j=0;j<int(temp.m_wCon);j++)
		{
			double total=0.0;
			for(int n=0;n<int(m_wCon);n++)
			{
				total+=m_Array[i*m_wCon+n]*m.m_Array[n*m.m_wCon+j];
			}
			temp.m_Array[i*temp.m_wCon+j]=total;
		}
	}
	
	return temp;	
}

CArrayMatrix CArrayMatrix::operator / (CArrayMatrix &m)
{
	CArrayMatrix tp;
	tp=(*this)*~m;
	return tp;
}
CArrayMatrix CArrayMatrix::operator *(double m)
{
	CArrayMatrix tp=*this;
	for(int i=0;i<int(m_wRow*m_wCon);i++)
	{
		tp.m_Array[i]*=m;
	}
	return tp;
}

CArrayMatrix CArrayMatrix::operator /(double m)
{
	CArrayMatrix tp=*this;
	if(m==0) return tp;
	for(int i=0;i<int(m_wRow*m_wCon);i++)
	{
		tp.m_Array[i]*=double(1.0/m);
	}
	return tp;
}



CArrayMatrix CArrayMatrix::operator ~()
//求逆矩阵
{
	CArrayMatrix tp;
	tp.m_wCon=1;
	tp.m_wRow=1;
	tp.m_Array.SetSize(1);
	tp.m_Array[0]=double(1.0/m_Array[0]);
	if(m_wRow==1&&m_wCon==1) return tp;
	tp.m_wCon=m_wCon;
	tp.m_wRow=m_wRow;
	tp.m_Array.SetSize(m_wCon*m_wRow);
	for(int i=0;i<int(m_wRow*m_wCon);i++)
		tp.m_Array[i]=m_Array[i];
    int *pnRow, *pnCol,i,j,k,l,u,v;
	double d = 0, p = 0;

	// 分配内存
	pnRow = new int[m_wCon];
	pnCol = new int[m_wCon];

    // 消元
	for (k=0; k<=int(m_wCon)-1; k++)
	{ 
		d=0.0;
		for (i=k; i<=int(m_wCon)-1; i++)
		{
			for (j=k; j<=int(m_wCon)-1; j++)
			{ 
				l=i*m_wCon+j; p=fabs(tp.m_Array[l]);
				if (p>d) 
				{ 
					d=p; 
					pnRow[k]=i; 
					pnCol[k]=j;
				}
			}
		}

		// 失败
		if (d == 0.0)
		{
			delete[] pnRow;
			delete[] pnCol;
			CArrayMatrix a;
			return a;
		}

		if (pnRow[k] != k)
		{
			for (j=0; j<=int(m_wCon)-1; j++)
			{ 
				u=k*m_wCon+j; 
				v=pnRow[k]*m_wCon+j;
				p=tp.m_Array[u]; 
				tp.m_Array[u]=tp.m_Array[v]; 
				tp.m_Array[v]=p;
			}
		}

		if (pnCol[k] != k)
		{
			for (i=0; i<=int(m_wCon)-1; i++)
			{ 
				u=i*m_wCon+k; 
				v=i*m_wCon+pnCol[k];
				p=tp.m_Array[u]; 
				tp.m_Array[u]=tp.m_Array[v]; 
				tp.m_Array[v]=p;
			}
		}

		l=k*m_wCon+k;
		tp.m_Array[l]=1.0/tp.m_Array[l];
		for (j=0; j<=int(m_wCon)-1; j++)
		{
			if (j != k)
			{ 
				u=k*m_wCon+j; 
				tp.m_Array[u]=tp.m_Array[u]*tp.m_Array[l];
			}
		}

		for (i=0; i<=int(m_wCon)-1; i++)
		{
			if (i!=k)
			{
				for (j=0; j<=int(m_wCon)-1; j++)
				{
					if (j!=k)
					{ 
						u=i*m_wCon+j;
						tp.m_Array[u]=tp.m_Array[u]-tp.m_Array[i*m_wCon+k]*tp.m_Array[k*m_wCon+j];
					}
				}
			}
		}

		for (i=0; i<=int(m_wCon)-1; i++)
		{
			if (i!=k)
			{ 
				u=i*m_wCon+k; 
				tp.m_Array[u]=-tp.m_Array[u]*tp.m_Array[l];
			}
		}
	}

	// 调整恢复行列次序
	for (k=int(m_wCon)-1; k>=0; k--)
	{ 
		if (pnCol[k]!=k)
		{
			for (j=0; j<=int(m_wCon)-1; j++)
			{ 
				u=k*m_wCon+j; 
				v=pnCol[k]*m_wCon+j;
				p=tp.m_Array[u]; 
				tp.m_Array[u]=tp.m_Array[v]; 
				tp.m_Array[v]=p;
			}
		}

		if (pnRow[k]!=k)
		{
			for (i=0; i<=int(m_wCon)-1; i++)
			{ 
				u=i*m_wCon+k; 
				v=i*m_wCon+pnRow[k];
				p=tp.m_Array[u]; 
				tp.m_Array[u]=tp.m_Array[v]; 
				tp.m_Array[v]=p;
			}
		}
	}

	// 清理内存
	delete[] pnRow;
	delete[] pnCol;

	// 成功返回
	return tp;
}