编写的一个矩阵类（C++描述）

头文件 Matrix.h

#ifndef _Matrix_
#define _Matrix_

#include <iostream>
#include <fstream>
#include <vector>
#include <string>
#include <sstream>
#include <cmath>
#include <iomanip>

#include "AllocateArray.h"
using namespace std;


// --------------------------------------------------------------------------

class Matrix
{
public:

	Matrix();
	Matrix(int, int);
	Matrix(const string&);
	
	Matrix Transpose();
	Matrix AlgebraicComplement(int,int);
	Matrix Triangulate();
	Matrix Inverse();
	
	double OneNorm();
	double Det();

	

	void GaussSeidel(vector<double>& X, vector<double>& B, double eps = 1e-6);
	void GaussElimination(vector<double>& X, vector<double>& B);

	~Matrix();


private:

	vector<vector<double> > vvElem_;
	int nRow_;
	int nCol_;
	bool minu_;

public:

	//"<<"overload
	inline friend std::ostream& operator << (std::ostream& out, const Matrix& rhs)
	{
		out << "Rows = " << rhs.nRow_ << "	Columns = "
			<< rhs.nCol_ << endl;
		for( int i = 0; i < rhs.nRow_; i++)
		{
			for( int j = 0; j < rhs.nCol_; j++)
			{
				out << rhs.vvElem_[i][j] << "\t";

			}
			out << endl;

		}
		return out;
	}



	//"+" overload
	Matrix operator + (const Matrix& rhs);

	//"-" overload
	Matrix operator - (const Matrix& rhs);

	//"*" overload
	Matrix operator * (const Matrix& rhs);

	//"*" overload
	Matrix operator * (const vector<double>& Vb);

	//"()" overload
	double& operator () (int i, int j);

	//"+=" overload
	void operator += (const Matrix& rhs);

	//"-=" overload
	void operator -= (const Matrix& rhs);

	//"=" overload
	Matrix& operator = (const Matrix& rhs);

};






#endif

实现文件

Matrix::Matrix():nRow_(0),nCol_(0)
{

}

Matrix::Matrix(int Row, int Col):nRow_(Row), nCol_(Col)
{
	AllocateDim(vvElem_, nRow_, nCol_);
}

Matrix::Matrix(const string& fileName)
{
	ifstream	infile(fileName.c_str(), ios::in);
	string		dataInput;
	while(getline( infile, dataInput ))
	{
		stringstream temp(dataInput);
		double d;
		vector<double> tempV;
		while ( temp >> d)
		{
			tempV.push_back(d);
		}
		vvElem_.push_back(tempV);
		this->nCol_ = tempV.size();
		tempV.clear();
	}
	this->nRow_ = vvElem_.size();
}

// --------------------------------------------------------------------------
/// Geting a transposion of a matrix
/// 
/// @return no
/// 
/// @author Wang li
/// @date 2016-08-09
/// @reviser 
Matrix Matrix::Transpose()
{
	Matrix m(this->nCol_,this->nRow_);
	for( int i = 0; i <  m.nRow_; i++)
	{
		for ( int j = 0; j < m.nCol_; j++)
		{
			m.vvElem_[i][j]=this->vvElem_[j][i];
		}
	}
	return m;
}

// --------------------------------------------------------------------------
/// Getting a algebraic complement
/// 
/// @param ii row to split
/// @param jj column to split
/// 
/// @return Matrix
/// 
/// @author Wang li
/// @date 2016-08-09
/// @reviser 
Matrix Matrix::AlgebraicComplement(int ii, int jj)
{
	Matrix m(this->nRow_-1,this->nCol_-1);
	for( int i = 0; i <  m.nRow_; i++)
	{
		for ( int j = 0; j < m.nCol_; j++)
		{
			if( i < (ii -1) && j < (jj -1))
			{
				m.vvElem_[i][j]=this->vvElem_[i][j];
			}
			else if( i >= (ii -1) && j < (jj -1))
			{
				m.vvElem_[i][j]=this->vvElem_[i+1][j];
			}
			else if( i < (ii -1) && j >= (jj -1))
			{
				m.vvElem_[i][j]=this->vvElem_[i][j+1];
			}
			else if( i >= (ii -1) && j >= (jj -1))
			{
				m.vvElem_[i][j]=this->vvElem_[i+1][j+1];
			}
		}
	}
	return m;
}

// --------------------------------------------------------------------------
/// Geting a transposion of a matrix
/// 
/// @return no
/// 
/// @author Wang li
/// @date 2016-08-09
/// @reviser 
Matrix Matrix::Triangulate()
{
	Matrix m(*this);
	m.minu_ = false;
	bool flag = true ;
	for ( int i = 1; i < nRow_ ; i++)
		if ( flag)
		{
			for ( int j = i; j < nCol_ ; j++)
				if (flag)
				{
					if ( m.vvElem_[i-1][i-1] == 0)
					{
						flag = false;
						for ( int ii = i; ii < nRow_; ii++)
						{
							if ( m.vvElem_[ii][i-1] != 0)
							{

								for(int k = 0 ; k < nCol_; k++)
								{
									double temp = m.vvElem_[ii][k];
									m.vvElem_[ii][k] = m.vvElem_[i-1][k];
									m.vvElem_[i-1][k] = temp;
								}
								m.minu_ = !(m.minu_);
								flag = true;
								break;
							}					
						}
						
						if (!flag) 
						{ 
							break;
						}
						
						double a1a2(0.0);
						a1a2 = m.vvElem_[j][i-1] / m.vvElem_[i-1][i-1];
						for(int k = 0 ; k < nCol_; k++)
						{						
							m.vvElem_[j][k] -= m.vvElem_[i-1][k] * a1a2;
						}

					}
					else
					{	
						double a1a2(0);
						a1a2 = m.vvElem_[j][i-1] / m.vvElem_[i-1][i-1];
						for(int k = 0 ; k < nCol_; k++)
						{						
							m.vvElem_[j][k] -= m.vvElem_[i-1][k] * a1a2;
						}
					}
				}

		}


		return m;
}

// --------------------------------------------------------------------------
/// Getting a inversion of a matrix
/// 
/// @return Matrix
/// 
/// @author Wang li
/// @date 2016-08-09
/// @reviser 
Matrix Matrix::Inverse()
{
	Matrix m(*this);
	double detT = m.Det();
	
	for (int i = 0; i < m.nRow_; i++)
	{
		for (int j = 0; j < m.nCol_; j++)
		{
			m.vvElem_[j][i] = this->AlgebraicComplement(i+1,j+1).Det() * pow(-1.0,i+j) / detT ;
		}
	}

	return m;
}

// --------------------------------------------------------------------------
/// Using gusss elimination to solve a linear matrix equation
/// 
/// @return no
/// 
/// @author Zhang shuai
/// @date 2016-08-09
/// @reviser 
void Matrix::GaussElimination( vector<double>& X, vector<double>& B)
{
	//The Finite Volume Method in Computational Fluid Dynamics An Advanced Introduction with OpenFOAM
	if (vvElem_.size() == X.size() && vvElem_.size() == B.size())
	{
		int N = vvElem_.size();

		for (int k = 0; k < N - 1; k++)
		{
			for (int i = k + 1; i < N; i++)
			{
				double Ratio = vvElem_[i][k] / vvElem_[k][k];

				for (int j = k + 1; j < N; j++)
				{
					vvElem_[i][j] = vvElem_[i][j] - Ratio * vvElem_[k][j];
				}

				B[i] = B[i] - Ratio * B[k];
			}
		}

		X[N - 1] = B[N - 1] / vvElem_[N - 1][N - 1];

		for (int i = N - 2; i >= 0; i--)
		{
			double Term(0);

			for (int j = i + 1; j < N; j++)
			{
				Term = Term + vvElem_[i][j] * X[j];
			}

			X[i] = (B[i] - Term) / vvElem_[i][i];

		}
	}
	else
	{
		cerr << "The input dimension of array is not the same.\nPlease Check the input array." << endl;
	}
}

// --------------------------------------------------------------------------
/// Using gusss seidel iteration method to solve a linear matrix equation
/// 
/// @return no
/// 
/// @author Huang chao
/// @date 2016-08-09
/// @reviser 
void Matrix::GaussSeidel(vector<double>& X, vector<double>& B,double eps)
{
	double error,sum,w(1.3),test(0);
	bool flag(true);

	if (vvElem_.size() == X.size() && vvElem_.size() == B.size())
	{
		int N = vvElem_.size();
		vector<double> Xold(N);
		do
		{
			for(int i = 0; i < N; i++)
			{
				Xold[i] = X[i];
				sum =0.0;
				for(int j = 0; j < N; j++)
				{
					if(j != i) {sum += vvElem_[i][j] * X[j];}
				}

				X[i] = (B[i]-sum)/vvElem_[i][i];

			}

			for(int i = 0; i < N; i++)
			{
				X[i]=w*X[i]+(1-w)*Xold[i];
			}

			for(int i = 0; i < N; i++)
			{
				error=fabs(Xold[i]-X[i])/X[i];
				if(error < eps) {flag = false;}
			}
			test++;

		} while (flag);

	}
	else
	{
		cerr << "The input dimension of array is not the same.\nPlease Check the input array." << endl;
	}


}

// --------------------------------------------------------------------------
/// Getting a determinate of a matrix
/// 
/// @return double
/// 
/// @author Wang li
/// @date 2016-08-09
/// @reviser 
double Matrix::Det()
{
	Matrix m(*this);	

	m = m.Triangulate();

	double mut(1);
	for ( int i = 0 ; i < this->nRow_; i++)
	{
		mut *= m.vvElem_[i][i];
	}
	if(m.minu_) {mut = -mut;}

	return mut;
}

// --------------------------------------------------------------------------
/// Getting the one norm of a matrix
/// 
/// @return double
/// 
/// @author Huang chao
/// @date 2016-08-09
/// @reviser 
double Matrix::OneNorm()
{	
	double mut1(0.0);
	for ( int i = 0 ; i < this->nRow_; i++)
	{
		for ( int j = 0 ; j < this->nCol_; j++)
		{
			mut1 += fabs(vvElem_[i][j]);
		}	
	}
	return mut1;
}

// --------------------------------------------------------------------------
/// overload operator + of matrix
/// 
/// @return Matrix
/// 
/// @author Wang jian
/// @date 2016-08-09
/// @reviser 
Matrix Matrix::operator + (const Matrix& rhs)
{
	Matrix mTemp(rhs);

	cout << " mTemp.vvElem_.size():\t" <<mTemp.vvElem_.size() << endl;

	for(int i = 0; i < rhs.vvElem_.size(); i++)
	{
		for(int j = 0; j < rhs.vvElem_[0].size(); j++)
		{
			mTemp.vvElem_[i][j] = vvElem_[i][j] + rhs.vvElem_[i][j];
		}
	}

	return mTemp;
}

// --------------------------------------------------------------------------
/// overload operator - of matrix
/// 
/// @return Matrix
/// 
/// @author Wang jian
/// @date 2016-08-09
/// @reviser 
Matrix Matrix::operator - (const Matrix& rhs)
{
	Matrix mTemp(rhs);

	for(int i = 0; i < rhs.vvElem_.size(); i++)
	{
		for(int j = 0; j < rhs.vvElem_[0].size(); j++)
		{
			mTemp.vvElem_[i][j] = this->vvElem_[i][j] - rhs.vvElem_[i][j];
		}
	}

	return mTemp;

}

// --------------------------------------------------------------------------
/// overload operator * of matrix
/// 
/// @return Matrix
/// 
/// @author Wang jian
/// @date 2016-08-09
/// @reviser 
Matrix Matrix::operator * (const Matrix& rhs)
{
	if(this->vvElem_[0].size() != rhs.vvElem_.size())
	{
		cerr << " The Matrix A cannot be multiplied by Matrix B" << endl;
	}
	Matrix mTemp(rhs);

	for(int i = 0; i < this->vvElem_.size(); i++)
	{
		for(int j = 0; j < rhs.vvElem_[0].size(); j++)
		{
			mTemp.vvElem_[i][j] = 0.0;
		}
	}

	for(int i = 0; i < this->vvElem_.size(); i++)
	{
		for(int j = 0; j < this->vvElem_[0].size(); j++)
		{
			for(int k = 0; k < rhs.vvElem_[0].size(); k++)
			{
				mTemp.vvElem_[i][k] += this->vvElem_[i][j] * rhs.vvElem_[j][k];
			}		
		}
	}

	return mTemp;

}

// --------------------------------------------------------------------------
/// overload operator () of matrix
/// 
/// @return Matrix
/// 
/// @author Wang li
/// @date 2016-08-09
/// @reviser 
double& Matrix::operator () (int i, int j)
{

	return this->vvElem_[i][j];

}

// --------------------------------------------------------------------------
/// overload operator += of matrix
/// 
/// @return Matrix
/// 
/// @author Wang jian
/// @date 2016-08-09
/// @reviser 
void Matrix::operator += (const Matrix& rhs)
{
	for(int i = 0; i < rhs.vvElem_.size(); i++)
	{
		for(int j = 0; j < rhs.vvElem_[0].size(); j++)
		{
			this->vvElem_[i][j] = this->vvElem_[i][j] + rhs.vvElem_[i][j];
		}
	}
}

// --------------------------------------------------------------------------
/// overload operator -= of matrix
/// 
/// @return Matrix
/// 
/// @author Wang jian
/// @date 2016-08-09
/// @reviser 
void Matrix::operator -= (const Matrix& rhs)
{
	for(int i = 0; i < rhs.vvElem_.size(); i++)
	{
		for(int j = 0; j < rhs.vvElem_[0].size(); j++)
		{
			this->vvElem_[i][j] = this->vvElem_[i][j] - rhs.vvElem_[i][j];
		}
	}
}

// --------------------------------------------------------------------------
/// overload operator = of matrix
/// 
/// @return Matrix
/// 
/// @author Wang jian
/// @date 2016-08-09
/// @reviser 
Matrix& Matrix::operator = (const Matrix& rhs)
{
	this->nRow_ = rhs.nRow_;
	this->nCol_ = rhs.nCol_;

	this->vvElem_.resize(nRow_);
	for(int i = 0; i < rhs.nRow_; i++)
	{
		this->vvElem_[i].resize(nCol_);
	}


	for(int i = 0; i < rhs.nRow_; i++)
	{

		for(int j = 0; j < rhs.nCol_; j++)
		{
			this->vvElem_[i][j] = rhs.vvElem_[i][j];
		}

	}

	this->minu_ = rhs.minu_;

	return *this;

}

Matrix::~Matrix()
{
}