matrix.cpp

原创C++的matrix类,有许多运算符重载,很简单,但是好用,好看

// Matrix.cpp: implementation of the CMatrix class.
//
//////////////////////////////////////////////////////////////////////

#include "Matrix.h"
#include <cmath>

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////


CMatrix::CMatrix()
{
	m_nr=0;
	m_nc=0;
	m_data=0;
}

CMatrix::~CMatrix()
{
	clear();
}

void CMatrix::create(int nr, int nc)
{
	m_nr=nr;
	m_nc=nc;
	m_data=new double*[m_nr];
	for(int i=0;i<m_nr;i++)
	{
		m_data[i]=new double[m_nc];
		memset(m_data[i],0,sizeof(double)*m_nc );
	}
}

void CMatrix::clear()
{
	if(m_data!=0)
	{
		for(int i=0;i<m_nr;i++)
		{
			delete[] m_data[i];
		}
		delete[] m_data;
		m_data=0;
		m_nr=0;
		m_nc=0;
	}
}

CMatrix::CMatrix(int nr, int nc)
{
	create(nr, nc);
}

CMatrix::CMatrix(double **data, int nr, int nc)
{
	create(nr,nc);
	for(int i=0;i<m_nr;i++)
	{
		memcpy(m_data[i],data[i],sizeof(double)*m_nc );
	}
}

CMatrix::CMatrix(const CMatrix &m)
{
	create(m.m_nr,m.m_nc);
	for(int i=0;i<m_nr;i++)
	{
		memcpy(m_data[i],m.m_data[i],sizeof(double)*m_nc );
	}
}

CMatrix CMatrix::operator =(const CMatrix &m)
{
	if(&m==this) return *this;
	clear();
	create(m.m_nr,m.m_nc);
	for(int i=0;i<m_nr;i++)
	{
		memcpy(m_data[i],m.m_data[i],m_nc*sizeof(double));
	}
	return *this;
}

CMatrix CMatrix::operator *(const CMatrix &m)
{
	CMatrix mp(m_nr,m.m_nc);
	for(int i=0;i<m_nr;i++)
	{
		for(int j=0;j<m.m_nc;j++)
		{
			mp[i][j]=0;
			for(int k=0;k<m_nc;k++)
			{
				mp[i][j]+=m_data[i][k]*m[k][j];
			}
		}
	}
	return mp;
}

CMatrix CMatrix::operator +(const CMatrix &m)
{
	CMatrix mp(*this);
	for(int i=0;i<m_nr;i++)
	{
		for(int j=0;j<m_nc;j++) 
		{
			mp[i][j]+=m[i][j];
		}
	}
	return mp;
}

CMatrix CMatrix::operator -(const CMatrix &m)
{
	CMatrix mp(*this);
	for(int i=0;i<m_nr;i++)
	{
		for(int j=0;j<m_nc;j++) 
		{
			mp[i][j]-=m[i][j];
		}
	}
	return mp;
}

CMatrix CMatrix::operator *(double x)
{
	CMatrix m(*this);
	for(int i=0;i<m_nr;i++)
	{
		for(int j=0;j<m_nc;j++)
		{
			m[i][j]=m_data[i][j]*x;
		}
	}
	return m;
}

CMatrix CMatrix::operator /(double x)
{
	CMatrix m(*this);
	for(int i=0;i<m_nr;i++)
	{
		for(int j=0;j<m_nc;j++)
		{
			m[i][j]=m_data[i][j]/x;
		}
	}
	return m;
}

CMatrix CMatrix::operator -()
{	
	CMatrix m(*this);
	for(int i=0;i<m_nr;i++) 
	{
		for(int j=0;j<m_nc;j++) 
		{
			m[i][j]=-m[i][j];
		}
	}
	return m;
}

CMatrix CMatrix::operator ~()
{
	CMatrix m(m_nc,m_nr);
	for(int i=0;i<m_nc;i++)
	{
		for(int j=0;j<m_nr;j++)
		{
			m[i][j]=m_data[j][i];
		}
	}
	return m;
}

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

CMatrix CMatrix::operator -=(const CMatrix &m)
{
	*this=*this-m;
	return *this;
}

CMatrix CMatrix::operator *=(const CMatrix &m)
{
	*this=*this*m;
	return *this;
}

CMatrix CMatrix::operator *=(double x)
{
	*this=*this*x;
	return *this;
}

CMatrix CMatrix::operator /=(double x)
{
	*this=*this/x;
	return *this;
}

bool CMatrix::operator ==(const CMatrix &m)
{
	double conv=0;
	double delta=0;
	for(int i=0;i<m_nr;i++)
	{
		for(int j=0;j<m_nc;j++)
		{
			if( fabs(m_data[i][j]-m[i][j])>1.0e-13) return false;
		}
	}
	return true;
}

CMatrix CMatrix::operator !()
{
	CMatrix AI(step(*this|identity(m_nr)));
	CMatrix invs(m_nr, m_nc);
	for(int i=0;i<m_nr;i++)
	{
		for(int j=0;j<m_nc;j++)
		{
			invs[i][j]=AI[i][j+m_nc];
		}
	}
	return invs;
}

CMatrix CMatrix::operator |(const CMatrix &m)
{
	CMatrix mp(m_nr, m_nc+m.m_nc);
	for(int i=0;i<m_nr;i++)
	{
		for(int j=0;j<m_nc;j++)
		{
			mp[i][j]=m_data[i][j];
		}
		for(j=0;j<m.m_nc;j++)
		{
			mp[i][j+m_nc]=m[i][j];
		}
	}
	return mp;
}



CMatrix operator *(double x, const CMatrix &m)
{	
	CMatrix mp(m);
	for(int i=0;i<m.m_nr;i++)
	{
		for(int j=0;j<m.m_nc;j++)
		{
			mp[i][j]=x*m[i][j];
		}
	}
	return mp;
}


istream & operator >>(istream & input, CMatrix &m)
{
	int nr=0;
	int nc=0;
	input>>nr>>nc;
	m=CMatrix(nr,nc);
	for(int i=0;i<nr;i++)
	{
		for(int j=0;j<nc;j++)
		{
			input>>m[i][j];
		}
	}
	return input;
}

ostream & operator <<(ostream & output, CMatrix &m)
{
	int nr=m.getnr();
	int nc=m.getnc();
	output<<nr<<"  "<<nc<<endl;
	for(int i=0;i<nr;i++)
	{
		for(int j=0;j<nc;j++)
		{
			output<<m[i][j]<<"  ";
		}
		output<<endl;
	}
	output<<endl;
	return output;
}


void householder(CMatrix &m, CMatrix &d, CMatrix &e)
{
	int l,k,j,i;
	double scale, hh, h, g, f;

	int n=m.m_nc;

	for(i=n-1;i>0;i--)
	{
		l=i-1;
		h=scale=0;
		if(l>0)
		{
			for(k=0;k<l+1;k++)
			{
				scale+=fabs(m[i][k]);
			}
			if(scale==0.0)
			{
				e[i][0]=m[i][l];
			}
			else
			{
				for(k=0;k<l+1;k++)
				{
					m[i][k]/=scale;					
					h+=m[i][k]*m[i][k];					
				}
				f=m[i][l];
				g=(f>=0.0?-sqrt(h):sqrt(h));
				e[i][0]=scale*g;
				h-=f*g;
				m[i][l]=f-g;
				f=0.0;

				for(j=0;j<l+1;j++)
				{
					m[j][i]=m[i][j]/h;
					g=0.0;
					for(k=0;k<j+1;k++)
					{
						g+=m[j][k]*m[i][k];
					}
					for(k=j+1;k<l+1;k++)
					{
						g+=m[k][j]*m[i][k];
					}
					e[j][0]=g/h;
					f+=e[j][0]*m[i][j];
				}
				hh=f/(h+h);
				for(j=0;j<l+1;j++)
				{
					f=m[i][j];
					e[j][0]=g=e[j][0]-hh*f;
					for(k=0;k<j+1;k++)
					{
						m[j][k]-=( f*e[k][0]+g*m[i][k] );
					}
				}
			}
		}
		else
		{
			e[i][0]=m[i][l];
		}
		d[i][0]=h;
	}

	d[0][0]=0.0;
	e[0][0]=0.0;

	for(i=0;i<n;i++)
	{
		l=i;
		if(d[i][0]!=0.0)
		{
			for(j=0;j<l;j++)
			{
				g=0.0;
				for(k=0;k<l;k++)
				{
					g+=m[i][k]*m[k][j];
				}
				for(k=0;k<l;k++)
				{
					m[k][j]-=g*m[k][i];
				}
			}
		}
		d[i][0]=m[i][i];
		m[i][i]=1.0;
		for(j=0;j<l;j++)
		{
			m[j][i]=m[i][j]=0.0;
		}
	}
}

void ql(CMatrix &z, CMatrix &d, CMatrix &e)
{
	int m,l,iter,i,k;
	double s,r,p,g,f,dd,c,b;

	int n=z.m_nc;
	for(i=1;i<n;i++)	
		e[i-1][0]=e[i][0];	
	e[n-1][0]=0.0;

	for(l=0;l<n;l++)
	{
		iter=0;
		do
		{
			for(m=l;m<n-1;m++)
			{
				dd=fabs(d[m][0])+fabs(d[m+1][0]);
				if(fabs(e[m][0])+dd==dd)  break;
			}
			if(m!=l)
			{
				if(iter++==30) 	return;
				g=(d[l+1][0]-d[l][0])/(e[l][0]+e[l][0]);
				r=sqrt(g*g+1.0);
				g=d[m][0]-d[l][0]+e[l][0]/( g + g>0?(r>0?r:-r):(r>0?-r:r) );
				s=c=1.0;
				p=0.0;
				for(i=m-1;i>=l;i--)
				{
					f=s*e[i][0];
					b=c*e[i][0];
					e[i+1][0]=(r=sqrt(f*f+g*g));
					if(r==0.0)
					{
						d[i+1][0]-=p;
						e[m][0]=0.0;
						break;
					}
					s=f/r;
					c=g/r;
					g=d[i+1][0]-p;
					r=(d[i][0]-g)*s+2.0*c*b;
					d[i+1][0]=g+(p=s*r);
					g=c*r-b;
					for(k=0;k<n;k++)
					{
						f=z[k][i+1];
						z[k][i+1]=s*z[k][i]+c*f;
						z[k][i]=c*z[k][i]-s*f;
					}
				}
				if(r==0.0 && i>=l)  continue;
				d[l][0]-=p;  e[l][0]=g;  e[m][0]=0.0;
			}
		}while(m!=l);
	}
}

CMatrix symmeig(const CMatrix &H, CMatrix &E)
{
	CMatrix m(H);
	CMatrix d(m.m_nc,1);
	CMatrix e(m.m_nc,1);
	householder(m, d, e);
	ql(m, d, e);
	E=CMatrix(m.m_nc,m.m_nc);
	for(int i=0;i<m.m_nc;i++) E[i][i]=d[i][0];
	return m;
}


CMatrix linequ(const CMatrix & A, const CMatrix & B)
{
	CMatrix a(A);
	CMatrix b(B);
	int n=a.m_nr;
	CMatrix x(n,1);
	for(int i=0;i<n-1;i++)// ok
	{
		int pos=i;
		double temp=a[i][i];
		for(int j=i+1;j<n;j++)
		{
			if(fabs(temp)<fabs(a[j][i]))
			{
				temp=a[j][i];
				pos=j;
			}
		}
		if(i!=pos)
		{
			for(int k=i;k<n;k++)
			{
				swap(a[i][k],a[pos][k]);
			}
			swap(b[i][0],b[pos][0]);
		}

		for(j=i+1;j<n;j++)
		{
			double temp=a[j][i]/a[i][i];
			for(int k=i+1;k<n;k++)
			{
				a[j][k]-=temp*a[i][k];
			}
			b[j][0]-=temp*b[i][0];
		}
	}

	for(i=n-1;i>=0;i--)
	{
		x[i][0]=b[i][0];
		for(int j=i+1;j<n;j++)
		{
			x[i][0]-=a[i][j]*x[j][0];
		}
		x[i][0]/=a[i][i];
	}
	return x;
}


double determinant(const CMatrix &m)
{
	CMatrix a(m);
	int n=m.m_nr;
	double det=1;

	for(int i=0;i<n-1;i++)// ok
	{
		int pos=i;
		double temp=a[i][i];
		for(int j=i+1;j<n;j++)
		{
			if(fabs(temp)<fabs(a[j][i]))
			{
				temp=a[j][i];
				pos=j;
			}
		}

		if(fabs(temp)<1e-18) return 0;

		if(i!=pos)
		{
			for(int k=i;k<n;k++)
			{
				swap(a[i][k],a[pos][k]);
			}
			det=-det;
		}
		det*=a[i][i];

		for(j=i+1;j<n;j++)
		{
			double temp=a[j][i]/a[i][i];
			for(int k=i;k<n;k++)
			{
				a[j][k]-=temp*a[i][k];
			}
		}
	}

	return det*a[n-1][n-1];
}

CMatrix identity(int n)
{
	CMatrix m(n,n);
	for(int i=0;i<n;i++)
	{
		m[i][i]=1;
	}
	return m;
}

CMatrix step(const CMatrix &m)
{
	CMatrix m1(m);
	int i=0;
	int j=0;
	while(i<m1.m_nr && j<m1.m_nc)
	{
		double max=m1[i][j];
		int pos=i;
		for(int k=i+1;k<m1.m_nr;k++)
		{
			if( fabs(max) < fabs(m1[k][j]) )
			{
				max=m1[k][j];
				pos=k;
			}
		}

		if(fabs(max)>1e-015)
		{
			for(int s=j;s<m1.m_nc;s++)
			{
				swap(m1[i][s],m1[pos][s]);
			}

			for(int t=0;t<m1.m_nr;t++)
			{
				double scale=m1[t][j]/m1[i][j];
				if(t!=i)
				{
					for(int s=j;s<m1.m_nc;s++)
					{
						m1[t][s]-=m1[i][s]*scale;
					} 
					
				}				
			}

			for(s=j+1;s<m1.m_nc;s++)
			{
				m1[i][s]/=m1[i][j];
			}

			m1[i][j]=1;
			
			i++;
		}

		j++;
	}
	return m1;
}

double trace(const CMatrix &m)
{
	double tr=0;
	for(int i=0;i<m.getnr();i++)
	{
		tr+=m[i][i];
	}
	return tr;
}

CMatrix rand(int nr, int nc)
{
	CMatrix m(nr, nc);
	for(int i=0;i<m.m_nr;i++)
	{
		for(int j=0;j<m.m_nc;j++)
		{
			m[i][j]=rand()%10;
		}
	}
	return m;
}

CMatrix linequ2(const CMatrix & A, const CMatrix & B)
{
	CMatrix a(A);
	CMatrix b(B);
	CMatrix AB(step(a|b));
	CMatrix x(a.m_nr,1);
	for(int i=0;i<a.m_nr;i++)
	{
		x[i][0]=AB[i][a.m_nr];
	}
	return x;
}