class CVector;//前向声明 //矩阵类 /**//* **author:phinecos **date:7/17/2008 */ class CMatrix { public: CMatrix(unsigned int r=0,unsigned int c=0); //由矩阵的行数和列数创建矩阵类对象,并为矩阵元素分配存储空间,将矩阵初始化为单位矩阵; CMatrix(constchar* pFileName); //由矩阵存储文件名创建矩阵类对象,文件格式可参考附录,但不限于采用此格式; CMatrix(const CMatrix& m); //复制构造函数,由已有矩阵类对象创建新的矩阵类对象; unsigned int GetRowsNum() const;//获取矩阵的行数 unsigned int GetColumnsNum() const;//获取矩阵的列数 double&operator ()(unsigned int r, unsigned int c);//重载运算符(),用于提取指定行(r)列(c)的元素值 CMatrix&operator=(const CMatrix& m); //重载运算符=,用于矩阵之间相互赋值 CMatrix&operator+() const;//重载一元运算符+,即取矩阵本身 CMatrix operator-() const;//重载一元运算符-,即矩阵元素取相反数 CMatrix operator+(const CMatrix& m) const;//重载二元运算符+,即两个矩阵求和 CMatrix operator-(const CMatrix& m) const;//重载二元运算符-,即两个矩阵求差 CMatrix operator*(const CMatrix& m) const;//重载二元运算符*,即两个矩阵求积 CMatrix operator*(constdouble& x) const;//重载二元运算符*,即矩阵与数相乘 CMatrix operator*(CVector& v) const;//重载二元运算符*,即矩阵与向量相乘 CMatrix operator/(constdouble& x ) const;//重载二元运算符/,即矩阵与数相除 CMatrix operator^(constint& t) const;//重载二元运算符^,即矩阵求t次幂 voidoperator-=(const CMatrix& m); //重载二元运算符-=,即自减运算 voidoperator+=(const CMatrix& m);//重载二元运算符-=,即自加运算 voidoperator*=(constdouble x);//重载二元运算符*=,即自乘运算(矩阵) CMatrix Tranpose() const;//求矩阵的转置的函数 CMatrix Invert() const;//求矩阵的逆的函数 void Zeros();//矩阵归零化,将当前矩阵的所有元素归零 void Unit();//矩阵单元化,将当前矩阵转换为单位矩阵 int AddRow(double* pe, unsigned int nr);//在矩阵的nr行位置插入一行,数据存放地址为pe,返回行标 int AddColumn(double* pe, unsigned int nc);//在矩阵的nc列位置插入一列,数据存放地址为pe,返回列标 double* DeleteRow(unsigned int nr);//删除矩阵nr行,返回该行元素值(临时存储地址) double* DeleteColumn(unsigned int nc);//删除矩阵nc列,返回该列元素值(临时存储地址) friend class CLinearEquation;//声明友元类CLinearEquation friend class CVector;//声明友元类CVector //析构函数,释放存储矩阵元素的空间 virtual~CMatrix(void); //void printMatrix()const; private: double* pElement;//元素存储区 unsigned int nRow;//行数 unsigned int nColumn;//列数 void destroy(); void alloc(unsigned int n=0);//分配空间 protected: virtualvoid InitFromFile(constchar* pFileName);//从文件中初始化 public: virtualint MatOut(constchar* pFileName); //将矩阵以pFileName为文件名进行文件输出 };
矩阵类实现 #include "stdafx.h" #include "Matrix.h" #include "Vector.h" #include <iostream> #include <cassert> #include <fstream> usingnamespace std; CMatrix::CMatrix(const CMatrix& m) {//复制构造函数,由已有矩阵类对象创建新的矩阵类对象; unsigned int i; this->nRow = m.GetRowsNum(); this->nColumn = m.GetColumnsNum(); //this->pElement = new double[nRow*nColumn]; this->alloc(nRow*nColumn); for (i=0;i<nRow*nColumn;++i) { this->pElement[i] = m.pElement[i]; } } CMatrix::CMatrix(constchar* pFileName) {//由矩阵存储文件名创建矩阵类对象,文件格式可参考附录,但不限于采用此格式; if(pFileName!=NULL) { this->InitFromFile(pFileName); } } void CMatrix::InitFromFile(constchar* pFileName) {//从文件中初始化 ifstream inFile( pFileName ); if ( !inFile ) { cerr <<"unable to open input file: "<< pFileName <<" -- bailing out!\n"; return; } inFile>>this->nRow>>this->nColumn; this->alloc(nRow*nColumn);//分配空间 unsigned int i; double num; for(i=0;i<nRow*nColumn;++i) { inFile>>num; this->pElement[i] = num; } inFile.close(); } CMatrix::CMatrix(unsigned int r,unsigned int c):nRow(r),nColumn(c) {//由矩阵的行数和列数创建矩阵类对象,并为矩阵元素分配存储空间,将矩阵初始化为单位矩阵; if(nRow>0&&nColumn>0) { this->nRow = r; this->nColumn = c; //this->pElement = new double[nRow*nColumn]; this->alloc(nRow*nColumn); for (unsigned int i=0;i<nRow*nColumn;++i) { this->pElement[i] =0; } } } CMatrix& CMatrix::operator=(const CMatrix& m) {//重载运算符=,用于矩阵之间相互赋值 if(this!=&m) { this->destroy(); this->nRow = m.GetRowsNum(); this->nColumn = m.GetColumnsNum(); this->alloc(nRow*nColumn); for (unsigned int i=0;i<nRow*nColumn;++i) { this->pElement[i] = m.pElement[i]; } } return*this; } double& CMatrix::operator ()(unsigned int r, unsigned int c) {//重载运算符(),用于提取指定行(r)列(c)的元素值 returnthis->pElement[r*nColumn+c]; } CMatrix& CMatrix::operator+()const {//重载一元运算符+,即取矩阵本身 return const_cast<CMatrix&>(*this); } CMatrix CMatrix::operator-() const {//重载一元运算符-,即矩阵元素取相反数 CMatrix result(nRow,nColumn); for (unsigned int i=0;i<nRow*nColumn;++i) { result.pElement[i] =-pElement[i]; } return result; } CMatrix CMatrix::operator+(const CMatrix& m) const {//重载二元运算符+,即两个矩阵求和 assert(this->nRow==m.GetRowsNum() &&this->nColumn==m.GetColumnsNum()); CMatrix result(*this); for (unsigned int i=0;i<nRow*nColumn;++i) { result.pElement[i]+=m.pElement[i]; } return result; } CMatrix CMatrix::operator-(const CMatrix& m) const {//重载二元运算符-,即两个矩阵求差 assert(this->nRow==m.GetRowsNum() &&this->nColumn==m.GetColumnsNum()); CMatrix result(*this); for (unsigned int i=0;i<nRow*nColumn;++i) { result.pElement[i]-=m.pElement[i]; } return result; } CMatrix CMatrix::operator*(const CMatrix& m) const {//重载二元运算符*,即两个矩阵求积 CMatrix result(*this); int ct =0, cm =0, cw =0; // compute w(i,j) for all i and j for (unsigned int i =1; i <= nRow; ++i) { // compute row i of result for (unsigned int j =1; j <= m.GetColumnsNum(); ++j) { // compute first term of w(i,j) double sum = pElement[ct] *m.pElement[cm]; // add in remaining terms for (unsigned int k =2; k <= nColumn; ++k) { ct++; // next term in row i of *this cm += m.GetColumnsNum(); // next in column j of m sum += pElement[ct] *m.pElement[cm]; }; result.pElement[cw++] = sum; // save w(i,j) // reset to start of row and next column ct -= nColumn -1; cm = j; }// reset to start of next row and first column ct += nColumn; cm =0; } return result; } CMatrix CMatrix::operator*(constdouble& x) const {//重载二元运算符*,即矩阵与数相乘 CMatrix result(*this); for (unsigned int i=0;i<nRow*nColumn;++i) { result.pElement[i] =this->pElement[i]*x; } return result; } CMatrix CMatrix::operator*(CVector& v) const {//重载二元运算符*,即矩阵与向量相乘 assert(this->nColumn==v.GetDegree()); CMatrix result(nRow,1); int ct =0, cm =0, cw =0; for (unsigned int i=0;i<nRow;++i) { double sum =0.0f; for (unsigned int j=0;j<v.GetDegree();++j) { sum += pElement[ct++]*v[j]; } result(i,0) = sum; } return result; } CMatrix CMatrix::operator^(constint& t) const {//重载二元运算符^,即矩阵求t次幂 CMatrix result(*this); CMatrix tmp(*this); for (int i=1;i<t;++i ) { result = result*(tmp); } return result; } void CMatrix::operator-=(const CMatrix& m) {//重载二元运算符-=,即自减运算 assert(nRow==m.GetRowsNum()&&nColumn==m.GetColumnsNum()); for (unsigned int i=0;i<nRow*nColumn;++i) { this->pElement[i] -= m.pElement[i]; } } void CMatrix::operator+=(const CMatrix& m) {//重载二元运算符-=,即自加运算 assert(nRow==m.GetRowsNum()&&nColumn==m.GetColumnsNum()); for (unsigned int i=0;i<nRow*nColumn;++i) { this->pElement[i] += m.pElement[i]; } } void CMatrix::operator*=(constdouble x) {//重载二元运算符*=,即自乘运算(矩阵) for (unsigned int i=0;i<nRow*nColumn;++i) { this->pElement[i] *= x; } } CMatrix CMatrix::operator/(constdouble& x ) const {//重载二元运算符/,即矩阵与数相除 CMatrix result(*this); for (unsigned int i=0;i<nRow*nColumn;++i) { result.pElement[i] =this->pElement[i]/x; } return result; } CMatrix CMatrix::Tranpose() const {//求矩阵的转置的函数 CMatrix result(*this); for (unsigned int i=0;i<nRow;++i) { for (unsigned int j=0;j<nColumn;++j) { result(i,j) =this->pElement[j*nColumn+i]; } } return result; } CMatrix CMatrix::Invert() const {//求矩阵的逆的函数 CMatrix result(nRow,nColumn*2); unsigned int i,j,p,q; for (i=0;i<nRow;++i) { for (j=0;j<nColumn;++j) { result(i,j) =this->pElement[i*nColumn+j]; } for (j=nColumn;j<nColumn*2;++j) { if ((j-i)==nColumn) { result(i,j) =1.0f; } else { result(i,j) =0.0f; } } } for(i=0;i<nRow;++i) { if(result(i,i)!=1.0f) { double bs = result(i,i); result(i,i)=1.0f; for(j=i+1;j<nColumn*2;++j) { result(i,j)/=bs; } } for(q=0;q<nRow;++q) { if(q!=i) { double bs = result(q,i); for(p=0;p<nColumn*2;++p) { result(q,p) -= bs*result(i,p); } } else { continue; } } } return result; } void CMatrix::Zeros() {//矩阵归零化,将当前矩阵的所有元素归零 for (unsigned int i=0;i<nRow*nColumn;++i) { this->pElement[i] =0.0f; } } void CMatrix::Unit() {//矩阵单位化,将当前矩阵转换为单位矩阵 for (unsigned int i=0;i<nRow;++i) { for (unsigned int j=0;j<nColumn;++j) { if (i==j) { (*this)(i,j) =1.0f; } else { (*this)(i,j) =0.0f; } } } } int CMatrix::AddRow(double* pe, unsigned int nr) {//在矩阵的nr行位置插入一行,数据存放地址为pe,返回行标 CMatrix tmp(*this); unsigned int oldRow,oldColumn; oldRow =this->nRow; oldColumn =this->nColumn; this->destroy(); this->nRow = oldRow+1; this->nColumn = oldColumn; this->alloc(nRow*nColumn); unsigned int i,j,k; for (i=0;i<nr*nColumn;++i) { this->pElement[i] = tmp.pElement[i]; } k = i; for (j=0;j<nColumn;++j) { this->pElement[i++] = pe[j]; } for (;i<nRow*nColumn;++i) { this->pElement[i] = tmp.pElement[k++]; } return nr; } int CMatrix::AddColumn(double* pe, unsigned int nc) {//在矩阵的nc列位置插入一列,数据存放地址为pe,返回列标 CMatrix tmp(*this); unsigned int oldRow,oldColumn; oldRow =this->nRow; oldColumn =this->nColumn; this->destroy(); this->nRow = oldRow; this->nColumn = oldColumn+1; this->alloc(nRow*nColumn); unsigned int i,j,k; for (i=0;i<nc;++i) { for (j=0;j<nRow;++j) { this->pElement[j*nColumn+i] = tmp.pElement[j*oldColumn+i]; } } k = i; for (j=0;j<nRow;++j) { this->pElement[j*nColumn+i] = pe[j]; } for (i=k+1;i<nColumn;++i) { for (j=0;j<nRow;++j) { this->pElement[j*nColumn+i] = tmp.pElement[j*oldColumn+k]; } k++; } return nc; } double* CMatrix::DeleteRow(unsigned int nr) {//删除矩阵nr行,返回该行元素值(临时存储地址) double* pTmp =newdouble[nColumn]; CMatrix tmp(*this); unsigned int oldRow,oldColumn; oldRow =this->nRow; oldColumn =this->nColumn; this->destroy(); this->nRow = oldRow-1; this->nColumn = oldColumn; this->alloc(nRow*nColumn); unsigned int i,j,ct=0; for(i=0;i<nr;++i) { for (j=0;j<nColumn;++j) { this->pElement[ct] = tmp(i,j); ct++; } } for (j=0;j<nColumn;++j) { pTmp[j] = tmp(nr,j); } for (i=nr+1;i<oldRow;++i) { for (j=0;j<nColumn;++j) { this->pElement[ct] = tmp(i,j); ct++; } } return pTmp; } double* CMatrix::DeleteColumn(unsigned int nc) {//删除矩阵nc列,返回该列元素值(临时存储地址) double* pTmp =newdouble[nRow]; CMatrix tmp(*this); unsigned int oldRow,oldColumn; oldRow =this->nRow; oldColumn =this->nColumn; this->destroy(); this->nRow = oldRow; this->nColumn = oldColumn-1; this->alloc(nRow*nColumn); unsigned int i,j,ct=0,cw=0; for (i=0;i<nRow;++i) { for (j=0;j<oldColumn;++j) { if(j!=nc) { this->pElement[ct] = tmp(i,j); ct++; } else { pTmp[cw] = tmp(i,j); cw++; } } } return pTmp; } unsigned int CMatrix::GetRowsNum() const {//获取矩阵的行数 returnthis->nRow; } unsigned int CMatrix::GetColumnsNum() const {//获取矩阵的列数 returnthis->nColumn; } void CMatrix::alloc(unsigned int n/**//* =0 */) { this->pElement =newdouble[n]; } void CMatrix::destroy() { if (this->pElement!=NULL) { delete[]this->pElement; this->pElement = NULL; this->nColumn =0; this->nRow =0; } } CMatrix::~CMatrix(void) { this->destroy(); } int CMatrix::MatOut(constchar* pFileName) {//将矩阵以pFileName为文件名进行文件输出 if(pFileName!=NULL) { ofstream outFile( pFileName ); if ( !outFile ) { cerr <<"unable to open output file: "<< pFileName <<" -- bailing out!\n"; return-1; } outFile<<this->nRow<<""<<this->nColumn<<endl; unsigned int i; for(i=0;i<nRow*nColumn;++i) { if((i+1)%nColumn==0) { outFile<<this->pElement[i]<<endl; } else outFile<<this->pElement[i]<<""; } outFile.close(); } return0; }
#include "matrix.h" class CVector; //线性方程组类 /**//* **author:phinecos **date:7/17/2008 */ class CLinearEquation :public CMatrix { public: CLinearEquation(void); CLinearEquation(CMatrix& coe, CVector& con );//通过系数矩阵和常数向量创建线性方程组 CLinearEquation(constchar* pFile);//通过数据文件创建线性方程组 CLinearEquation(unsigned int ne, unsigned int nv);//通过方程个数和未知数创建线性方程组 ~CLinearEquation(void); public: int AddVariable(double*pcoe, unsigned int pc);//为方程组增加一个变量,其系数矩阵增加一列,pcoe为增加系数地址,pc为增加系数列的序号,为最前面,默认追加在尾部,成功返回未知数个数,否则为-1 int AddEquation(double*pcoe, double con, unsigned int pr);//为方程组增加一个方程,其系数矩阵增加一行,pcoe为增加系数地址,con为增加方程的常数项,pr为增加系数行的序号,为最前面,默认追加在尾部,成功返回方程个数,否则为-1 int DeleteVariable(unsigned int pc);//删除方程组中第pc个未知数,并删除其系数列 int DeleteEquation(unsigned int pr);//删除方程组中第pr个方程 CVector Gaussian();//高斯消元法解线性方程组 int CLinearEquation::MatOut(constchar* pFileName); public: CVector* Constant;//方程组的常数向量,新增成员 //方程组系数矩阵,继承成员,不可见 //方程组中方程个数,继承成员,矩阵行数,不可见 //方程组中未知数个数,继承成员,矩阵列数,不可见 private: void InsertConstant(double num,int pos); void DeleteConstant(unsigned int pos); void InitFromFile(constchar* pFileName); };
线性方程组类实现 #include "stdafx.h" #include "LinearEquation.h" #include "Vector.h" #include <cmath> #include <iostream> #include <fstream> usingnamespace std; /**//* **author:phinecos **date:7/17/2008 */ CLinearEquation::CLinearEquation(void) { } CLinearEquation::CLinearEquation(constchar* pFile) {//通过数据文件创建线性方程组 if(pFile!=NULL) { this->InitFromFile(pFile); } } void CLinearEquation::InitFromFile(constchar* pFileName) { ifstream inFile( pFileName ); if ( !inFile ) { cerr <<"unable to open input file: "<< pFileName <<" -- bailing out!\n"; return; } inFile>>this->nRow>>this->nColumn; this->alloc(nRow*nColumn);//分配空间 unsigned int i; double num; for(i=0;i<nRow*nColumn;++i) { inFile>>num; this->pElement[i] = num; } this->Constant =new CVector(nRow); for(i=0;i<nRow;++i) { inFile>>num; (*this->Constant)[i] = num; } inFile.close(); } CLinearEquation::~CLinearEquation(void) { this->destroy(); if(this->Constant!=NULL) { delete Constant; this->Constant = NULL; } } CLinearEquation::CLinearEquation(unsigned int ne, unsigned int nv) {//通过方程个数和未知数创建线性方程组 this->nRow = ne; this->nColumn = nv; this->pElement =newdouble[nRow*nColumn]; this->Constant =new CVector(nRow); } CLinearEquation::CLinearEquation(CMatrix& coe, CVector& con ) {//通过系数矩阵和常数向量创建线性方程组 this->nRow = coe.GetRowsNum(); this->nColumn = coe.GetColumnsNum(); this->pElement =newdouble[nRow*nColumn]; for(unsigned int i=0;i<nRow*nColumn;++i) { this->pElement[i] = coe.pElement[i]; } this->Constant =new CVector(con); } int CLinearEquation::AddVariable(double*pcoe, unsigned int pc) {//为方程组增加一个变量,其系数矩阵增加一列,pcoe为增加系数地址,pc为增加系数列的序号,为最前面,默认追加在尾部,成功返回未知数个数,否则为-1 if(pc>=0) this->AddColumn(pcoe,pc); else this->AddColumn(pcoe,this->nColumn);//默认追加在尾部 returnthis->nColumn; } int CLinearEquation::AddEquation(double*pcoe, double con, unsigned int pr) {//为方程组增加一个方程,其系数矩阵增加一行,pcoe为增加系数地址,con为增加方程的常数项,pr为增加系数行的序号,为最前面,默认追加在尾部,成功返回方程个数,否则为-1 if(pr>=0) { this->AddRow(pcoe,pr); this->InsertConstant(con,pr); } else { this->InsertConstant(con,this->nRow); this->AddRow(pcoe,this->nRow); } returnthis->nRow; } int CLinearEquation::DeleteVariable(unsigned int pc) {//删除方程组中第pc个未知数,并删除其系数列 this->DeleteColumn(pc); return0; } int CLinearEquation::DeleteEquation(unsigned int pr) {//删除方程组中第pr个方程 this->DeleteRow(pr); this->DeleteConstant(pr); return0; } void CLinearEquation::DeleteConstant(unsigned int pos) { CVector tmp(*this->Constant); unsigned int oldDegree =this->Constant->GetDegree(); this->Constant->destory(); this->Constant =new CVector(oldDegree-1); int i,k; for(i=0;i<(int)pos;++i) { (*this->Constant)[i] = tmp[i]; } for(k=pos+1;k<(int)oldDegree;++k,++i) { (*this->Constant)[i] = tmp[k]; } } void CLinearEquation::InsertConstant(double num,int pos) { CVector tmp(*this->Constant); unsigned int oldDegree =this->Constant->GetDegree(); this->Constant->destory(); this->Constant =new CVector(oldDegree+1); int i,k; for(i=0;i<pos;++i) { (*this->Constant)[i] = tmp[i]; } k = i; (*this->Constant)[i] = num; ++i; for(;k<(int)oldDegree;++k,++i) { (*this->Constant)[i] = tmp[k]; } } CVector CLinearEquation::Gaussian() {//高斯消元法解线性方程组 int*js,l,k,i,j,is,p,q; int n = (int)this->nRow;//方程组阶数 double d,t; js =newint[n]; l=1; for (k=0;k<=n-2;k++) { d=0.0f; for (i=k;i<=n-1;i++) for (j=k;j<=n-1;j++) { t=fabs(this->pElement[i*n+j]); if (t>d) { d=t; js[k]=j; is=i;} } if (d+1.0==1.0) l=0; else { if (js[k]!=k) { for (i=0;i<=n-1;i++) { p=i*n+k; q=i*n+js[k]; t=this->pElement[p]; this->pElement[p]=this->pElement[q]; this->pElement[q]=t; } } if (is!=k) { for (j=k;j<=n-1;j++) { p=k*n+j; q=is*n+j; t=this->pElement[p]; this->pElement[p]=this->pElement[q]; this->pElement[q]=t; } t=(*this->Constant)[k]; (*this->Constant)[k]=(*this->Constant)[is]; (*this->Constant)[is]=t; } } if (l==0) { delete [] js; return(0); } d=this->pElement[k*n+k]; for (j=k+1;j<=n-1;j++) { p=k*n+j; this->pElement[p]=this->pElement[p]/d; } (*this->Constant)[k]=(*this->Constant)[k]/d; for (i=k+1;i<=n-1;i++) { for (j=k+1;j<=n-1;j++) { p=i*n+j; this->pElement[p]=this->pElement[p]-this->pElement[i*n+k]*this->pElement[k*n+j]; } (*this->Constant)[i]=(*this->Constant)[i]-this->pElement[i*n+k]*(*this->Constant)[k]; } } d=this->pElement[(n-1)*n+n-1]; if (fabs(d)+1.0==1.0) { delete [] js; return(0); } (*this->Constant)[n-1]=(*this->Constant)[n-1]/d; for (i=n-2;i>=0;i--) { t=0.0; for (j=i+1;j<=n-1;j++) { t=t+this->pElement[i*n+j]*(*this->Constant)[j]; } (*this->Constant)[i]=(*this->Constant)[i]-t; } js[n-1]=n-1; for (k=n-1;k>=0;k--) { if (js[k]!=k) { t=(*this->Constant)[k]; (*this->Constant)[k]=(*this->Constant)[js[k]]; (*this->Constant)[js[k]]=t; } } delete [] js; CVector result(*this->Constant); return result; } int CLinearEquation::MatOut(constchar* pFileName) {//将矩阵以pFileName为文件名进行文件输出 if(pFileName!=NULL) { ofstream outFile( pFileName ); if ( !outFile ) { cerr <<"unable to open output file: "<< pFileName <<" -- bailing out!\n"; return-1; } unsigned int i; for(i=0;i<this->nRow;++i) { if(i==nRow-1) { outFile<<(*this->Constant)[i]<<endl; } else outFile<<(*this->Constant)[i]<<""; } outFile.close(); } return0; }