matrix.h

一个矩阵类

/************************************************************************/
/*file name  :    matrix.h                                              */
/*function   :    basic operation with regard to matrix                 */
/*author     :    daihaijun                                             */
/*version    :    version 1.0                                           */
/*finish time:    2007.04.29                                            */
/************************************************************************/

#ifndef _MATRIX_H_
#define _MATRIX_H_

#include <math.h>
#include <assert.h>
#include <iostream>
using namespace std;

template<class T> class matrix
{
private:
	int m_nRows;
	int m_nCols;
	T** m_pData;//data space;
public:
	matrix();
	matrix(int nRows,int nCols);
	matrix(const matrix& Right);
	explicit matrix(int nSteps);//construct a nSteps-rank identity matrix
	~matrix();

public:
	friend const matrix<T> operator+ (const matrix<T>& Left,const matrix<T>& Right);
	friend const matrix<T> operator+ (const matrix<T>& Left,T Right);
	friend const matrix<T> operator+ (T Left,const matrix<T>& Right);
	friend const matrix<T> operator- (const matrix<T>& Left,const matrix<T>& Right);
	friend const matrix<T> operator- (const matrix<T>& Left,T Right);
	friend const matrix<T> operator- (T Left,const matrix<T>& Right);
	friend const matrix<T> operator* (const matrix<T>& Left,const matrix<T>& Right);
	friend const matrix<T> operator* (const matrix<T>& Left,T Right);
	friend const matrix<T> operator* (T Left,const matrix<T>& Right);
	friend ostream& operator<<(ostream& output,const matrix<T>& Obj);
	friend istream& operator>>(istream& input,matrix<T>& Obj);

	const matrix<T>& operator= (const matrix<T>& Right);
	//const matrix<T>& operator- ();
	bool operator== (const matrix<T> Right) const;
	bool operator!= (const matrix<T> Right) const;
	T* operator[] (int nRows);//read & write
	const T* operator[] (int nRows) const;//read only
	bool RowsSwap(int nk,int nl);
	bool ColsSwap(int nk,int nl);
	bool Inverse();
	void Tranpose();
private:
	void DestroyObject();

public:
	int GetRows(void) const{ return m_nRows };
	int GetCols(void) const{ return m_nCols };
	void Export(T** pData) const;//export data to 2D Array.memory must be allocated to parameter by client first
	void Import(const T** pData);//import data form 2D Array

public:
	static void Make2DArray(T** &Array,int nRows,int nCols);
	static void Delete2DArray(T** &Array,int nRows);
};

#define TMPLT(T) template <class T>

TMPLT(T) matrix<T>::matrix(void):m_nRows(0),m_nCols(0),m_pData(NULL)
{
}

TMPLT(T) matrix<T>::matrix(int nRows,int nCols):m_nRows(nRows),m_nCols(nCols)
{
	matrix<T>::Make2DArray(m_pData,nRows,nCols);
	for( int i=0;i<m_nRows;i++ )
	{
		for( int j=0;j<m_nCols;j++ )
		{
			m_pData[i][j] = 0;
		}
	}
}

TMPLT(T) matrix<T>::~matrix()
{
	DestroyObject();
}

TMPLT(T) matrix<T>::matrix(const matrix<T>& Right)
{
	m_nRows = Right.m_nRows;
	m_nCols = Right.m_nCols;
	matrix<T>::Make2DArray(m_pData,m_nRows,m_nCols);
	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];
		}
	}
}

TMPLT(T) matrix<T>::matrix(int nSteps)
{
	m_nRows = nSteps;
	m_nCols = nSteps;
	matrix<T>::Make2DArray(m_pData,m_nRows,m_nCols);
	for( int i=0;i<m_nRows;i++ )
	{
		for( int j=0;j<m_nCols;j++ )
		{
			if(i==j)
			{
				m_pData[i][j] = 1;
			}
			else
			{
				m_pData[i][j] = 0;
			}
		}
	}
}

TMPLT(T) void matrix<T>::DestroyObject()
{
	if( m_pData!=NULL )
	{
		matrix<T>::Delete2DArray(m_pData,m_nRows);
	}
	m_pData = NULL;
	m_nRows = 0;
	m_nCols = 0;
}

TMPLT(T) void matrix<T>::Make2DArray(T** &Array,int nRows,int nCols)
{
	//ceate rows pointer
	Array = new T* [nRows];
	//allocate memery for evry row
	for( int i=0;i<nRows;i++ )
	{
		Array[i] = new T [nCols];
	}
}

TMPLT(T) void matrix<T>::Delete2DArray(T** &Array,int nRows)
{
	for( int i=0;i<nRows;i++ )
	{
		//free memery for evry row
		delete [] Array[i];
	}
	//destroy row pointer
	delete [] Array;
	Array = NULL;
}

TMPLT(T) istream& operator>>(istream& input,matrix<T>& Obj)
{
	input>>ws;//Skips white space in the stream
	for(int i=0;i<Obj.m_nRows;i++)
	{
		for(int j=0;j<Obj.m_nCols;j++)
		{
			input>>Obj.m_pData[i][j];
		}
	}
	return input;
}

TMPLT(T) ostream& operator<<(ostream& output,const matrix<T>& Obj)
{
	for(int i=0;i<Obj.m_nRows;i++)
	{
		for(int j=0;j<Obj.m_nCols;j++)
		{
			output<<Obj.m_pData[i][j]<<",";
		}
		output<<endl;
	}
	return output;
}

TMPLT(T) const matrix<T> operator+ (const matrix<T>& Left,const matrix<T>& Right)
{
	assert(Left.m_nRows==Right.m_nRows && Left.m_nCols==Right.m_nCols);
	matrix<T> Result(Left.m_nRows,Left.m_nCols);
	for( int i=0;i<Left.m_nRows;i++ )
	{
		for( int j=0;j<Left.m_nCols;j++ )
		{
			Result.m_pData[i][j] = Left.m_pData[i][j] + Right.m_pData[i][j];
		}
	}
	return Result;
}

TMPLT(T) const matrix<T> operator+ (const matrix<T>& Left,T Right)
{
    matrix<T> Result(Left.m_nRows,Left.m_nCols);
	for( int i=0;i<Left.m_nRows;i++ )
	{
		for( int j=0;j<Left.m_nCols;j++ )
		{
			Result.m_pData[i][j] = Left.m_pData[i][j] + Right;
		}
	}
	return Result;
}

TMPLT(T) const matrix<T> operator+ (T Left,const matrix<T>& Right)
{
	return Right+Left;
}

TMPLT(T) const matrix<T> operator- (const matrix<T>& Left,const matrix<T>& Right)
{
	assert(Left.m_nRows==Right.m_nRows && Left.m_nCols==Right.m_nCols);
	matrix<T> Result(Left.m_nRows,Left.m_nCols);
	for( int i=0;i<Left.m_nRows;i++ )
	{
		for( int j=0;j<Left.m_nCols;j++ )
		{
			Result.m_pData[i][j] = Left.m_pData[i][j] - Right.m_pData[i][j];
		}
	}
	return Result;
}

TMPLT(T) const matrix<T> operator- (const matrix<T>& Left,T Right)
{
	Right = -Right;
	return Left+Right;
}

TMPLT(T) const matrix<T> operator- (T Left,const matrix<T>& Right)
{
	matrix<T> Result(Right.m_nRows,Right.m_nCols);
	for( int i=0;i<Right.m_nRows;i++ )
	{
		for( int j=0;j<Right.m_nCols;j++ )
		{
			Result.m_pData[i][j] = Left - Right.m_pData[i][j];
		}
	}
	return Result;
}

TMPLT(T) const matrix<T> operator* (const matrix<T>& Left,const matrix<T>& Right)
{
	assert(Left.m_nCols==Right.m_nRows);
	matrix<T> Result(Left.m_nRows,Right.m_nCols);
	int i,j;
	int p = Left.m_nCols/2;
	T* t1;
	t1 = new T [Left.m_nRows];//存放乘数每一行的奇偶积
	//计算乘数奇偶积
	for(i=0;i<Left.m_nRows;i++)
	{
		t1[i] = 0;
		for(j=0;j<p;j++)
		{
			t1[i] += Left.m_pData[i][2*j] * Left.m_pData[i][2*j+1];
		}
	}
	T* t2;
	t2 = new T [Right.m_nCols];//存放被乘数每一列的奇偶积
	//计算被乘数奇偶积
	for(i=0;i<Right.m_nCols;i++)
	{
		t2[i] = 0;
		for(j=0;j<p;j++)
		{
			t2[i] += Right.m_pData[2*j][i] * Right.m_pData[2*j+1][i];
		}
	}
	
	//计算结果
	for(i=0;i<Left.m_nRows;i++)
	{
		for(j=0;j<Right.m_nCols;j++)
		{
			Result.m_pData[i][j] = 0;
			for(int k=0;k<p;k++)
			{
				Result.m_pData[i][j] += (Left.m_pData[i][2*k]+Right.m_pData[2*k+1][j])
					*(Left.m_pData[i][2*k+1]+Right.m_pData[2*k][j]);

			}
			Result.m_pData[i][j] = Result.m_pData[i][j] - t1[i] - t2[j];
		}
	}
	
	delete [] t1;
	delete [] t2;
	
	if( Left.m_nCols%2 )//如果m_nCols为奇数
	{
		for(i=0;i<Left.m_nRows;i++)
		{
			for(j=0;j<Right.m_nCols;j++)
			{
				Result.m_pData[i][j] += Left.m_pData[i][Left.m_nCols-1]*Right.m_pData[Left.m_nCols-1][j];
			}
		}
	}	
	return Result;
}

TMPLT(T) const matrix<T> operator* (const matrix<T>& Left,T Right)
{
	matrix<T> Result(Left.m_nRows,Left.m_nCols);
	for(int i=0;i<Left.m_nRows;i++)
	{
		for(int j=0;j<Left.m_nCols;j++)
		{
			Result.m_pData[i][j] = Right*Left.m_pData[i][j];
		}
	}
	return Result;
}

TMPLT(T) const matrix<T> operator* (T Left,const matrix<T>& Right)
{
	return Right*Left;
}

TMPLT(T) const matrix<T>& matrix<T>::operator= (const matrix<T>& Right)
{
	if( this->m_pData!=NULL||m_nRows!=Right.m_nRows||m_nCols!=Right.m_nCols )
	{
		DestroyObject();
		m_nCols = Right.m_nCols;
		m_nRows = Right.m_nRows;
		Make2DArray(m_pData,m_nRows,m_nCols);
	}
	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;
}

TMPLT(T) bool matrix<T>::operator ==(const matrix<T> Right) const
{
	if( (m_nRows!=Right.m_nRows)||(m_nCols!=Right.m_nCols) )
	{
		return false;
	}
	for(int i=0;i<m_nRows;i++)
	{
		for(int j=0;j<m_nCols;j++)
		{
			if( m_pData[i][j] != Right.m_pData[i][j] )
			{
				return false;
			}
		}
	}
	return true;
}

TMPLT(T) bool matrix<T>::operator !=(const matrix<T> Right) const
{
	return !(*this==Right);
}

TMPLT(T) T* matrix<T>::operator[] (int nRows)
{
	if( nRows<m_nRows )
	{
		return m_pData[nRows];
	}
	else
	{
		return NULL;
	}
}

TMPLT(T) const T* matrix<T>::operator[] (int nRows) const
{
	if( nRows<m_nRows )
	{
		return m_pData[nRows];
	}
	else
	{
		return NULL;
	}
}

TMPLT(T) bool matrix<T>::RowsSwap(int nk,int nl)
{
	if( nk<m_nRows&&nl<m_nRows )
	{
		T m;
		for(int i=0;i<m_nCols;i++)
		{
			m = m_pData[nk][i];
			m_pData[nk][i] = m_pData[nl][i];
			m_pData[nl][i] = m;
		}
		return true;
	}
	else 
		return false;
}

TMPLT(T) bool matrix<T>::ColsSwap(int nk,int nl)
{
	if( nk<m_nCols&&nl<m_nCols )
	{
		T m;
		for(int i=0;i<m_nRows;i++)
		{
			m = m_pData[i][nk];
			m_pData[i][nk] = m_pData[i][nl];
			m_pData[i][nl] = m;
		}
		return true;
	}
	else
		return false;
}

TMPLT(T) bool matrix<T>::Inverse()
{
	bool bFlag = true;
	if( m_nRows!=m_nCols )
	{
		bFlag = false;
	}
	else
	{
		int i,j,k;//循环控制变量
		int n = m_nRows;
		int *IS;//保存行交换信息
		int *JS;//保存列交换信息
		IS = new int [n];
		JS = new int [n];
		for(k=0;k<n;k++)
		{
			//全选主元
			T Max = 0;
			for(i=k;i<n;i++)
			{
				for(j=k;j<n;j++)
				{
					if( fabs(m_pData[i][j])>Max )
					{
						Max = fabs(m_pData[i][j]);
						IS[k] = i;
						JS[k] = j;
					}
				}
			}
			if( Max+1.0==1.0 )
			{
				bFlag = false;
				return bFlag;//是奇异矩阵
			}
			if( IS[k]!=k )
			{
				RowsSwap(k,IS[k]);//交换两行
			}
			if( JS[k]!=k )
			{
				ColsSwap(k,JS[k]);//交换两列
			}
			m_pData[k][k] = 1/m_pData[k][k];
			for(j=0;j<n;j++)
			{
				if(j!=k)
				{
					m_pData[k][j] = m_pData[k][j]*m_pData[k][k];
				}
			}
			for(i=0;i<n;i++)
			{
				if(i!=k)
				{
					for(j=0;j<n;j++)
					{
						if(j!=k)
						{
							m_pData[i][j] = m_pData[i][j] - m_pData[i][k]*m_pData[k][j];
						}
					}
				}
			}
			for(i=0;i<n;i++)
			{
				if(i!=k)
				{
					m_pData[i][k] = -m_pData[i][k]*m_pData[k][k];
				}
			}
		}
		for(k=n-1;k>=0;k--)
		{
			if( JS[k]!=k )
			{
				RowsSwap(k,JS[k]);//交换两行
			}
			if( IS[k]!=k )
			{
				ColsSwap(k,IS[k]);
			}
		}
		delete [] IS;
		delete [] JS;
	}
	return bFlag;
}

TMPLT(T) void matrix<T>::Tranpose()
{
	int nNewRows = m_nCols;
	int nNewCols = m_nRows;
	T** pNewData = NULL;
	matrix<T>::Make2DArray(pNewData,nNewRows,nNewCols);
	for(int i=0;i<nNewRows;i++)
	{
		for(int j=0;j<nNewCols;j++)
		{
			pNewData[i][j] = m_pData[j][i];
		}
	}
	matrix<T>::Delete2DArray(m_pData,m_nRows);
	m_pData = pNewData;
	m_nCols = nNewCols;
	m_nRows = nNewRows;
}

TMPLT(T) void matrix<T>::Export(T** pData) const
{
	assert(pData!=NULL);
	for(int i=0;i<m_nRows;i++)
	{
		for(int j=0;j<m_nCols;j++)
		{
			pData[i][j] = m_pData[i][j];
		}
	}
}

TMPLT(T) void matrix<T>::Import(const T** pData)
{
	assert(pData!=NULL);
	for(int i=0;i<m_nRows;i++)
	{
		for(int j=0;j<m_nCols;j++)
		{
			m_pData[i][j] = pData[i][j];
		}
	}
}


#endif