matrix.h

C++语言的复数与矩阵运算库

#include "afxwin.h"
#include "strstrea.h"
#include "fstream.h"
#include "iomanip.h"
#include "stdlib.h"
#include "malloc.h"
#include "math.h"

#define CDMatrix (CMatrix<double>)
#define CIMatrix (CMatrix<int>)

template<class Type>
class CMatrix :public CObject
{
	protected:		
			Type    **f;
			int     row;
			int     col;
	public:
		int  GetRow(){return row;}
		int  GetCol(){return col;}
		void SetRow(int r){row=r;}
		void SetCol(int l){col=l;}
		Type** GetF(){return f;}		
		CMatrix();
		CMatrix(int,int);
		CMatrix(CMatrix &copy);
		~CMatrix();
		BOOL Ones(int);
		BOOL Ones(int,int);
		BOOL Zeros(int);
		BOOL Zeros(int,int);
		CMatrix Turn();
		CMatrix Inverse();
		CMatrix Construct(CMatrix &mat,CString style);
		CMatrix GetChild(int,int,int,int);
		BOOL Diag(Type *,int);
		BOOL Malloc(int,int);
		CMatrix   operator*(CMatrix  &m);
		friend CMatrix  operator*(double,CMatrix);
		friend CMatrix  operator*(CMatrix,double);
		CMatrix   operator/(Type x);
		CMatrix   operator+(CMatrix  &m);
		CMatrix   operator+(Type x);
		CMatrix   operator-(CMatrix  &m);
		CMatrix   operator-(Type x);
		CMatrix   operator-();
		Type* &operator[](int);
		friend istream  &operator>>(istream  &,CMatrix  &);
		friend ostream  &operator<<(ostream  &,CMatrix  &);
		void operator=(const CMatrix  &m);
		void operator=(Type);
		Type ** MatrixAlloc(int,int);
		void    MatrixFree(Type **);
		void    Hessenberg();
		BOOL    Cholesky(); 
		CMatrix QREigen(double,int);// return eigenvalue in two row matrix of a real matrix 		
		CMatrix Jacobi(double eps,int jt);
};

template<class Type>
Type* &CMatrix<Type>::operator[](int i)
{
	Type *tp;
	tp=(Type *)f[i];
	return tp;
}
template<class Type>
CMatrix<Type>::CMatrix()
{
	f=NULL;
	row=0;
	col=0;
}
template<class Type>
CMatrix<Type>::CMatrix(int m,int n)
{	
	f=MatrixAlloc(m,n);
	row=m;
	col=n;  
}
template<class Type>
CMatrix<Type>::~CMatrix()
{
	if(f!=NULL) MatrixFree(f);
}

template<class Type>
BOOL CMatrix<Type>::Ones(int m)
{ 	
	f=MatrixAlloc(m,m);
	for(int i=0;i<m;i++)
	for(int j=0;j<m;j++)
	{
		if(i!=j)f[i][j]=0;
		else f[i][j]=1;
	}
	row=m;
	col=m;
	return TRUE;
}
template<class Type>
BOOL CMatrix<Type>::Ones(int m,int n)
{ 	
	f=MatrixAlloc(m,n);
	for(int i=0;i<m;i++)
	for(int j=0;j<n;j++)
	{
		f[i][j]=1;
	}
	row=m;
	col=n;
	return TRUE;
}
template<class Type>
BOOL CMatrix<Type>::Zeros(int m)
{ 	
	f=MatrixAlloc(m,m);
	for(int i=0;i<m;i++)
	for(int j=0;j<m;j++)
	{
		f[i][j]=0;
	}
	row=m;
	col=m;
	return TRUE;
}
template<class Type>
BOOL CMatrix<Type>::Zeros(int m,int n)
{ 	
	f=MatrixAlloc(m,n);
	if(f==NULL) return FALSE;
	for(int i=0;i<m;i++)
	for(int j=0;j<n;j++)
	{
		f[i][j]=0;
	}
	row=m;
	col=n;
	return TRUE;
}

template<class Type>
CMatrix<Type> CMatrix<Type>::operator*(CMatrix<Type>  &mat)
{
	int l,m,n;
	l=GetRow();
	m=GetCol();
	n=mat.GetCol();
	CMatrix tmp(l,n);
	for(int k=0;k<l;k++)
	for(int j=0;j<n;j++)
	{
		for(int i=0;i<m;i++)
		tmp.f[k][j]+=f[k][i]*mat.f[i][j];
	}	
	return tmp;
}

template<class Type>
CMatrix<Type> CMatrix<Type>::operator/(Type x)
{
	int l,n;
	l=GetRow();
	n=GetCol();
	CMatrix<Type> tmp(l,n);
	for(int j=0;j<l;j++)
	for(int i=0;i<n;i++)
		tmp.f[j][i]=f[j][i]/x;
	return tmp;
}

template<class Type>
CMatrix<Type> CMatrix<Type>::operator+( CMatrix<Type>  &mat)
{
	int l1,m1;
	l1=GetRow();
	m1=GetCol();
	CMatrix<Type> tmp(l1,m1);
	for(int j=0;j<l1;j++)
	for(int i=0;i<m1;i++)
		tmp.f[j][i]=f[j][i]+mat.f[j][i];
	return tmp;
}

template<class Type>
CMatrix<Type> CMatrix<Type>::operator+(Type x)
{
	int l,m;
	l=GetRow();
	m=GetCol();
	CMatrix<Type> tmp(l,m);
	for(int j=0;j<l;j++)
	for(int i=0;i<m;i++)
		tmp.f[j][i]=f[j][i]+x;
	return tmp;
}

template<class Type>
CMatrix<Type> CMatrix<Type>::operator-(CMatrix<Type>  &mat)
{
	int l1,m1;
	l1=GetRow();
	m1=GetCol();
	CMatrix<Type> tmp(l1,m1);
	for(int j=0;j<l1;j++)
	for(int i=0;i<m1;i++)
		tmp.f[j][i]=f[j][i]-mat.f[j][i];
	return tmp;
}
template<class Type>
CMatrix<Type> CMatrix<Type>::operator-()
{
	for(int i=0;i<row;i++)
	for(int j=0;j<col;j++)
		f[i][j]=-f[i][j];
	return *this;
}	

template<class Type>
CMatrix<Type> CMatrix<Type>::operator-( Type x)
{
	int l,m;
	l=GetRow();
	m=GetCol();
	CMatrix<Type> temp(l,m);
	for(int j=0;j<l;j++)
	for(int i=0;i<m;i++)
		temp.f[j][i]=f[j][i]-x;
	return temp;
}

template<class Type>
istream &operator>>(istream &in,CMatrix<Type> &mat)
{    
	int i,j;
	for(i=0;i<mat.row;i++)
	for(j=0;j<mat.col;j++)
		 in>>mat.f[i][j];
	cout<<endl;
	return in;
};

template<class Type>
ostream  & operator<<(ostream& out,CMatrix<Type> &v1)
{
	if(v1.GetF()==NULL) 
	{
		out<<"This Matrix cannot be output!";
		return out<<endl;
	}
	out<<setiosflags(ios::right||ios::fixed);
	out<<setprecision(5); 
	out<<"["<<endl;
	int mr=v1.GetRow();
	int mc=v1.GetCol();
	for(int j=0;j<mr;j++)
	{
		for(int i=0;i<mc-1;i++)
			out<<setw(12)<<v1.f[j][i];
		out<<setw(12)<<v1.f[j][mc-1]<<";"<<endl;
	}
	return out<<"]"<<endl;
	
}
template<class Type>
CMatrix<Type>::CMatrix(CMatrix<Type> &copy)
{
	f=NULL;
	*this=copy;
}
	
template<class Type>
void CMatrix<Type>::operator=(const CMatrix<Type> &copy)
{
	
	if(this==&copy) return;
	if(f!=NULL) MatrixFree(f);
	int m,n;
	m=copy.row;
	n=copy.col;
	f=MatrixAlloc(m,n);
  	for(int i=0;i<m;i++)
	for(int j=0;j<n;j++)
		f[i][j]=copy.f[i][j];
	row=m;
	col=n;  
}

template<class Type>
void CMatrix<Type>::operator=(Type d)
{
	
	if(f==NULL) 
	{
		cout<<"The Matrix is not be allociated"<<endl;
		return;
	}
	for(int i=0;i<row;i++)
	for(int j=0;j<col;j++)
		f[i][j]=d;
}

template<class Type>
Type ** CMatrix<Type>::MatrixAlloc(int r,int c)
{
    Type *x,**y;
    int   n;
    x=(Type *)calloc(r*c,sizeof(Type));
    y=(Type **)calloc(r,sizeof(Type *));
	for(n=0;n<=r-1;++n)
	y[n]=&x[c*n];
	return(y);
}

template<class Type>
void CMatrix<Type>::MatrixFree(Type **x)
{
   free(x[0]);
   free(x);
}
template<class Type>
BOOL CMatrix<Type>::Malloc(int m,int n)
{
	f=MatrixAlloc(m,n);
	if(f==NULL) return FALSE;
	row=m;
	col=n;
	return TRUE;
}

	
template<class Type>
CMatrix<Type> CMatrix<Type>::Turn()
{
	CMatrix<Type> tmp(col,row);
	for(int i=0;i<row;i++)
	for(int j=0;j<col;j++)
		tmp.f[j][i]=f[i][j];
	return tmp;
}

template<class Type>
BOOL CMatrix<Type>::Diag(Type *array,int m)
{
	f=MatrixAlloc(m,m);
	for(int i=0;i<m;i++)
	for(int j=0;j<m;j++)
		if(i==j) f[i][j]=array[i];
	row=m;
	col=m;
	return TRUE;
}

template<class Type>
CMatrix<Type> CMatrix<Type>::Construct(CMatrix<Type> &mat,CString style)
{
	int i,j;
	CMatrix<Type> tmp;
	if(style=="LR"||style=="lr")
	{
		if(row!=mat.row)    return tmp;
		if(!tmp.Malloc(row,col+mat.col)) return tmp;
		for(i=0;i<tmp.row;i++)
		{
			for(j=0;j<tmp.col;j++)
			{
				if(j<col) tmp.f[i][j]=f[i][j];
				else      tmp.f[i][j]=mat.f[i][j-col];
			}
		}
	}
	return tmp;
}

template<class Type>	
CMatrix<Type> CMatrix<Type>::GetChild(int sr,int sc,int er,int ec)
{
	int i,j;
	CMatrix<Type> tmp(er-sr+1,ec-sc+1);
	for(i=0;i<tmp.row;i++)
	for(j=0;j<tmp.col;j++)
		tmp[i][j]=f[sr+i][sc+j];
	return tmp;