matrix.cpp

模式识别HK分类算法 vc下dos程序

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

#include "stdafx.h"
//#include "Mat.h"
#include "Matrix.h"
#include "math.h"
#include "iostream.h"


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

IMPLEMENT_SERIAL(CMatrix,CObject,0)

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

CMatrix::CMatrix(unsigned int row,unsigned int col,double x)
{
	RowSize=Row=row;
	ColSize=Col=col;
	Val=new double * [RowSize];
	for (unsigned int k=0;k<RowSize;k++)
	{
		Val[k]=new double [ColSize];
	}
	int i,j;

	for(i=0;i<Row;i++)
		for(j=0;j<Col;j++)
			Val[i][j]=x;
		
}

CMatrix::~CMatrix()
{
	for (unsigned int k=0;k<RowSize;k++)
	{
		delete []Val[k];
	}
	delete []Val;


}

CMatrix::CMatrix(const CMatrix &m)
{
	RowSize=Row=m.Row;
	ColSize=Col=m.Col;
	Val=new double * [RowSize];

	for (unsigned int k=0;k<RowSize;k++)
	{
		Val[k]=new double [ColSize];
        //将内存按块拷贝
		memcpy(Val[k],m.Val[k],ColSize * sizeof(double));
	}

}

void CMatrix::realloc(unsigned int row,unsigned int col)
{
	if (row==RowSize&&col==ColSize)
	{
		Row=RowSize;
		Col=ColSize;
		return;
	}

	double ** Val1=new double * [row];

	for (unsigned int i=0;i<row;i++)
	{
		Val1[i]=new double  [col];

	}

	unsigned int ColSize=min(Col,col)*sizeof(double);
	unsigned int minRow=min(Row,row);

	for (i=0;i<minRow;i++)
	{
		memcpy(Val1[i],Val[i],ColSize);
	}

	for (i=0;i<RowSize;i++)
	{
		delete [] Val[i];
	}

	delete []Val;

	RowSize=Row=row;
	ColSize=Col=col;
	Val=Val1;

	return;
}

int  CMatrix::pivot(unsigned int row)
{
	int k=int(row);
	double t,amax=-1;

	for (unsigned int i=row;i<Row;i++)
	{
		if ((t=fabs(Val[i][row]))>amax&&t!=0.0)
		{
			amax=t;
			k=i;
		}
	}

	if (Val[k][row]==double(0))
	{
		return -1;
	}

	if (k!=int (row))
	{
		double * rowptr=Val[k];
		Val[k]=Val[row];
		Val[row]=rowptr;
		return k;
	}

	return 0;
}

double & CMatrix::operator ()(unsigned int row, unsigned int col)
{
	if (row>=Row+1||col>=Col+1||row==0||col==0)
	{
		AfxMessageBox("CMatrix::operator ():Index out of range!");
	}
	return Val[row-1][col-1];
	
}

CMatrix CMatrix::operator + ()
{
	return * this;
	
}

CMatrix CMatrix::operator - ()
{
	CMatrix T(Row,Col);

	for (unsigned int i=0;i<Row;i++)
	{
		for (unsigned int j=0;j<Col;j++)
		{
			T.Val[i][j]=-Val[i][j];
		}
	}
	return T;
}

CMatrix & CMatrix::operator = (const CMatrix &m)
{
	if(Row!=m.Row||Col!=m.Col)
	{
		realloc(m.Row,m.Col);
	}

	unsigned int colbyte=m.Col*sizeof(double);

	for (unsigned int k=0;k<m.Row;k++)
	{
		memcpy(Val[k],m.Val[k],colbyte);
	}

	return * this;
}

CMatrix & CMatrix::operator += (const CMatrix &m)
{
	if(Row!=m.Row||Col!=m.Col)
	{
		AfxMessageBox("CMatrix::operator + =:Inconsistent CMatrix Size in addition!");
	}

	for (unsigned int i=0;i<m.Row;i++)
	{
		for (unsigned int j=0;j<m.Col;j++)
		{
			Val[i][j]+=m.Val[i][j];
		}
	}

	return * this;
}

CMatrix & CMatrix::operator -= (const CMatrix &m)
{
	if(Row!=m.Row||Col!=m.Col)
	{
		AfxMessageBox("CMatrix::operator - =:Inconsistent CMatrix Size in subtration!");
	}

	for (unsigned int i=0;i<m.Row;i++)
	{
		for (unsigned int j=0;j<m.Col;j++)
		{
			Val[i][j]-=m.Val[i][j];
		}
	}

	return * this;
}

CMatrix & CMatrix::operator *= (const double& c)
{
	for (unsigned int i=0;i<Row;i++)
	{
		for (unsigned int j=0;j<Col;j++)
		{
			Val[i][j]*=c;
		}
	}

	return * this;
}

CMatrix & CMatrix::operator *= (const CMatrix &m)
{
	if(Col!=m.Row)
	{
		AfxMessageBox("CMatrix::operator * =:Inconsistent CMatrix Size in multiplication!");
	}
	
	*this=* this * m;

	return * this;
}

CMatrix & CMatrix::operator /= (const double &c)
{
	for (unsigned int i=0;i<Row;i++)
	{
		for (unsigned int j=0;j<Col;j++)
		{
			Val[i][j]/=c;
		}
	}

	return * this;
}

CMatrix & CMatrix::operator ^= (const unsigned int & pow)
{
	for (unsigned int i=2;i<=pow;i++)
	{
		* this =* this * * this;
	}

	return * this;
}

CMatrix CMatrix::operator + (const CMatrix &m)
{
	if(Row!=m.Row||Col!=m.Col)
	{
		AfxMessageBox("CMatrix::operator +:Inconsistent CMatrix Size in addition!");
	}

	CMatrix T(Row,Col);

	for (unsigned int i=0;i<Row;i++)
	{
		for (unsigned int j=0;j<Col;j++)
		{
			T.Val[i][j]=Val[i][j]+m.Val[i][j];
		}
	}
	return T;
}

CMatrix CMatrix::operator - (const CMatrix &m)
{
	if(Row!=m.Row||Col!=m.Col)
	{
		AfxMessageBox("CMatrix::operator -:Inconsistent CMatrix Size in subtraction!");
	}

	CMatrix T(Row,Col);

	for (unsigned int i=0;i<Row;i++)
	{
		for (unsigned int j=0;j<Col;j++)
		{
			T.Val[i][j]=Val[i][j]-m.Val[i][j];
		}
	}
	return T;
}

CMatrix CMatrix::operator * (const double &no)
{
	CMatrix T(Row,Col);

	for (unsigned int i=0;i<Row;i++)
	{
		for (unsigned int j=0;j<Col;j++)
		{
			T.Val[i][j]=no*Val[i][j];
		}
	}
	return T;
}

CMatrix CMatrix::operator * (const CMatrix &m)
{
	if(Col!=m.Row)
	{
		AfxMessageBox("CMatrix::operator *:Inconsistent CMatrix Size in multiplicationg!");
	}

	CMatrix T(Row,m.Col);

	for (unsigned int i=0;i<Row;i++)
	{
		
		for (unsigned int j=0;j<m.Col;j++)
		{
			T.Val[i][j]=double(0);
			for (unsigned int k=0;k<Col;k++)
			{
				T.Val[i][j]+=Val[i][k]*m.Val[k][j];
			}
		}
	}
	return T;
}

CMatrix CMatrix::operator ^ (const unsigned int & pow)
{
	CMatrix T(*this);
	for (unsigned int i=2;i<=pow;i++)
	{
		T=T* T;
	}

	return T;
}

CMatrix CMatrix::operator ~ ()
{
	CMatrix T(Col,Row);

	for (unsigned int i=0;i<Row;i++)
	{
		for (unsigned int j=0;j<Col;j++)
		{
			T.Val[j][i]=Val[i][j];
		}
	}
	return T;
}

CMatrix CMatrix::operator !()
{
	unsigned int i,j,k;
	double a1,a2,* rowptr;

	if(Row!=Col)
	{
		AfxMessageBox("CMatrix::operator !:Inversion of a non-square CMatrix");
	}

	CMatrix T(Row,Col);
	T.Unit();
	CMatrix M(*this);

	for(k=0;k<Row;k++)
	{
		int indx=M.pivot(k);

		if (indx==-1)
		{
			AfxMessageBox("CMatrix::operator !:Inversion of a singular CMatrix");
		}
		if (indx!=0)
		{
			rowptr=T.Val[k];
			T.Val[k]=T.Val[indx];
			T.Val[indx]=rowptr;
		}
		a1=M.Val[k][k];

		for (j=0;j<Row;j++)
		{
			M.Val[k][j]/=a1;
			T.Val[k][j]/=a1;
		}
		
		for (i=0;i<Row;i++)
		{
			if (i!=k)
			{
				a2=M.Val[i][k];
				for (j=0;j<Col;j++)
				{
					M.Val[i][j]-=a2*M.Val[k][j];
					T.Val[i][j]-=a2*T.Val[k][j];
				}
			}
		}
	}
	return T;
}

void CMatrix::Null(const unsigned int & row,const unsigned int & col)
{
	if(row!=Row||col!=Col)
	{
		realloc(row,col);
	}

	for (unsigned int i=0;i<Row;i++)
	{
		for (unsigned int j=0;j<Col;j++)
		{
			Val[i][j]=double (0);
		}
	}
	return;
	
}

void CMatrix::Null()
{
	for (unsigned int i=0;i<Row;i++)
	{
		for (unsigned int j=0;j<Col;j++)
		{
			Val[i][j]=double (0);
		}
	}
	return;
}

void CMatrix ::Unit(const unsigned int & row)
{
	if(row!=Row||row!=Col)
	{
		realloc(row,row);
	}

	for (unsigned int i=0;i<Row;i++)
	{
		for (unsigned int j=0;j<Col;j++)
		{
			Val[i][j]=i==j?double (1):double(0);
		}
	}
	return;
	
}

void CMatrix ::Unit()
{
	unsigned int row=min(Row,Col);
	Row=Col=row;

	for (unsigned int i=0;i<Row;i++)
	{
		for (unsigned int j=0;j<Col;j++)
		{
			Val[i][j]=i==j?double (1):double(0);
		}
	}
	return;
	
}

double CMatrix::Det()
{
	unsigned int i,j,k;
	double piv,detVal=double(1);

	if(Row!=Col)
	{
		AfxMessageBox("CMatrix::Det():Determinant a non-square CMatrix!");
	}

	CMatrix T(* this);

	for (k=0;k<Row;k++)
	{
		int indx=T.pivot(k);
		if (indx==-1)
		{
			return 0;
		}
		if (indx!=0)
		{
			detVal=-detVal;
		}

		detVal=detVal*T.Val[k][k];

		for (i=k+1;i<Row;i++)
		{
			piv=T.Val[i][k]/T.Val[k][k];
			
			for (j=k+1;j<Row;j++)
			{
				T.Val[i][j]-=piv*T.Val[k][j];
			}
		}
	}
	return detVal;
}

void CMatrix::Serialize(CArchive &ar)
{
	CObject::Serialize(ar);

	if(ar.IsStoring())
	{
		ar<<Row<<Col;
		for (unsigned i=0;i<Row;i++)
		{
			for (unsigned j=0;j<Col;j++)
			{
				ar<<Val[i][j];
			}
		}
	}
	else 
	{
		unsigned int row,col;
		ar>>row>>col;
		SetSize(row,col);

		for (unsigned i=0;i<Row;i++)
		{
			for (unsigned j=0;j<Col;j++)
			{
				ar>>Val[i][j];
			}
		}
	}

}

void CMatrix::OnDraw(CDC *pDC, CPoint Pos)
{
	char ch[50];

	for (unsigned int i=0;i<Row;i++)
	{
		for (unsigned int j=0;j<Col;j++)
		{
			if(fabs(Val[i][j])>1e3||fabs(Val[i][j])<1e-3)
			{
				sprintf(ch,"%0.3e",Val[i][j]);
			}
			else
			{
				sprintf(ch,"%0.3f",Val[i][j]);
			}

			pDC->TextOut(Pos.x+90*j,Pos.y+20*i,ch);
		}
	}
}

double CMatrix::Max()
{
	double t=-1e308;

	for (unsigned int i=0;i<Row;i++)
	{
		for (unsigned int j=0;j<Col;j++)
		{
			t=max(t,Val[i][j]);
		}
	}
	return t;
}

double CMatrix::Min()
{
	double t=1e308;

	for (unsigned int i=0;i<Row;i++)
	{
		for (unsigned int j=0;j<Col;j++)
		{
			t=min(t,Val[i][j]);
		}
	}
	return t;
}

void CMatrix::SetSize(unsigned int row,unsigned int col)
{
	unsigned int i,j;
	unsigned int oldRow=Row;
	unsigned int oldCol=Col;

	if (row!=RowSize||col!=ColSize)
	{
		realloc(row,col);
	}

	if (row>oldRow)
	{
		for(i=oldRow;i<row;i++)
		{
			for (j=0;j<col;j++)
			{
				Val[i][j]=double(0);
			}
		}
	}

	if (col>oldCol)
	{
		for(i=0;i<row;i++)
		{
			for (j=oldCol;j<col;j++)
			{
				Val[i][j]=double(0);
			}
		}
	}

	return;
}

bool CMatrix::operator == (const CMatrix &m)
{
	bool retVal=false;

	if (Row!=m.Row||Col!=m.Col)
	    return retVal;

    for (unsigned int i=0;i<Row;i++)
			for (unsigned int j=0;j<Col;j++)
				if(Val[i][j]!=m.Val[i][j])
			      return retVal;

	return true;
}

bool CMatrix::operator != (const CMatrix &m)
{

	return(*this==m)?false:true;
	
}

bool CMatrix::operator >= (const CMatrix &m)
{
	bool retVal=false;

	if (Row!=m.Row||Col!=m.Col)
	   {
		AfxMessageBox("CMatrix::>=():can not be compare");
//		return retVal;
		}

    for (unsigned int i=0;i<Row;i++)
			for (unsigned int j=0;j<Col;j++)
				if(Val[i][j]<m.Val[i][j])
					return retVal;

	return true;
}

bool CMatrix::operator > (const CMatrix &m)
{
	bool retVal=false;

	if (Row!=m.Row||Col!=m.Col)
	   {
		AfxMessageBox("CMatrix::>():can not be compare");
//		return retVal;
		}

    for (unsigned int i=0;i<Row;i++)
			for (unsigned int j=0;j<Col;j++)
				if(Val[i][j]<=m.Val[i][j])
					return retVal;

	return true;
}

bool CMatrix::operator <= (const CMatrix &m)
{
	bool retVal=false;

	if (Row!=m.Row||Col!=m.Col)
	   {
		AfxMessageBox("CMatrix::<=():can not be compare");
//		return retVal;
		}

    for (unsigned int i=0;i<Row;i++)
			for (unsigned int j=0;j<Col;j++)
				if(Val[i][j]>m.Val[i][j])
					return retVal;

	return true;
}

bool CMatrix::operator < (const CMatrix &m)
{
	bool retVal=false;

	if (Row!=m.Row||Col!=m.Col)
	   {
		AfxMessageBox("CMatrix::<():can not be compare");
//		return retVal;
		}

    for (unsigned int i=0;i<Row;i++)
			for (unsigned int j=0;j<Col;j++)
				if(Val[i][j]>=m.Val[i][j])
					return retVal;

	return true;
}

CMatrix Abs(const CMatrix &m)
{

	CMatrix T(m.Row,m.Col);
		for (unsigned int i=0;i<m.Row;i++)
	{
		for (unsigned int j=0;j<m.Col;j++)
		{
			if(m.Val[i][j]<0)
				T.Val[i][j]=-m.Val[i][j];
			else
				T.Val[i][j]=m.Val[i][j];
		}
	}
	return T;

}

//friend CMatrix operator Abs(CMatrix &m);

void CMatrix::Display()
{
	int i,j;
	for( j=0;j<Row;j++)
		{
			for( i=0;i<Col;i++)
			cout<<Val[j][i]<<"    ";
			cout<<endl;
		}
	 cout<<endl;



}