-Good tutorial on C++ Template for beginnersC++寫的通用矩陣類CMatrix
關(guān)鍵詞： 矩陣    矩陣運(yùn)算    數(shù)值方法                                           

     在編寫數(shù)值計(jì)算的時(shí)候經(jīng)常會(huì)用到矩陣和矩陣相關(guān)的操作，每次都從重新寫，煩不勝煩，在前些日子在編寫程序計(jì)算高速飛行器再入大氣層時(shí)翼粱受熱的溫度場(chǎng)時(shí)，寫了一個(gè)通用的矩陣類CMatrix，封裝了常用的一些算法，還有很多算法，暫時(shí)沒空寫。

CMatrix類的基本功能:

1.內(nèi)存管理一律根據(jù)實(shí)際需要的大小在堆中動(dòng)態(tài)分配內(nèi)存。

2.邊界檢查。如果數(shù)組下標(biāo)超越了數(shù)組大小界限，會(huì)給出警告信息，可以防止非法內(nèi)存訪問以及方便程序的調(diào)試。

3.重載了+,-,*,+=,-=,數(shù)乘等常見運(yùn)算符；

4.可以保存數(shù)組為二進(jìn)制數(shù)據(jù)文件和文本文件兩種形式，也可以從二進(jìn)制數(shù)據(jù)文件和文本文件讀取數(shù)據(jù)到數(shù)組。

5.實(shí)現(xiàn)了和矩陣相關(guān)的線性代數(shù)方程組求解算法。一是高斯選主元消去法二是針對(duì)三對(duì)角矩陣的追趕法。

6.靜態(tài)函數(shù)Bspline3()：3次B樣條曲線插值算法。

頭文件Matrix.h;

// Matrix.h: interface for the CMatrix class.
//
//////////////////////////////////////////////////////////////////////

#if !defined(AFX_MATRIX_H__E7F9EF3A_5FEB_41C3_B3A4_C909D4E45F02__INCLUDED_)
#define AFX_MATRIX_H__E7F9EF3A_5FEB_41C3_B3A4_C909D4E45F02__INCLUDED_

#if _MSC_VER > 1000
#pragma once
#endif // _MSC_VER > 1000

#include "math.h"
#include "afxtempl.h"

#define DOUBLE double

class CMatrix : public CObject  
{
 //構(gòu)造函數(shù)；
public:
 CMatrix(const CMatrix& m);
 CMatrix();
 CMatrix(int row,int col);
 CMatrix(DOUBLE dbArrary[],int n,BOOL bAsRow);
 
 //運(yùn)算符重載；
 DOUBLE& operator()(int i, int j);
 CMatrix operator+(const CMatrix& M);
 CMatrix operator+(const DOUBLE& x);
 void operator+=(const CMatrix& M);
 CMatrix operator-(const CMatrix& M);
 CMatrix operator-(const DOUBLE& x);
 void operator-=(const CMatrix& M);
 CMatrix operator=(const CMatrix& M);
 CMatrix operator*(const CMatrix& M);
 CMatrix operator*(const DOUBLE& x);
 
 void Diag(DOUBLE a, DOUBLE b, DOUBLE c);
 void Draw(CDC* pDC,CPoint Pos=CPoint(0,0));
 virtual ~CMatrix();
 void Serialize(CArchive& ar);
 
protected:
 int nRow;//行數(shù)；
 int nCol;//列數(shù)；
 DOUBLE** Val;
 void allocMemory(int m,int n);
 void freeMemory(void);
 void reallocm(int row,int col);
private:
 //void GetParam(DOUBLE Vx[],DOUBLE Vy[],int n,DOUBLE* h,DOUBLE* f,DOUBLE* u,DOUBLE* r,DOUBLE* ff,DOUBLE*d);
 int findMainElem(int row);
 void exchangeRow(int i,int j);
 
public:
 BOOL LoadMatrix();
 BOOL LoadMatrix(CString strPath);
 int GetColSize();
 int GetRowSize();
 void GetSize(int& m,int& n);
 void SetSize(int row,int cor);
 BOOL SaveMatrix(CString strPath,BOOL bAscii=true); 
 BOOL SaveMatrix(BOOL bAscii=true);
 
 static DOUBLE* Bspline3(DOUBLE xi[], int m,DOUBLE Vx[], DOUBLE Vy[], int n);
 static DOUBLE Bspline3(DOUBLE xi,DOUBLE Vx[],DOUBLE Vy[],int n);
 
 BOOL ReadMatrix(char c=' ');
 BOOL ReadMatrix(CString strFilePath,char c=' ');
 void chase_EQ(DOUBLE* pb);
 void GSElimin_EQ(DOUBLE* pb,DOUBLE TOL=1e-6);
 void Multiply(DOUBLE a);
 void SetValue(DOUBLE a);
 

    //DOUBLE max(DOUBLE a,DOUBLE b);
 
protected:
 
 DECLARE_SERIAL(CMatrix)
};

#endif // !defined(AFX_MATRIX_H__E7F9EF3A_5FEB_41C3_B3A4_C909D4E45F02__INCLUDED_)


 

實(shí)現(xiàn)文件:Matrix.cpp；

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

#include "stdafx.h"
#include "Matrix.h"
#include "math.h"
#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

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

void GetA(DOUBLE* u,DOUBLE* r,DOUBLE* d,CMatrix& A)
{
 int i,j;
 i=A.GetRowSize();
 
 A(0,0)=2;A(0,1)=r[1];
 for(j=1;j<=i-2;j++)
 {
  A(j,j)=2;
  A(j,j-1)=u[j+1];
  A(j,j+1)=r[j+1];
 }
 A(i-1,i-1)=2;A(i-1,i-2)=u[i];
 
}
void GetParam(DOUBLE Vx[], DOUBLE Vy[], int n, DOUBLE *h, DOUBLE *f, DOUBLE *u, DOUBLE *r, DOUBLE *ff, DOUBLE *d)
{
 
 int i;
 
 for(i=0;i<=n-2;i++)
 {
  h[i]=Vx[i+1]-Vx[i];
  f[i]=(Vy[i+1]-Vy[i])/h[i];
 }
 
 for(i=1;i<=n-2;i++)
 {
  u[i]=h[i-1]/(h[i-1]+h[i]);
  r[i]=1-u[i];
  ff[i]=(f[i-1]-f[i])/(Vx[i-1]-Vx[i+1]);
  d[i]=6*ff[i];
 }
}

CMatrix::CMatrix()
{
 nRow=0;
 nCol=0;
 Val=NULL;
 //這一句初始化非常重要，如果不初始化，Val的值是隨機(jī)值；
 //初始化之后，可以方便以后判斷是否已經(jīng)分配過內(nèi)存；
}

CMatrix::CMatrix(int row,int col)
{
 nRow=row;
 nCol=col;
 Val=new DOUBLE* [nRow];
 VERIFY(Val!=NULL);
 for(int i=0;i<nRow;i++)
 {
  Val[i]=new DOUBLE [nCol]; 
  VERIFY(Val!=NULL);
 }
}

CMatrix::CMatrix(DOUBLE dbArray[],int n,BOOL bAsRow)
{
 int i;
 if(bAsRow)
 {
  nRow=1;
  nCol=n;
  
  Val=new DOUBLE* [nRow];
  VERIFY(Val!=NULL);
  for(i=0;i<nRow;i++)
   Val[i]=new DOUBLE [nCol];
  
  for(i=0;i<=n-1;i++)
   Val[0][i]=dbArray[i];
 }
 else
 {

  nRow=n;
  nCol=1;
  
  Val=new DOUBLE* [nRow];
  for( i=0;i<nRow;i++)
   Val[i]=new DOUBLE [nCol];
  
  for(i=0;i<=n-1;i++)
   Val[i][0]=dbArray[i];
 }
}

CMatrix::~CMatrix()
{
 if(Val!=NULL)
 {
  for(int i=0;i<nRow;i++)
  {
   delete [] Val[i];
   Val[i]=NULL;
  }
  delete [] Val;
  Val=NULL;
 }
 
}

//拷貝構(gòu)造函數(shù)；
CMatrix::CMatrix(const CMatrix &m)
{
 nRow=m.nRow;
 nCol=m.nCol;
 
 Val=new DOUBLE* [nRow];
 for(int i=0;i<nRow;i++)
 {
  Val[i]=new DOUBLE [nCol];
  memcpy(Val[i],m.Val[i],nCol*sizeof(DOUBLE));
 } 
}

void CMatrix::Serialize(CArchive &ar)
{
 CObject::Serialize(ar);
 if(ar.IsStoring())
 {
  ar<<nRow<<nCol;
  for(int i=0;i<nRow;i++)
  {
   for(int j=0;j<nCol;j++)
   {
    ar<<Val[i][j];
   }
  }
 }
 else
 {
  int row,col;
  ar>>row>>col;
  SetSize(row,col);
  for(int i=0;i<row;i++)
  {
   for(int j=0;j<col;j++)
   {
    ar>>Val[i][j];
   }
  }
 }
}

void CMatrix::Draw(CDC *pDC, CPoint Pos)
{
 char ch[200];
 CString str;
 for(int i=0;i<nRow;i++)
 {
  for(int j=0;j<nCol;j++)
  {
   
   //if(fabs(Val[i][j])>1e3 || fabs(Val[i][j]<1e-3))
   //sprintf(ch,"%0.3e",Val[i][j]);
   //else
   sprintf(ch,"%8.3lf",Val[i][j]);
   pDC->TextOut(Pos.x+90*j,Pos.y+20*i,ch);
  } 
 } 
}

DOUBLE& CMatrix::operator()(int i, int j)
{
 if(i>=nRow||j>=nCol)
 {
  AfxMessageBox("** 警告：數(shù)組越界! **");
  exit(0);
 }
 return Val[i][j];
}

void CMatrix::SetSize(int row, int col)
{
 //重新設(shè)置數(shù)組的大小；原有的保持不變；增加的置0；
 int i,j,oldRow,oldCol;
 oldRow=nRow;
 oldCol=nCol;
 if(nRow!=row || nCol!=col)
  reallocm(row,col);
 
 if((row>oldRow)&&(col<=oldCol))//行增大，列減小；
 {
  for(int i=oldRow;i<row;i++)
   for(int j=0;j<col;j++)
    Val[i][j]=0.0;
 }
 if((col>oldCol)&&(row<=oldRow))//列增大，行減小；
 {
  for(i=0;i<row;i++)
   for(j=oldCol;j<col;j++)
    Val[i][j]=0.0;
 }
 
 if(row>oldRow && col>oldCol)//行，列均增大；
 {
  for(i=oldRow;i<row;i++)
  {
   for(j=0;j<col;j++)
    Val[i][j]=0.0;
  }
  for(i=0;i<oldRow;i++)
   for(j=oldCol;j<col;j++)
    Val[i][j]=0.0;
 }
}

void CMatrix::reallocm(int row, int col)
{
 //重新調(diào)整矩陣大小；
 //如果矩陣擴(kuò)大，原矩陣與新矩陣對(duì)應(yīng)位置的值不變，擴(kuò)大部分的值為隨機(jī)值；
 //如果矩陣減小，對(duì)應(yīng)位置的值保持不變；
 if((this->nRow)==row && (this->nCol)==col)
  return;
 
 DOUBLE** Val1=new DOUBLE* [row];
 for(int i=0;i<row;i++)
  Val1[i]=new DOUBLE[col];
 
 int colSize=min(nCol,col)*sizeof(DOUBLE);
 int minRow=min(nRow,row);
 for(i=0;i<minRow;i++)
  memcpy(Val1[i],Val[i],colSize);
 
 for(i=0;i<nRow;i++)
  delete [] Val[i];
 delete [] Val;
 
 nRow=row;
 nCol=col;
 Val=Val1;
 return;
}

CMatrix CMatrix::operator+(const CMatrix& M)
{
 
 if(nRow==M.nRow && nCol==M.nCol)
 {
  CMatrix s(nRow,nCol);
  for(int i=0;i<nRow;i++)
   for(int j=0;j<nCol;j++)
    s.Val[i][j]=this->Val[i][j]+M.Val[i][j]; 
   return s; 
 }
 else
 {
  AfxMessageBox("不滿足數(shù)組相加的條件");
  exit(0);
 }
}

CMatrix CMatrix::operator+(const DOUBLE& x)
{
 
 CMatrix s(nRow,nCol);
 for(int i=0;i<nRow;i++)
  for(int j=0;j<nCol;j++)
   s.Val[i][j]+=x; 
  return s;
}

void CMatrix::operator+=(const CMatrix& M)
{
 if(M.nRow==nRow && M.nCol==nCol)
 {
  for(int i=0;i<nRow;i++)
   for(int j=0;j<nCol;j++)
    Val[i][j]+=M.Val[i][j];
 }
 else
 {
  AfxMessageBox("數(shù)組維數(shù)不相等，不能相加");
  exit(0);
 }
 
}

CMatrix CMatrix::operator-(const CMatrix& M)
{
 
 if(nRow==M.nRow && nCol==M.nCol)
 {
  CMatrix s(nRow,nCol);
  for(int i=0;i<nRow;i++)
   for(int j=0;j<nCol;j++)
    s.Val[i][j]=this->Val[i][j]-M.Val[i][j]; 
   return s; 
 }
 else
 {
  int r=max(nRow,M.nRow);
  int c=max(nCol,M.nCol);
  
  CMatrix a(*this);
  CMatrix b(M);
  
  a.SetSize(r,c);
  b.SetSize(r,c);
  CMatrix s(r,c);
  for(int i=0;i<r;i++)
  {
   for(int j=0;j<c;j++)
    s.Val[i][j]=a.Val[i][j]-b.Val[i][j]; 
  }
  return s;
 }
}

CMatrix CMatrix::operator-(const DOUBLE& x)
{
 CMatrix s(nRow,nCol);
 for(int i=0;i<nRow;i++)
  for(int j=0;j<nCol;j++)
   s.Val[i][j]-=x;
  return s;
}
void CMatrix::operator-=(const CMatrix& M)
{
 if(M.nRow==nRow && M.nCol==nCol)
 {
  for(int i=0;i<nRow;i++)
   for(int j=0;j<nCol;j++)
    Val[i][j]-=M.Val[i][j];
 }
 else
 {
  AfxMessageBox("數(shù)組維數(shù)不相等，不能相減");
  exit(0);
 }
}

CMatrix CMatrix::operator*(const CMatrix& M)
{
 
 if(nCol==M.nRow)
 {
  CMatrix s(nRow,M.nCol);
  
  for(int i=0;i<nRow;i++)
  {
   for(int j=0;j<M.nCol;j++)
   {
    s.Val[i][j]=0;
    for(int k=0;k<nCol;k++)
     s.Val[i][j]+=Val[i][k]*M.Val[k][j];
   }
  }
  return s;
 }
 else
 {
  AfxMessageBox("數(shù)組維數(shù)不滿足相乘條件");
  exit(0);
 }
}


//重載矩陣乘以常數(shù)的運(yùn)算；
//注意：不能寫成“常數(shù)乘以矩陣”的形式；
//即：常數(shù)必須在矩陣后面；
CMatrix CMatrix::operator*(const DOUBLE& x)
{
 CMatrix s(nRow,nCol);
 for(int i=0;i<nRow;i++)
  for(int j=0;j<nCol;j++)
   s.Val[i][j]=Val[i][j]*x;
  return s;
  
}

CMatrix CMatrix::operator=(const CMatrix& M)
{
 
 if(this==&M)
  return *this;
 if(this->nCol==M.nCol && this->nRow==M.nRow)
 {
  for(int i=0;i<nRow;i++)
   for(int j=0;j<nCol;j++)
    Val[i][j]=M.Val[i][j];
   return *this;
   
 }
 else
 {
  AfxMessageBox("數(shù)組維數(shù)不相等，不能用＝賦值");
  exit(0);
 }
}

void CMatrix::SetValue(DOUBLE a)
{
 for(int i=0;i<nRow;i++)
  for(int j=0;j<nCol;j++)
   Val[i][j]=a;
}

void CMatrix::Multiply(DOUBLE a)
{
 for(int i=0;i<nRow;i++)
  for(int j=0;j<nCol;j++)
   Val[i][j]=Val[i][j]*a;
}

 


void CMatrix::Diag(DOUBLE a, DOUBLE b, DOUBLE c)
{
 int i;
 if(nRow==nCol)
 {
  this->Multiply(0.0);
  Val[0][0]=b;Val[0][1]=c;
  Val[nRow-1][nRow-1]=b;Val[nRow-1][nRow-2]=a;
  for(i=1;i<=nRow-2;i++)
  {
   Val[i][i]=b;
   Val[i][i-1]=a;
   Val[i][i+1]=c;
  }
 }
 else
 {
  AfxMessageBox("** 只有方陣才能轉(zhuǎn)換為Diag標(biāo)準(zhǔn)形  **",MB_ICONSTOP);
  exit(0);
 }
}

//釋放數(shù)組中Val被分配的內(nèi)存；
void CMatrix::freeMemory()
{
 if(Val!=NULL)
 {
  for(int i=0;i<nRow;i++)
  {
   delete [] Val[i];
   Val[i]=NULL;
  }
  delete [] Val;
  Val=NULL;
 }
}

void CMatrix::chase_EQ(DOUBLE *pb)
{ 
 //只有三對(duì)角矩陣才可以使用此方法求解；
 ASSERT(nRow==nCol);
 if(nRow!=nCol)
 {
  AfxMessageBox("只有方陣才能用追趕法求解");
  exit(0);
 }
 
 DOUBLE* a=new DOUBLE[nRow];
 DOUBLE* b=new DOUBLE[nRow];
 DOUBLE* c=new DOUBLE[nRow];
 
 DOUBLE* u=new DOUBLE[nRow];
 DOUBLE* l=new DOUBLE[nRow];
 
 //a[0],c[nRow-1] 不能被使用
 b[0]=Val[0][0];
 c[0]=Val[0][1];
 int i;
 for(i=1;i<nRow-1;i++)
 {
  b[i]=Val[i][i];
  a[i]=Val[i][i-1];
  c[i]=Val[i][i+1];
 }
 a[nRow-1]=Val[nRow-1][nRow-1-1];
 b[nRow-1]=Val[nRow-1][nRow-1];
 u[0]=b[0];
 
 //l[0] useless;
 for(i=1;i<nRow;i++)
 {