matrix.cpp

一个最小二乘法的拟合数据的小程序

// Matrix.cpp: implementation of the CMatrix class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "Matrix.h"
#include "math.h"

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

CMatrix::CMatrix()
{
	m_nNumColumns=1;
	m_nNumRows=1;
	m_pData=NULL;
	Init(m_nNumRows,m_nNumColumns);
}

CMatrix::~CMatrix()//释放矩阵所占有内存
{
	if (m_pData) 
	{
		delete[] m_pData;
		m_pData=NULL;
	}
}

CMatrix::CMatrix(int nRows,int nCols)//指定行列构造函数
{
	m_nNumRows=nRows;
	m_nNumColumns=nCols;
	m_pData=NULL;
	Init(m_nNumRows,m_nNumColumns);
}
CMatrix::CMatrix(int nRows,int nCols,double value[])//指定数据构造函数
{
	m_nNumRows=nRows;
	m_nNumColumns=nCols;
	m_pData=NULL;
	Init(m_nNumRows,m_nNumColumns);
	
	SetData(value);
}
CMatrix::CMatrix(int nSize)//方阵构造函数
{
	m_nNumColumns=nSize;
	m_nNumRows=nSize;
	m_pData=NULL;
	Init(m_nNumRows,m_nNumColumns);
}
CMatrix::CMatrix(int nSize,double value[])//指定数据方阵构造函数
{
	m_nNumColumns=nSize;
	m_nNumRows=nSize;
	m_pData=NULL;
	Init(nSize,nSize);
	
	SetData(value);
}
CMatrix::CMatrix(const CMatrix& other)//拷贝构造函数
{
	m_nNumColumns=other.GetNumColums();
	m_nNumRows=other.GetNumRows();
	m_pData=NULL;
	Init(m_nNumRows,m_nNumColumns);
	
	
	//copy the pointer
	memcpy(m_pData,other.m_pData,
		sizeof(double)*m_nNumColumns*m_nNumRows);
}

BOOL CMatrix::Init(int nRows,int nCols)//初始化矩阵
{
	if(m_pData)
	{
		delete[] m_pData;
		m_pData=NULL;
	}
	m_nNumColumns=nCols;
	m_nNumRows=nRows;
	int nSize=nCols*nRows;
	if (nSize<0) {return FALSE;}
	
	//分配内存
	m_pData=new double[nSize];
	if(m_pData==NULL)return FALSE;
	if(IsBadReadPtr(m_pData,sizeof(double)*nSize))return FALSE;
	
	//将各元素值置0
	memset(m_pData,0,sizeof(double)*nSize);
	return TRUE;
}
BOOL CMatrix::MakeUnitMatrix(int nSize)//将方阵初始化为单位矩阵
{
	if (!Init(nSize,nSize) )return FALSE;
	for(int i=0;i<nSize;i++)
		for(int j=0;j<nSize;j++)
			if(i==j)SetElement(i,j,1);
			
			return TRUE;
}


//元素与值操作

BOOL CMatrix::SetElement(int nRow,int nCol,double value)//指定元素的值
{
	if(nCol<0||nCol>=m_nNumColumns||nRow<0||nRow>=m_nNumRows)
		return FALSE;//array bounds error
	if(m_pData==NULL)return FALSE;//bad pointer error
	
	m_pData[nCol+nRow*m_nNumColumns]=value;
	
	return TRUE;
}
double CMatrix::GetElement(int nRow,int nCol)const//获取元素的值
{
	 ASSERT(nCol>=0&&nCol<m_nNumColumns&&nRow>=0&&nRow<m_nNumRows);////array bounds error
	ASSERT(m_pData);
	return m_pData[nCol+nRow*m_nNumColumns];
}
void CMatrix::SetData(double value[])//设定矩阵的值
{
    //empty the memory
	memset(m_pData,0,sizeof(double)*m_nNumColumns*m_nNumRows);
	//copy data
	memcpy(m_pData,value,sizeof(double)*m_nNumColumns*m_nNumRows);
	
}
double * CMatrix::GetData()const//获取矩阵的值
{
	return m_pData;
}
int CMatrix::GetNumColums()const//获取矩阵的列数
{
	return m_nNumColumns;
}
int CMatrix::GetNumRows()const//获取矩阵的行数
{
	return m_nNumRows;
}


///////////////////////
//数学操作
CMatrix& CMatrix::operator=(const CMatrix & other)
{
	if (&other!=this) {
		m_nNumColumns=other.GetNumColums();
		m_nNumRows=other.GetNumRows();
		
		Init(m_nNumRows,m_nNumColumns);
		
		//copy the pointer
		memcpy(m_pData,other.m_pData,sizeof(double)*m_nNumColumns*m_nNumRows);
		
	}
	
	//finally return a referrence to ourselves
	return *this;
}
BOOL CMatrix::operator==(const CMatrix& other)const
{
//首先检查行列数是否相等
	if(m_nNumColumns!=other.GetNumColums()||m_nNumRows!=other.GetNumRows())
		return FALSE;

	for(int i=0;i<m_nNumRows;++i)
	{
		for(int j=0;j<m_nNumColumns;++j)
		{
			if(GetElement(i,j)!=other.GetElement(i,j))
				return FALSE;
		}
	}
	return TRUE;
}
BOOL CMatrix::operator!=(const CMatrix& other)const
{
return !(*this==other);
}
CMatrix CMatrix::operator+(const CMatrix& other)const
{

	//检查行列数是否相等
	ASSERT( m_nNumColumns==other.GetNumColums()&&m_nNumRows==other.GetNumRows());

	//构造结果矩阵
	CMatrix result(*this);//拷贝构造
	//矩阵加法
	for(int i=0;i<m_nNumRows;++i)
	{
		for(int j=0;j<m_nNumColumns;++j)
		{
			result.SetElement(i,j,result.GetElement(i,j)+other.GetElement(i,j));
		}
		
	}
	return result;
}
CMatrix CMatrix::operator-(const CMatrix& other)const
{
//首先检查行列数是否相等
	ASSERT(m_nNumColumns==other.GetNumColums()&&m_nNumRows==other.GetNumRows());

	//构造目标矩阵
	CMatrix result(*this);
	//进行减法操作
	for(int i=0;i<m_nNumRows;++i)
	{
		for(int j=0;j<m_nNumColumns;++j)
		
			result.SetElement(i,j,result.GetElement(i,j)-other.GetElement(i,j));
			}
	return result;
}
CMatrix CMatrix::operator*(const CMatrix& other)const
{
//首先检查行列数是否符合要求
	ASSERT(m_nNumColumns==other.GetNumRows());

	//construct the object we are going to return
	CMatrix result(m_nNumRows,other.GetNumColums());

	//矩阵乘法，即
	double value;
	for(int i=0;i<result.GetNumRows();++i)
	{
		for(int j=0;j<other.GetNumColums();++j)
		{
			value=0.0;
			for(int k=0;k<m_nNumColumns;++k)
			{
				value+=GetElement(i,k)*other.GetElement(k,j);
			}
			result.SetElement(i,j,value);

		}
	}
	return result;
}
CMatrix CMatrix::operator*(double value)const
{
	//构造目标矩阵
	CMatrix result(*this);//copy ourselves
	//进行数乘
	for(int i=0;i<m_nNumRows;++i)
	{
		for(int j=0;j<m_nNumColumns;++j)
			result.SetElement(i,j,result.GetElement(i,j)*value);
	}
	return result;
}


//转置
CMatrix CMatrix::Transpose()const
{
//构造目标矩阵
	CMatrix Trans(m_nNumColumns,m_nNumRows);

	//转置各元素
	for(int i=0;i<m_nNumRows;++i)
	{
		for(int j=0;j<m_nNumColumns;++j)
			Trans.SetElement(j,i,GetElement(i,j));
	}
	return Trans;
}

//实矩阵求逆——全选主元高斯-约旦法
BOOL CMatrix::InvertGaussJordan()
{
	int *pnRow,*pnCol,i,j,k,l,u,v;
	double d=0,p=0;
	
	//分配内存
	pnRow=new int[m_nNumColumns];
	pnCol=new int[m_nNumColumns];
	if(pnRow==NULL||pnCol==NULL)return FALSE;
	
	//消元
	for(k=0;k<=m_nNumColumns-1;k++)
	{
		d=0.0;
		for(i=k;i<=m_nNumColumns-1;i++)
		{
			for(j=k;j<=m_nNumColumns-1;j++)
			{
				l=i*m_nNumColumns+j;
				p=fabs(m_pData[l]);
				if(p>d)
				{
					d=p;
					pnRow[k]=i;
					pnCol[k]=j;
					
				}
			}
		}

		//失败
		if (d==0.0) {
			delete[]pnRow;
			delete[]pnCol;
			return FALSE;
		}
		if (pnRow[k]!=k) {
			for(j=0;j<=m_nNumColumns-1;j++)
			{
				u=k*m_nNumColumns+j;
				v=pnRow[k]*m_nNumColumns+j;
				p=m_pData[u];
				m_pData[u]=m_pData[v];
				m_pData[v]=p;
			}
		}
		if (pnCol[k]!=k) {
			for(i=0;i<=m_nNumColumns-1;i++)
			{
				u=i*m_nNumColumns+k;
				v=pnCol[k]+i*m_nNumColumns;
				p=m_pData[u];
				m_pData[u]=m_pData[v];
				m_pData[v]=p;
			}
		}

		l=k*m_nNumColumns+k;
		m_pData[l]=1.0/m_pData[l];
		for(j=0;j<=m_nNumColumns-1;j++)
		{
			if (j!=k) {
				u=k*m_nNumColumns+j;
				m_pData[u]=m_pData[u]*m_pData[l];
			}
		}
		for(i=0;i<=m_nNumColumns-1;i++)
		{
			if (i!=k) {
				for(j=0;j<=m_nNumColumns-1;j++)
				{
					if (j!=k) {
						u=i*m_nNumColumns+j;

						m_pData[u]=m_pData[u]-m_pData[i*m_nNumColumns+k]*
							m_pData[k*m_nNumColumns+j];					}
				}
			}
		}
		for(i=0;i<=m_nNumColumns-1;i++)
		{
			if (i!=k) {
				u=i*m_nNumColumns+k;
				m_pData[u]=-m_pData[u]*m_pData[l];
			}
		}
	}
	//调整恢复行列次序
	for(k=m_nNumColumns-1;k>=0;k--)
	{
		if (pnCol[k]!=k) {
			for(j=0;j<=m_nNumColumns-1;j++)
			{
				u=k*m_nNumColumns+j;
				v=pnCol[k]*m_nNumColumns+j;
				p=m_pData[u];
				m_pData[u]=m_pData[v];
				m_pData[v]=p;
			
			}
		}
		if (pnRow[k]!=k) {
			for(i=0;i<=m_nNumColumns-1;i++)
			{
				u=i*m_nNumColumns+k;
				v=i*m_nNumColumns+pnRow[k];
				p=m_pData[u];
				m_pData[u]=m_pData[v];
				m_pData[v]=p;
			}
		}
	}

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

	//成功返回
	return TRUE;
}

//方阵行列式的值

double CMatrix::DetGauss()
{
	int i,j,k,is,js,l,u,v;
	double f,det,q,d;
	
	//初值
	f=1.0;
	det=1.0;
	
	//消元
	for(k=0;k<=m_nNumColumns-2;k++)
	{
		q=0.0;
		for(i=k;i<=m_nNumColumns-1;i++)
		{
			for(j=k;j<=m_nNumColumns-1;j++)
			{
				l=i*m_nNumColumns+j;
				d=fabs(m_pData[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<=m_nNumColumns-1;j++)
			{
				u=k*m_nNumColumns+j;
				v=is*m_nNumColumns+j;
				d=m_pData[u];
				m_pData[u]=m_pData[v];
				m_pData[v]=d;
			}
		}
		if (js!=k) {
			f=-f;
			for(i=k;i<=m_nNumColumns-1;i++)
			{
				u=i*m_nNumColumns+js;
				v=i*m_nNumColumns+k;
				d=m_pData[u];
				m_pData[u]=m_pData[v];
				m_pData[v]=d;
			}
		}
		l=k*m_nNumColumns+k;
		det=det*m_pData[l];
		for(i=k+l;i<=m_nNumColumns-1;i++)
		{
			d=m_pData[i*m_nNumColumns+k]/m_pData[l];
			for(j=k+1;j<=m_nNumColumns-1;j++)
			{
				u=i*m_nNumColumns+j;
				m_pData[u]=m_pData[u]-d*m_pData[k*m_nNumColumns+j];
			}
			
		}
		
	}	
	
		//求值
		det=f*det*m_pData[m_nNumColumns*m_nNumColumns-1];
		return(det);
		
		

}



//求秩

int CMatrix::RankGauss()
{
	
	int i,j,k,nn,is,js,l,ll,u,v;
	double q,d;
	
	//秩小于等于行列数
	nn=m_nNumRows;
	if (m_nNumRows>=m_nNumColumns) 
		nn=m_nNumColumns;
	
	
	k=0;
	
	//消元求解
	for(l=0;l<=nn-1;l++)
	{
		q=0.0;
		for(i=1;i<=m_nNumRows-1;i++)
		{
			for(j=l;j<=m_nNumColumns-1;j++)
			{
				ll=i*m_nNumColumns+j;
				d=fabs(m_pData[ll]);
				if (d>q)
				{
					q=d;
					is=i;
					js=j;
				}
			}
		}
		if (q==0.0) return(k);
		
		k=k+l;
		if (is!=l)
		{
			for(j=l;j<=m_nNumColumns-1;j++)
			{
				u=l*m_nNumColumns+j;
				v=is*m_nNumColumns+j;
				d=m_pData[u];
				m_pData[u]=m_pData[v];
				m_pData[v]=d;
			}
		}
		if (js!=l) 
		{
			for(i=l;i<=m_nNumColumns-1;i++)
			{
				u=i*m_nNumColumns+js;
				v=i*m_nNumColumns+l;
				d=m_pData[u];
				m_pData[u]=m_pData[v];
				m_pData[v]=d;
			}
		}
		
		ll=l*m_nNumColumns+l;
		for(i=l+1;i<=m_nNumColumns-1;i++)
		{
			d=m_pData[i*m_nNumColumns+l]/m_pData[ll];
			for(j=l+1;j<=m_nNumColumns-1;j++)
			{
				u=i*m_nNumColumns+j;
				m_pData[u]=m_pData[u]-d*m_pData[l*m_nNumColumns+j];
			}
		}
	}
	return(k);
}