goodtutorialonc++templateforbeginners.txt

给c++初学者一个用于科学计算的好例子

-Good tutorial on C++ Template for beginnersC++写的通用矩阵类CMatrix
关键词： 矩阵    矩阵运算    数值方法                                           

     在编写数值计算的时候经常会用到矩阵和矩阵相关的操作，每次都从重新写，烦不胜烦，在前些日子在编写程序计算高速飞行器再入大气层时翼粱受热的温度场时，写了一个通用的矩阵类CMatrix，封装了常用的一些算法，还有很多算法，暂时没空写。

CMatrix类的基本功能:

1.内存管理一律根据实际需要的大小在堆中动态分配内存。

2.边界检查。如果数组下标超越了数组大小界限，会给出警告信息，可以防止非法内存访问以及方便程序的调试。

3.重载了+,-,*,+=,-=,数乘等常见运算符；

4.可以保存数组为二进制数据文件和文本文件两种形式，也可以从二进制数据文件和文本文件读取数据到数组。

5.实现了和矩阵相关的线性代数方程组求解算法。一是高斯选主元消去法二是针对三对角矩阵的追赶法。

6.静态函数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  
{
 //构造函数；
public:
 CMatrix(const CMatrix& m);
 CMatrix();
 CMatrix(int row,int col);
 CMatrix(DOUBLE dbArrary[],int n,BOOL bAsRow);
 
 //运算符重载；
 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;//行数；
 int nCol;//列数；
 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_)


 

实现文件: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的值是随机值；
 //初始化之后，可以方便以后判断是否已经分配过内存；
}

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;
 }
 
}

//拷贝构造函数；
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("** 警告：数组越界! **");
  exit(0);
 }
 return Val[i][j];
}

void CMatrix::SetSize(int row, int col)
{
 //重新设置数组的大小；原有的保持不变；增加的置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)
{
 //重新调整矩阵大小；
 //如果矩阵扩大，原矩阵与新矩阵对应位置的值不变，扩大部分的值为随机值；
 //如果矩阵减小，对应位置的值保持不变；
 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("不满足数组相加的条件");
  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("数组维数不相等，不能相加");
  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("数组维数不相等，不能相减");
  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("数组维数不满足相乘条件");
  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]=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("数组维数不相等，不能用＝赋值");
  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("** 只有方阵才能转换为Diag标准形  **",MB_ICONSTOP);
  exit(0);
 }
}

//释放数组中Val被分配的内存；
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)
{ 
 //只有三对角矩阵才可以使用此方法求解；
 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++)
 {