matrix.cpp

一个简单的矩阵类

//#ifndef _MATRIX_H_
//#define _MATRIX_H_
#include "StdAfx.h"
#include "math.h"

class Matrix{
public:
	double *Data;
   
private:
	int m_nRow;
	int m_nCol;
public:
	Matrix();                       //构造函数
	Matrix(int nRow,int nCol);
	Matrix(int nDim,double nDiag);   //构造单位阵
	~Matrix();
	Matrix(const Matrix &m1);              //拷贝构造函数
    Matrix operator + (const Matrix &m1);
	Matrix operator - (const Matrix &m1);
	Matrix operator * (const Matrix &m1);
	Matrix operator * (double f);
	Matrix operator / ( const Matrix &m1);   
	Matrix& operator = (const Matrix &m1);
	Matrix operator ~();
	void Set(int nRow,int nCol);
	friend double Mod(const Matrix &m1);      //求模
	friend Matrix Chol(const Matrix &m1);    //切罗斯基变换
};

//#endif


//#include "Matrix.h"

Matrix::Matrix()
{
	Data = NULL;
}


Matrix::Matrix(int nRow,int nCol)
{
    ASSERT(nRow>0 && nCol>0);
	m_nRow = nRow;
	m_nCol = nCol;
	Data = new double[m_nRow*m_nCol];
	for(int i=0;i<m_nRow*m_nCol;i++)
		Data[i]=0;
}

Matrix::Matrix(int nDim,double nDiag)
{
	ASSERT(nDim>0);
	m_nRow = m_nCol = nDim;
	Data = new double[nDim*nDim];
	for(int i=0;i<nDim*nDim;i++)
		Data[i] = 0;
	for(i=0;i<nDim;i++)
		Data[i*nDim+i] = nDiag;
}

Matrix::Matrix(const Matrix &m1)
{
	m_nRow = m1.m_nRow;
    m_nCol = m1.m_nCol;
    Data = new double[m_nRow*m_nCol];
	for(int i=0;i<m_nRow*m_nCol;i++)
		Data[i] = m1.Data[i];
}

Matrix::~Matrix()
{
	if(Data!=NULL)
		delete[] Data;
}

Matrix Matrix :: operator +(const Matrix &m1)
{
	Matrix m2(m_nRow,m_nCol);
	for(int i=0;i<m_nRow*m_nCol;i++)
		m2.Data[i] = Data[i]+m1.Data[i];
    return m2;
}

Matrix Matrix :: operator -(const Matrix &m1)
{
	Matrix m2(m_nRow,m_nCol);
	for(int i=0;i<m_nRow*m_nCol;i++)
		m2.Data[i] = Data[i]-m1.Data[i];
    return m2;
}

Matrix Matrix :: operator *(const Matrix &m1)
{

	ASSERT(m_nCol == m1.m_nRow);

	Matrix m2(m_nRow,m1.m_nCol);
	int i,j,k;
	double temp = 0;
    for( i=0;i<m_nRow;i++)
	{
    	 for(j=0;j<m1.m_nCol;j++)
		 {
    		 temp = 0;
	    	 for(k=0;k<m_nCol;k++)
    			 temp += Data[i*m_nCol+k]*m1.Data [k*m1.m_nCol+j];
			 m2.Data[i*m1.m_nCol+j] = temp;
		 }
	}
	return m2;
}

Matrix Matrix :: operator *(double f)
{
	Matrix m2(m_nRow,m_nCol);
    for(int i=0;i<m_nRow*m_nCol;i++)
	{
    	m2.Data [i] = Data[i]*f;
	}
	return m2;
}
////////////////////////////////////////////////////////////////////////////
//普通求逆，可以作为除法用(例如: M1 = M2*M3'可表为m1 = m2/m3),但不能求伪逆
/////////////////////////////////////////////////////////////////////////////

Matrix Matrix :: operator / (const Matrix &m1)    
{
    ASSERT(m1.m_nRow == m1.m_nCol);   // 只能方阵求逆
	ASSERT(m1.m_nRow == m_nCol);      // 要满足矩阵乘法要求
	int nDim;
	int i,j,k,s,dVec,dCol;
	double f;
	nDim = m1.m_nRow ;
	if(nDim == 1)
	{
		ASSERT(m1.Data[0]>1e-12 || m1.Data[0]<-1e-12);
		Matrix m5(m_nRow,m_nCol);
		for(i=0;i<m_nRow;i++)
			for(j=0;j<m_nCol;j++)
				m5.Data[i*m_nCol+j] = Data[i*m_nCol+j]/m1.Data[0];
		return m5;
	}

    double nMod;
    nMod = Mod(m1);	
	ASSERT(nMod>1e-12 || nMod<-1e-12);

	Matrix m2(nDim-1,nDim-1),m3(nDim,nDim),m4(m_nRow,m_nCol);
 
    for(i=0;i<nDim;i++)
      for(j=0;j<nDim;j++)
	  {    
		  for(k=0;k<nDim-1;k++)       //构造伴随矩阵       
			  for(s=0;s<nDim-1;s++)   
			  {
				  if(k<i)
					  dVec=0;
				  else 
					  dVec=1;
				  if(s<j)
					  dCol=0;
				  else
					  dCol=1;
				  m2.Data [k*(nDim-1)+s]=m1.Data [(k+dVec)*nDim+s+dCol]; 
			  }
		  if((i+1+j+1)%2)           //判断余子式符号
		  {
			  f=Mod(m2);
			  m3.Data [j*nDim+i]=-f/nMod;
		  }
		  else
		  {
			  f=Mod(m2);
			  m3.Data [j*nDim+i]=f/nMod;
		  }
	}
 
	double temp = 0;
    for( i=0;i<m_nRow;i++)
	{
    	 for(j=0;j<m_nCol;j++)
		 {
    		 temp = 0;
	    	 for(k=0;k<m_nCol;k++)
    			 temp += Data[i*m_nCol+k]*m3.Data [k*m_nCol+j];
			 m4.Data[i*m_nCol+j] = temp;
		 }
	}
	return m4;
}


double Mod(const Matrix &m1)      //求模
{
	int i,j,k;
	ASSERT(m1.m_nRow == m1.m_nCol);

	int nDim = m1.m_nRow;
	double temp = 0;
	if(nDim == 1)
		return m1.Data [0];
	
	Matrix m2(nDim-1,nDim-1);
	for(i=0;i<nDim;i++)
	{  
       for(j=0;j<nDim-1;j++)
	      if(j<i)
			  for(k=0;k<nDim-1;k++)
		    	   m2.Data [k*nDim-k+j]=m1.Data [k*nDim+nDim+j];
		  else
              for(k=0;k<nDim-1;k++)
		    	   m2.Data [k*nDim-k+j]=m1.Data [k*nDim+nDim+j+1];
       if(i%2)
		   temp -= m1.Data [i]*Mod(m2);
	   else
		   temp += m1.Data [i]*Mod(m2);
   
	}
    return temp;
}

Matrix Matrix :: operator ~()
{
	Matrix m2(m_nCol,m_nRow);
	for(int i=0;i<m_nRow;i++)
		for(int j=0;j<m_nCol;j++)
			m2.Data[j*m_nRow+i]=Data[i*m_nCol+j];
	return m2;
}


Matrix& Matrix :: operator =(const Matrix &m1)
{
	if(this==&m1)
		return *this;
	if(Data!=NULL)
	{
		delete[] Data;
	}
	m_nRow = m1.m_nRow;
    m_nCol = m1.m_nCol;
    Data=new double[m_nRow*m_nCol];
	for(int i=0;i<m_nRow*m_nCol;i++)
		Data[i] = m1.Data[i];
	return *this;
}

void Matrix :: Set(int nRow,int nCol)
{
	ASSERT(nRow>0&&nCol>0);
	if(Data != NULL)
		delete[] Data;
	m_nRow = nRow;
	m_nCol = nCol;
	Data = new double[m_nRow*m_nCol];
	for(int i=0;i<m_nRow*m_nCol;i++)
		Data[i] = 0;
}

Matrix Chol(const Matrix &m1)
{
	double temp=0;
	ASSERT(m1.m_nRow == m1.m_nCol);
	Matrix m2(m1.m_nRow,m1.m_nCol);
	
	for(int j=0;j<m1.m_nCol;j++)
	{
		temp = m1.Data [j*m1.m_nCol+j];
		for(int k=0;k<j;k++)
			temp -= m2.Data [j*m1.m_nCol+k]*m2.Data [j*m1.m_nCol+k];
		m2.Data [j*m1.m_nCol+j] = sqrt(temp);
		for(int i=j+1;i<m1.m_nRow;i++)
		{
			temp = m1.Data [i*m1.m_nCol+j];
			for(k=0;k<j;k++)
			{
				temp -= m2.Data [i*m1.m_nCol+k]*m2.Data [j*m1.m_nCol+k];
			}
			ASSERT(m2.Data [j*m1.m_nCol+j]>1e-11 || m2.Data [j*m1.m_nCol+j]<-1e-11);
			m2.Data [i*m1.m_nCol+j]=temp/m2.Data [j*m1.m_nCol+j];
		}
	}
	return m2;
}