matrix.cpp

矩阵类 可实现转置

// CMatrix.cpp : implementation file
//

#include "stdafx.h"
#include "MatrixTemp.h"
#include "Matrix.h"
#include "Math.h"

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

/////////////////////////////////////////////////////////////////////////////
// CCMatrix

IMPLEMENT_DYNCREATE(CMatrix, CDocument)



BOOL CMatrix::OnNewDocument()
{
	if (!CDocument::OnNewDocument())
		return FALSE;
	return TRUE;
}



BEGIN_MESSAGE_MAP(CMatrix, CDocument)
	//{{AFX_MSG_MAP(CCMatrix)
		// NOTE - the ClassWizard will add and remove mapping macros here.
	//}}AFX_MSG_MAP
END_MESSAGE_MAP()

/////////////////////////////////////////////////////////////////////////////
// CCMatrix diagnostics

#ifdef _DEBUG
void CMatrix::AssertValid() const
{
	CDocument::AssertValid();
}

void CMatrix::Dump(CDumpContext& dc) const
{
	CDocument::Dump(dc);
}
#endif //_DEBUG

/////////////////////////////////////////////////////////////////////////////
// CCMatrix serialization

void CMatrix::Serialize(CArchive& ar)
{
	if (ar.IsStoring())
	{
		// TODO: add storing code here
	}
	else
	{
		// TODO: add loading code here
	}
}

/////////////////////////////////////////////////////////////////////////////
// CCMatrix commands

CMatrix::CMatrix(void)
{
 p=new double[1];
 width=0;
 height=0;
}

CMatrix::CMatrix(long n)
{
 height=width=n;
 p=new double[n*n];
 for(long i=0;i<n;i++){
  for(long j=0;j<n;j++){
   if(i==j) *(p+n*i+j)=1;
   else *(p+n*i+j)=0;
  }
 }
}

CMatrix::CMatrix(double * arrAddress,long arrWidth)
{
 long arrHeight=1;
 p=new double[arrWidth*arrHeight];
 for(long i=0;i<arrHeight;i++){
  for(long j=0;j<arrWidth;j++){
   *(p+arrWidth*i+j)=*(arrAddress+arrWidth*i+j);
  }
 }
 width=arrWidth;
 height=arrHeight;
}

CMatrix::CMatrix(double * arrAddress,long arrHeight,long arrWidth)
{
 p=new double[arrWidth*arrHeight];
 for(long i=0;i<arrHeight;i++){
  for(long j=0;j<arrWidth;j++){
   *(p+arrWidth*i+j)=*(arrAddress+arrWidth*i+j);
  }
 }
 width=arrWidth;
 height=arrHeight;
}

CMatrix::CMatrix(const CMatrix & m)//copy constructor
{
 height=m.height;
 width=m.width;
 p=new double[height*width];
 for(long i=0;i<height;i++){
   for(long j=0;j<width;j++){
    *(p+width*i+j)=*(m.p+width*i+j);
   }
 }
}

long CMatrix::getWidth(){
 return width;
}

long CMatrix::getHeight(){
 return height;
}

bool CMatrix::isVector()
{
 return !(width==1 && height==1);
}
CMatrix CMatrix::subCMatrix(long offset)
{
 ASSERT(height==width && offset<=width && offset>=0);
 double * t=new double[offset*offset];
 for(long i=0;i<offset;i++){
  for(long j=0;j<offset;j++){
   *(t+offset*i+j)=*(p+width*i+j);
  }
 }
 CMatrix m(t,offset,offset);
 delete [] t;
 return m;
}

double CMatrix::Arg()//矩阵的行列式
{
ASSERT(width==height);
 double result=1;
 CMatrix m = this->GetTopTrian();
 for(long i=0;i<height;i++)
  result*=m[i][i];
  return result;
}

bool CMatrix::isPtv()
{
 ASSERT(width==height);//是方阵才可以计算
 bool result=true;
 CMatrix m;
 for(long i=1;i<=height;i++){
   m=this->subCMatrix(i);
   if(m.Arg()<=0){
    result=false;
    break;
   }
 }
 return result;
}

CMatrix  CMatrix::T()
{
 double * t=new double[width*height];
 for(long i=0;i<height;i++){
  for(long j=0;j<width;j++){
   *(t+height*j+i)=*(p+width*i+j);
  }
 }
 CMatrix m(t,height,width);
 delete [] t;
 return m;
}

CMatrix  CMatrix:: addMatrix(CMatrix & m,bool sign)
{
	if(sign == 1)
	{
     long tempW = m.getWidth();
	 long tempH = m.getHeight();
     ASSERT(height == tempH);
     long w = tempW + width;
     long h = height;
     double* t;
	 t = new double[w * h];
	 for(long i = 0; i < h; i++)
	 {
	   for(long j = 0;j < width; j++)
	   {
		   *(t + w * i + j) = *(p + width * i+j);
	   }
	   for(long k = 0;k < tempW;k++)
	   {
		   *(t + w * i + width + k) = m[i][k];
	   }
	 }
     CMatrix matrix(t,h,w);
     return matrix;
	}//(A,b)形式

	else if(sign == 0)
	{
     long tempW = m.getWidth();
	 long tempH = m.getHeight();
     ASSERT(width == tempW);
     long w = width;
     long h = height + tempH;
     double* t;
	 t = new double[w * h];
	 for(long i = 0; i < height; i++)
	 {
	  for(long j = 0; j < width; j++)
	  {
	    *(t + w * i + j) = *(p + width * i+j);
	  }
	 }

	 for(i = height; i < h; i++)
	 {
      for(long j = 0;j < width;j++)
	  {
	   	*(t + w * i + j) = m[i - height][j];
	  }
	 }
     CMatrix matrix(t,h,w);
     return matrix;   
	}//(A
//	    b)形式
}

long CMatrix::findMainElem(long row,long col)
{
 long i;
 double max;
 long n;
 max= *(p+width*row+col);
 n=row;
 for(i=row+1;i<=height-1;i++)
 {
  
  if(*(p+width*i+col)>fabs(max))
  {
   max=*(p+width*i+col);
   n=i;
  } 
 }
 return n;
}//选主元

void CMatrix::exchangeRow(long i, long j)
{
 //仅仅用于高斯消元；
 long k;
 double t;
 for(k=0;k<=width-1;k++)
 {
  t=*(p+width*i+k);
  double z1 = *(p+width*i+k)=*(p+width*j+k);
  double z2 = *(p+width*j+k)=t;
 }

}//更换ij行

CMatrix  CMatrix:: GetTrian()
{
 double k1;
 CMatrix m=*this;

 for(long row = 0,col = 0;(row<height-1)&&(col<width);)
 {
   long sign = m.findMainElem(row,col);
   if(sign != row)
   m.exchangeRow(row,sign);

   while( (m[row][col] == 0)&&(col<width) )
   {
	  col++;
	  long s = m.findMainElem(row,col);
      if(s != row)
	  m.exchangeRow(row,s);
   }
  if(col<width)
  {
    for(long j=row+1;j<height;j++)
    {
     k1=m[j][col]/m[row][col];
     m[j][col]=0;
     for(long n=col+1;n<width;n++)
     {
     m[j][n]=m[j][n]-k1*m[row][n];
	 }
    }
	row++;
	col++;
  }
  else;
 }
 return m;
}//获得上三角形阵，用选行主元消去法

CMatrix  CMatrix:: GetTopTrian()
{
 double k1;
 CMatrix m=*this;

 for(long row = 0,col = 0;(row<height-1)&&(col<width);)
 {
   while( (m[row][col] == 0)&&(col<width) )
   {
	  col++;
   }
  if(col<width)
  {
    for(long j=row+1;j<height;j++)
    {
     k1=m[j][col]/m[row][col];
     m[j][col]=0;
     for(long n=col+1;n<width;n++)
     {
     m[j][n]=m[j][n]-k1*m[row][n];
	 }
    }
	row++;
	col++;
  }
  else;
 }
 return m;
}//获得上三角形阵，用顺序消元法

CMatrix  CMatrix::GetBotTrian()
 {
 double k1;
 CMatrix m=*this;
 for(long i=height-1;i>=0;i--)
 {//Gauss消去法，变换成下三角矩阵
  for(long j=i-1;j>=0;j--)
  {
   k1=m[j][i]/m[i][i];
   for(long n=0;n<width;n++)
   {
    m[j][n]=m[j][n]-k1*m[i][n];
   }
  }
 }

 return m;
 }//获得下三角形阵，用顺序消元法

double  CMatrix:: order()
{
	double order = 0;
	CMatrix m=this->GetTrian();
    long h;
    long w;
    for(h = 0;h < m.height;h++)
	{
     double sub = 0;
     for(w = 0;w < m.width;w++)
	 {
	 sub += fabs(m[h][w]); 
	 }
	 if(sub != 0)
	  order++;
	}
	return order;
}//求秩

CMatrix  CMatrix:: rever()
{
 ASSERT(width==height);//是方阵才可以计算
 ASSERT(order() == width);//满秩才可以计算
 CMatrix e(width);
 CMatrix m = addMatrix(e,1);//增广矩阵
 CMatrix mTop = m.GetTrian();//行变换为上三角阵
 CMatrix mBot = mTop.GetBotTrian();//行变换为下三角阵
 double* t = new double[width * width];
 for(long i = 0;i < height; i++)
	 for(long j = width; j<mBot.getWidth(); j++)
	 {
	  double k = mBot[i][i];
	  mBot[i][j] /= k;
	  *(t + width * i + j - width) = mBot[i][j];
	 }
 CMatrix matrix(t,width,width);
 return matrix;
}//求逆阵

CMatrix CMatrix::operator +(CMatrix &m1)
{
 ASSERT(m1.height==height && m1.width==width);
 long tmpHeight=m1.height;
 long tmpWidth=m1.width;
 double * t=new double[tmpWidth*tmpHeight];
 for(long i=0;i<tmpHeight;i++){
  for(long j=0;j<tmpWidth;j++){
   *(t+tmpWidth*i+j)=*(m1.p+tmpWidth*i+j)+*(p+tmpWidth*i+j);
  }
 }
 CMatrix m(t,tmpHeight,tmpWidth);
 delete [] t;
 return m;
}

CMatrix CMatrix::operator -(CMatrix &m1)
{
 ASSERT(m1.height==height && m1.width==width);
 long tmpHeight=m1.height;
 long tmpWidth=m1.width;
 double * t=new double[tmpWidth*tmpHeight];
 for(long i=0;i<tmpHeight;i++){
  for(long j=0;j<tmpWidth;j++){
   *(t+tmpWidth*i+j)=*(p+tmpWidth*i+j)-*(m1.p+tmpWidth*i+j);
  }
 }
 CMatrix m(t,tmpHeight,tmpWidth);
 delete [] t;
 return m;
}

CMatrix CMatrix:: operator-()
{
	double *t = new double[height * width];
	for(long i = 0; i < height; i++)
	 for(long j = 0; j < width; j++)
		 *(t + width * i + j ) = - *(p + width * i + j );
	 CMatrix matrix(t,height,width);
	 return matrix;
}//负号

CMatrix CMatrix::operator *(CMatrix &m1)
{
 ASSERT(m1.height==width);
 if(!this->isVector() && m1.isVector()){//左为数,右为矩阵
  CMatrix m;
   m=p[0]*m1;
  return m;
 }else if(this->isVector() && !m1.isVector()){//左为矩阵,右为数
  CMatrix m;
   m=*this*m1[0][0];
  return m;
 }else if(!this->isVector() && m1.isVector()){//左右都为数
  double * t=new double[1];
  t[0]=p[0]*m1[0][0];
  CMatrix m(t,1,1);
  delete [] t;
  return m;
 }else if(this->isVector() && m1.isVector() && width==m1.height){//左为矩阵,右为矩阵
  double sum;
  double * t=new double[height*m1.width];
  for(long i=0;i<height;i++){
   for(long j=0;j<m1.width;j++){
    sum=0;
    for(long k=0;k<width;k++){
     sum+=(*(p+width*i+k))*(m1[k][j]);
    }
    *(t+m1.width*i+j)=sum;
   }
  }
  CMatrix m(t,height,m1.width);
  delete [] t;
  return m;
 }else{
  ASSERT(0);//未知运算
  return *this;
 }
}

CMatrix operator*(double alpha,CMatrix & m1)
{
 CMatrix m=m1;
 for(long i=0;i<m.height;i++){
  for(long j=0;j<m.width;j++){
   m[i][j]=alpha*m1[i][j];
  }
 }
 return m;
}

CMatrix CMatrix::operator*(double alpha)
{
 return alpha*(*this);
}

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

CMatrix CMatrix::operator-=(CMatrix & m)
{
 return *this-m;
}

CMatrix CMatrix::operator *=(double alpha)
{
 return *this*alpha;
}

CMatrix sqrt(CMatrix m)
{
 m[0][0]=sqrt(m[0][0]);
 return m;
}

double abs(CMatrix & m)
{
 double sum=0;
 for(long i=0;i<m.height;i++){
  for(long j=0;j<m.width;j++){
   sum+=m[i][j]*m[i][j];
  }
 }
 return sqrt(sum);
}

CMatrix CMatrix::operator /(CMatrix &m1)
{
 ASSERT(m1.width==1 && m1.height==1);
 return *this/m1[0][0];
}

CMatrix CMatrix::operator /(double sub)
{
 CMatrix m=*this;
 for(long i=0;i<height;i++){
  for(long j=0;j<width;j++){
   m[i][j]=*(p+width*i+j)/sub;
  }
 }
 return m;
}

CMatrix & CMatrix::operator =(CMatrix & m)
{
 if(&m==this) return *this;
 height=m.height;
 width=m.width;
 delete [] p;
 p=new double[height*width];
 for(long i=0;i<height;i++){
    for(long j=0;j<width;j++){
     *(p+width*i+j)=*(m.p+width*i+j);
    }
  }
 return *this;
}

double * CMatrix::operator [](long heightPos)
{
 ASSERT(heightPos>=0 && heightPos<height);
 return p+width*heightPos;
}
ostream & operator<<(ostream & os,CMatrix & m)
{
 os<<"CMatrix:height="<<m.height<<endl;
 os<<"CMatrix:width="<<m.width<<endl;
 for(long i=0;i<m.height;i++){
  for(long j=0;j<m.width;j++){
   os<<m[i][j]<<" "; 
  } 
  os<<endl;
 }
 return os;
}

CMatrix::~CMatrix(void)
{
 delete [] p;
}
void CMatrix::Test(){
 double arr[4][4]={{1,1,1,1},{1,2,3,4},{1,3,6,10},{1,4,10,20}};
 CMatrix m1,m(&arr[0][0],4,4),m2(&arr[0][0],4,4);
 m.isPtv();
 /*            
 cout<<"arranger="<<m.Arg()<<endl;
 cout<<m;*/
}