matrix.h

矩阵计算的c++代码

/******************************************
 * Copyright (c) 2002-2004, LHD, IMech, CAS
 * All rights reserved.
 * 
 * Filename: Matrix.h
 * File title:
 * Abstract: In this file, a __MatrixType__-type matrix is defined,
 *			together wiht some functions operating the matrix.
 *			__MatrixType__ is 'double' in this file.
 
 本矩阵类实现了对矩阵的一些基本操作，比如 +, -, *，求逆等。
 实现了矩阵的列主元LU分解，并利用它实现了线性代数方程组的直接解法；还有线性方程组的迭代解法。
 另外定义了一些非数学上但经常用到的一些操作
 比如两个矩阵对应元素的相乘、相除，对矩阵的每个元素求abs,sqrt等
 *			
 * 
 * Current Version:	1.4
 * Current Author:		Liu Changli (刘长礼)
 * Modification date:	April 03, 2004

 * Version 1.4.1, On Apr. 6, 2004, modified by Liu Changli (刘长礼)
 *    add ItJacobi() & SorSolver() functions to solve linear equations
 * Version 1.4, On Apr. 3, 2004, modified by Liu Changli (刘长礼)
 *    save all the non-class function to another file(xxMatrix)
 * Version 1.3.1, On Mar. 25, 2004, modified by Liu Changli (刘长礼)
 *    change Write function and add Read() function
 * Version 1.3, On Jan. 03, 2004, modified by Liu Changli (刘长礼)
 *    change vector to CVector and reserve vector
 * Version 1.2, On Jul. 05, 2003, modified by Liu Changli (刘长礼)
 *    add LUDmcp, LUSolver, Inverse, Cond, Det.
 * Version 1.1, On Mar. 03, 2003, modified by Liu Changli (刘长礼)
 * Version 1.0, On Dec. 23, 2002, modified by Liu Changli (刘长礼)
 ********************************************/

#pragma once
#ifndef __MATRIX_H
#define __MATRIX_H

////////////////////////////////////////////////////////////////////////
// fmxl is abbreviation of 'formal matrix library'
////////////////////////////////////////////////////////////////////////

#ifndef __FMXL_BEGIN_NAMESPACE
#define __FMXL_BEGIN_NAMESPACE namespace fmxl   {
#endif
#ifndef __FMXL_END_NAMESPACE
#define __FMXL_END_NAMESPACE   }
#endif

// namespace begin
__FMXL_BEGIN_NAMESPACE

#include "MyHeader.h"
#include "Vector.h"

using namespace std;

typedef double __MatrixType__;

class CMatrix
{
public:
	typedef std::vector<__MatrixType__>	DVector;
	typedef basic_filebuf<__MatrixType__, char_traits<__MatrixType__> > _Mx_filebuf;
	typedef basic_ifstream<__MatrixType__, char_traits<__MatrixType__> > _Mx_ifstream;
	typedef basic_ofstream<__MatrixType__, char_traits<__MatrixType__> > _Mx_ofstream;
	typedef basic_fstream<__MatrixType__, char_traits<__MatrixType__> > _Mx_fstream;
public:
	// constructor and destructor
	CMatrix(void);
	virtual ~CMatrix(void);
	CMatrix(int nRows, int nCols);
	CMatrix(const CMatrix& src);
	
	void SetSize(int nRows, int nCols);	// create a new matrix
	void Clear(void); // clear the matrix

	// operation of matrix
	CMatrix& operator =  (__MatrixType__ src);	// set a __MatrixType__ to the matrix
	CMatrix& operator =  (const CMatrix& src);	// set a matrix to the matrix
	CMatrix operator~()const;	// get inverse matrix

	// implemention
	// solve linear equations
	int LuDcmp(CMatrix &lower, CMatrix &upper, vector<int>& index)const;	// make LU Decomposition
	int LuBksb(CVector &xx, const CVector &bb)const;	// LU Back Substitution 
	int LuBksb(DVector &xx, const DVector &bb)const;	// solve a linear equations using LU Decomposition
	int Sor(CVector &xx, const CVector &bb, double omega=1.0, bool bInitlized=false, int nMaxIter=100, double dError=ERROR8)const;	// solve a linear equations using sor method
	int Sor(DVector &xx, const DVector &bb, double omega=1.0, bool bInitlized=false, int nMaxIter=100, double dError=ERROR8)const;	// solve a linear equations using sor method
	int ItJacobi(CVector & xx, const CVector &bb, bool bInitlized=false, int nMaxIter=100, double dError=ERROR8)const;	// Jacobi iterative method

	CMatrix Transpose(void)const;	// get the transpose of the matrix
	CMatrix Inverse(void)const;	// get inverse matrix
	__MatrixType__ Norm(void)const; // return the matrix's Norm
	__MatrixType__ Det(void)const; // return the matrix's Determinent
	__MatrixType__ Cond(void)const; // return the matrix's Conditions
	
	//+++++++++++++++++++++++++++++++++++++++++++++++
	const CMatrix&  operator+()const;
	CMatrix& operator+=(__MatrixType__ right);
	CMatrix& operator+=(const CMatrix& right);
	//-----------------------------------------------
	CMatrix  operator-()const;
	CMatrix& operator-=(__MatrixType__ right);
	CMatrix& operator-=(const CMatrix& right);
	//***********************************************	
	CMatrix& operator*=(__MatrixType__ right);

	// get Hilbert matrix, a famous bad matrix
	CMatrix& Hilbert(int nDim);

	//eye, unit matrix
	CMatrix& Eye(void); // unit the matrix
	CMatrix& Eye(int nDim); // unit the matrix with the dimension of nDim

	//ones, all matrix elements is ONE
	CMatrix& Ones(void); // set all the elements of the matrix as 1
	CMatrix Ones(int nRows, int nCols); // set all the elements of a new matrix as 1

	//zeros, zero matrix
	CMatrix& Zeros(void); // set all the elements of the matrix as 0
	CMatrix Zeros(int nRows, int nCols);	 // set all the elements of a new matrix as 1
	
	CMatrix TriL(bool bDiagonal=true)const; // get the lower triangle matrix
	CMatrix TriU(bool bDiagonal=true)const; // get the upper triangle matrix
	
	//access
	const	__MatrixType__& operator()(int nRow, int nCol)const;
			__MatrixType__& operator()(int nRow, int nCol);
	const	__MatrixType__& At(int nRow, int nCol)const;
			__MatrixType__& At(int nRow, int nCol);
	
	bool IsEmpty(void)const;	// test whether the matrix contains no elements
	int Rows(void)const;	// get row size
	int Cols(void)const;	// get col size

// member data
private:
	int m_nRows, m_nCols;
	__MatrixType__  **m_pData;

public:
	void Print(void)const;	// print to screen
	void Write(const string & sFilename, bool bApp=true, bool bBinAscii=true)const;	// write to a file
	void Read(const string & sFilename, bool bBinAscii=true)const;	// read from a file

public:
	friend const CVector operator * (const CMatrix& left, const CVector& right);
	friend const CMatrix operator * (const CMatrix& left, const CMatrix& right);
	friend const CMatrix operator + (const CMatrix& left, const CMatrix& right);
	friend const CMatrix operator - (const CMatrix& left, const CMatrix& right);
};

CMatrix:: CMatrix(void) : m_nRows(0), m_nCols(0), m_pData(NULL)
{}
CMatrix:: ~CMatrix(void)
{	Clear();}

CMatrix:: CMatrix(int nRows, int nCols) : m_pData(NULL)
{	SetSize(nRows, nCols);	}
CMatrix:: CMatrix(const CMatrix& src) : m_pData(NULL)
{
	SetSize(src.m_nRows, src.m_nCols);
	for(int i=0; i<m_nRows; ++i)
	{	memcpy(*(m_pData+i), *(src.m_pData+i), sizeof(__MatrixType__)*m_nCols);	}
}

// inline access
inline const __MatrixType__& CMatrix::operator()(int nRow, int nCol)const
{
	assert(nRow>=0 && nRow<m_nRows && nCol>=0 && nCol<m_nCols);
	return *(*(m_pData+nRow)+nCol);
}
inline __MatrixType__& CMatrix::operator()(int nRow, int nCol)
{
	assert(nRow>=0 && nRow<m_nRows && nCol>=0 && nCol<m_nCols);
	return *(*(m_pData+nRow)+nCol);
}
inline const __MatrixType__& CMatrix::At(int nRow, int nCol)const
{
	assert(nRow>=0 && nRow<m_nRows && nCol>=0 && nCol<m_nCols);
	return *(*(m_pData+nRow)+nCol);
}
inline __MatrixType__& CMatrix::At(int nRow, int nCol)
{
	assert(nRow>=0 && nRow<m_nRows && nCol>=0 && nCol<m_nCols);
	return *(*(m_pData+nRow)+nCol);
}
inline int CMatrix::Rows(void)const
{	return m_nRows;	}

inline int CMatrix::Cols(void)const
{	return m_nCols;	}

inline bool CMatrix::IsEmpty(void)const
{	return (m_pData == NULL);	}

// non-inline access

// create the matrix
void CMatrix::SetSize(int nRows, int nCols)
{	// 分三种情况来分配内存
	assert(nRows>=0 && nCols>=0);
	if(nRows==0 || nCols==0)
	{	// clear the matrix & the Row Vector
		Clear();
	}
	else if(IsEmpty())
	{	// create a new matrix
		m_pData = new __MatrixType__ * [nRows];
		for(int i=0; i<nRows; ++i)
		{	// must use this form, never change it
			m_pData[i] = new __MatrixType__ [nCols];
		}
		m_nRows = nRows;
		m_nCols = nCols;
	}	
	else
	{	// create a new matrix and copy the orginal elements to new ones
		__MatrixType__** pData;
		pData = new __MatrixType__  * [nRows];
		for(int i=0; i<nRows; ++i)
		{	// must use this form, never change it
			pData[i] = new __MatrixType__ [nCols];
		}

		int _M = min(m_nRows, nRows), _N = min(m_nCols, nCols);
		for(i=0; i<_M; ++i)	// copy older ones to the newer
		{	memcpy(*(pData+i), *(m_pData+i), sizeof(__MatrixType__)*_N);	}
		Clear();	// clear the matrix & the Row Vector
		m_pData = pData;
		m_nRows = nRows;
		m_nCols = nCols;
	}
}

// destroy the matrix
void CMatrix:: Clear(void)
{
	for(int i=0; i<m_nRows; ++i)
	{	// free each row's Pointer
		delete [] m_pData[i];
		m_pData[i] = NULL;
	}

	delete [] m_pData; // free the header
	m_pData = NULL;
	m_nRows = 0;
	m_nCols = 0;
}

//	print to screen
void CMatrix:: Print(void)const
{
	if(IsEmpty()) 
	{	cout<<"The matrix is a NULL matrix"<<endl;	return;	}
	for(int i=0; i<m_nRows; ++i)
	{
		for(int j=0; j<m_nCols; ++j)			
		{	cout<<*(*(m_pData+i)+j)<<"  ";	}
		cout<<endl;
	}
	cout<<endl;
}

// read from a file
void CMatrix::Read(const string & sFilename, bool bBinAscii)const
{	// bBinAscii==true means using ascii
	if(IsEmpty()) 
	{	cout<<"The matrix is a NULL matrix"<<endl;	return;	}

	if(bBinAscii)	// bBinAscii==true means that ascii is used
	{
		ifstream inFMatrix;
		inFMatrix.open(sFilename.c_str());
		assert(inFMatrix);
		// seek to begin of file
		inFMatrix.seekg(ios_base::beg);

		for(int i=0; i<m_nRows; ++i)
		{
			for(int j=0; j<m_nCols; ++j)			
			{	inFMatrix>>*(*(m_pData+i)+j);	}
		}
	}
	else	// save matrix as a binary file
	{
		_Mx_ifstream inFMatrix;
		inFMatrix.open(sFilename.c_str(), ios_base::binary|ios_base::in);
		assert(inFMatrix);
		// seek to begin of file
		inFMatrix.seekg(ios_base::beg);

		for(int i=0; i<m_nRows; ++i)
		{
			for(int j=0; j<m_nCols; ++j)			
			{	inFMatrix.read( (*(m_pData+i)+j), 1 );	}
		}
	}
}
// write to a file
void CMatrix::Write(const string & sFilename, bool bApp, bool bBinAscii)const
{	// bApp==true means appending a file; bBinAscii==true means using ascii
	if(IsEmpty()) 
	{	cout<<"The matrix is a NULL matrix"<<endl;	return;	}

	if(bBinAscii)	// bBinAscii==true means that ascii is used
	{
		ofstream outFMatrix;
		if(bApp)
		{	// open a file and seek to end of file
			outFMatrix.open(sFilename.c_str(), ios_base::app);
			assert(outFMatrix);
			outFMatrix.seekp(ios_base::end);
		}
		else
		{	// open a file and seek to begin of file
			outFMatrix.open(sFilename.c_str(), ios_base::out);
			assert(outFMatrix);
			outFMatrix.seekp(ios_base::beg);	// 必须到文件头部位置，否则出错
		}

		for(int i=0; i<m_nRows; ++i)
		{
			for(int j=0; j<m_nCols; ++j)			
			{	outFMatrix<<*(*(m_pData+i)+j)<<"  ";	}
			outFMatrix<<endl;
		}
		outFMatrix<<endl;
	}
	else	// save matrix as a binary file
	{
		_Mx_ofstream outFMatrix;
		if(bApp)
		{	// open a file and seek to end of file
			outFMatrix.open(sFilename.c_str(), ios_base::binary|ios_base::app);
			assert(outFMatrix);
			outFMatrix.seekp(ios_base::end);
		}
		else
		{	// open a file and seek to begin of file
			outFMatrix.open(sFilename.c_str(), ios_base::binary|ios_base::out);
			assert(outFMatrix);
			outFMatrix.seekp(ios_base::beg);	// 必须到文件头部位置，否则出错
		}

		for(int i=0; i<m_nRows; ++i)
		{
			for(int j=0; j<m_nCols; ++j)			
			{	outFMatrix.write( (*(m_pData+i)+j), 1 );	}
		}
	}
}

// operations
CMatrix& CMatrix:: operator=(__MatrixType__ src)
{	// no m_pData is needed because this object has been constructed
	if(IsEmpty())
	{// if the source matrix is null, return a null matrix
		return (*this);
	}

	for(int i=0; i<m_nRows; ++i)
	{
		for(int j=0; j<m_nCols; ++j)
		{	*(*(m_pData+i)+j) = src;	}
	}
	return (*this);
}
CMatrix& CMatrix:: operator=(const CMatrix& src)
{	// no m_pData is needed because this object has been constructed
	assert(src.m_nRows==m_nRows && src.m_nCols==m_nCols);

	if(IsEmpty())
	{// if the source matrix is null, return a null matrix
		return (*this);
	}

	// 检查自赋值
	if(this == &src)
	{	return (*this);	}

	for(int i=0; i<m_nRows; ++i)
	{	memcpy(*(m_pData+i), *(src.m_pData+i), sizeof(__MatrixType__)*m_nCols);	}
	return (*this);
}


// transpose the matrix
CMatrix CMatrix::Transpose(void)const
{
	if(IsEmpty())
	{// if the source matrix is null, return a null matrix
		return CMatrix();
	}

	CMatrix _result(m_nCols, m_nRows);
	for(int i=0; i<m_nRows; ++i)
	{
		for(int j=0; j<m_nCols; ++j)
		{	*(*(_result.m_pData+j)+i) = *(*(m_pData+i)+j);	}
	}
	return _result;
}

// infinity norm
__MatrixType__ CMatrix::Norm(void)const
{
	if(IsEmpty())
	{// if the source matrix is null, return a minus value
		return (__MatrixType__)(-1);
	}
	__MatrixType__ _result=0, tSum=0;
	for(int i=0; i<m_nRows; ++i)
	{
		tSum = 0;
		for(int j=0; j<m_nCols; ++j)
		{	tSum += std::fabs( *(*(m_pData+i)+j) );	}
		_result = max(_result, tSum);
	}
	return _result;
}

// operator +, +=
const CMatrix& CMatrix::operator+()const
{	return (*this);	}

CMatrix& CMatrix::operator+=(__MatrixType__ right)
{	
	if(IsEmpty())
	{// if the matrix is null, return a null matrix
		return (*this);
	}

	for(int i=0; i<m_nRows; ++i)
	{
		for(int j=0; j<m_nCols; ++j)
		{	*(*(m_pData+i)+j) += right;	}
	}
	return (*this);
}
CMatrix& CMatrix::operator+=(const CMatrix& right)
{	
	assert(right.m_nRows==m_nRows && right.m_nCols==m_nCols);

	if(IsEmpty())
	{// if the matrix is null, return a null matrix
		return (*this);
	}

	for(int i=0; i<m_nRows; ++i)
	{
		for(int j=0; j<m_nCols; ++j)
		{	*(*(m_pData+i)+j) += *(*(right.m_pData+i)+j);	}
	}
	return (*this);
}


// operator -, -=
CMatrix CMatrix::operator-()const
{	// 0 - value = -value
	if(IsEmpty())
	{// if the matrix is null, return a null matrix
		return (*this);
	}

	CMatrix _result(m_nRows, m_nCols);
	for(int i=0; i<m_nRows; ++i)
	{
		for(int j=0; j<m_nCols; ++j)
		{	*(*(_result.m_pData+i)+j) = - *(*(m_pData+i)+j);	}
	}
	return _result;