matrix.cpp

非线形最小二乘程序,希望对大家有用哈,我现在急需APRIOR算法的C++程序,有的请给我发啊一份,

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

#include "Matrix.h"
#include "math.h"
#include "Vector.h"
#include <memory.h>

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

CMatrix::CMatrix(const CMatrix &M)
{
    m_nv=M.m_nv;
	m_ne=M.m_ne;
	m_browvector=M.m_browvector;
	m_val=new CVector[m_nv];
	for(int i=0;i<m_nv;i++)
		m_val[i]=M.m_val[i];
}
CMatrix::CMatrix(int nv ,int ne,bool browvector)
{
	m_nv=nv;
	m_ne=ne;
	m_val=new CVector[nv];
	for(int i=0;i<nv;i++)
		m_val[i].Resize(ne);
	m_browvector=browvector;
}

CMatrix &CMatrix::operator =(const CVector &M)
{
    delete []m_val;
	/*
	m_browvector=M.m_brow;
	if(M.m_brow)
	{
		m_nv=1;
		m_ne=M.m_p;
	}
	else
	{
		m_nv=M.m_p;
		m_ne=1;
	}
	m_val=new CVector[m_nv];
	if(M.m_brow)
	{
		m_val[0].m_p=M.m_p;
		for(int i=0;i<m_ne;i++)
			m_val[0][i]=M.m_val[i];
	}
	else
	{
		for(int i=0;i<m_nv;i++)
		{
			m_val[i].m_p=1;
			m_val[i][0]=M.m_val[i];
		}
	}*/
	m_browvector=1;
	m_nv=1;
	m_ne=M.m_p;
	m_val=new CVector[m_nv];
	m_val[0]=M;
	return *this;
}

CMatrix &CMatrix::operator =(const CMatrix &M)
{
    delete []m_val;
	m_nv=M.m_nv;
	m_ne=M.m_ne;
	m_browvector=M.m_browvector;
	m_val=new CVector[m_nv];
	for(int i=0;i<m_nv;i++)
		m_val[i]=M.m_val[i];
	return *this;

}

CMatrix CMatrix::operator *( const double &n)
{
	CMatrix C(m_nv,m_ne);
		for(int i=0;i<m_nv;i++)
			for(int j=0;j<m_ne;j++)
			{
				C[i][j]=0;
				for(int k=0;k<m_ne;k++)
					C[i][j]+=m_val[i][k]*n;

			}
		return C;
}

CMatrix CMatrix::operator *( const CMatrix  &M)
{
	CMatrix B=M;
   	if(m_ne==B.m_nv)
	{
		int n=m_nv;
		int p=B.m_ne;
		int m=m_ne;
		CMatrix C(n,p);
		for(int i=0;i<n;i++)
			for(int j=0;j<p;j++)
			{
				C[i][j]=0;
				for(int k=0;k<m;k++)
					C[i][j]+=m_val[i][k]*B[k][j];
			}
		return C;

	}
	else
		throw("no match size");

	return B;
}

CMatrix::~CMatrix()
{
		delete[] m_val;
}

CVector &CMatrix::operator [](int index)
{
	if(index<m_nv &&index>=0)
		return m_val[index];
	else
		throw("index out of range");
}
CMatrix CMatrix::Inverse()//(int n, float **a,float**m)
{
	if(m_nv!=m_ne)
		throw(" rows!=cols eror");
	else
	{
		int	n=m_nv;
		register i,j,k;
	//Set m as Unit Matrix
	/*	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
		  m[i][j]= (i==j) ? 1.0f : 0.0f;*/
		CMatrix m(n,n);
		m.Unit();
		double w,y;
		CVector c(n);
	//float *c=new float [n];

		for (i=1;i<=n;++i)
		{
			k=i;
			y=m_val[i-1][i-1];
			for (j=i+1;j<=n;++j)
			{
				w=m_val[j-1][i-1];
				if(fabs(w)>fabs(y))
				{
					k=j;
					y=w;
				}
			}
			if(fabs(y)<1e-8)
				throw("data abnormal ");		
      		y=1.0f/y;
			for(j=1;j<=n;++j) 
			{
				c[j-1]=m_val[k-1][j-1];
				m_val[k-1][j-1]=m_val[i-1][j-1];
				m_val[i-1][j-1]=c[j-1];
			
				c[j-1]=m[k-1][j-1];
				m[k-1][j-1]=m[i-1][j-1];
				m[i-1][j-1]=c[j-1];
			}

			for (j=1;j<=n;++j)  
			{
				m_val[i-1][j-1]=m_val[i-1][j-1]*y;
				m[i-1][j-1]=m[i-1][j-1]*y;
			}
			for (k=1;k<=n;++k)
				if (k!=i) 
				{
					y=m_val[k-1][i-1];
					for(j=n;j>=1;--j)
					{
						m[k-1][j-1]-=m[i-1][j-1]*y;
						m_val[k-1][j-1]-=m_val[i-1][j-1]*y;
					}
				}
	 
		}
		return m;  
	}
}

CMatrix CMatrix::Transpos()
{
	CMatrix m;
	m.Resize(m_ne,m_nv);
	for(int i=0;i<m_ne;i++)
		for(int j=0;j<m_nv;j++)
			m[i][j]=m_val[j][i];
	return m;
}

CMatrix CMatrix::Standard()
{
	CMatrix m;
	m.Resize(m_nv,m_ne);
	CVector s(m_ne);
	CVector aver(m_ne);
	for(int i=0;i<m_ne;i++)
	{
		aver[i]=0.0;
		s[i]=0.0;
	}
	for( i=0;i<m_nv;i++)
		for(int j=0;j<m_ne;j++)
			aver[j]+=m_val[i][j]/m_nv;
    for(i=0;i<m_nv;i++)
		for(int j=0;j<m_ne;j++)
			s[j]+=(m_val[i][j]-aver[j])*(m_val[i][j]-aver[j]);
	for(i=0;i<m_ne;i++)
	{
		s[i]=sqrt((s[i])/m_nv);
	}
    for( i=0;i<m_nv;i++)
		for(int j=0;j<m_ne;j++)
			m[i][j]=(m_val[i][j]-aver[j])/s[j];
    return m;
}

CMatrix CMatrix::StandardNonLinear()
{
	CMatrix M(*this) ;
	for ( int i = 0; i < M.m_nv; i++ ) 
	{
		for ( int j = 0; j < M.m_nv; j++ ) 
		{
			M[i][j] = M[i][j]/sqrt(M[i][i])/sqrt(M[j][j]) ;
		}
	}
	return M ;
}
	    	
void CMatrix::Unit()
{
	if(m_nv!=m_ne)
		throw("rows!=cols errro");
	else
	{
		int n=m_nv;
		for(register i=0;i<n;i++)
			for(register j=0;j<n;j++)
				m_val[i][j]= (i==j) ? 1.0 : 0.0;
	}

}

CMatrix CMatrix::CaclCov()
{
	register int i,j,k;
	CMatrix S(m_ne,m_ne);
	double *x=new double[m_ne];

    memset(x,0,m_ne*sizeof(float));
	for(i=0;i<m_ne;i++)
    {
		for(j=0;j<m_nv;j++)
			x[i]+=m_val[j][i];
		x[i]=x[i]/m_nv;
	}

	for(i=0;i<m_ne;i++)
    {
		for(j=0;j<m_ne;j++)
        {
			S[i][j]=0;
			for( k=0;k<m_nv;k++)
				S[i][j]+=(m_val[k][i]-x[i])*(m_val[k][j]-x[j]);
			S[i][j]=S[i][j]/(m_nv-1);
		}
	}
	delete[] x;
	return S;
}

int CMatrix::GetRows()
{
	return m_nv;
}

int CMatrix::GetCols()
{
	return m_ne;
}

void CMatrix::Resize(int nv,int ne)
{
	m_nv=nv;
	m_ne=ne;
	delete []m_val;
	m_val=new CVector[nv];
	for(int i=0;i<nv;i++)
		m_val[i].Resize(ne);

}

/*
CMatrix::CMatrix(CMatrix &M)
{
	m_nv=M.m_nv;
	m_ne=M.m_ne;
	m_browvector=M.m_browvector;
	m_val=new CVector[m_nv];
	for(int i=0;i<m_nv;i++)
		m_val[i]=M.m_val[i];

}
*/


CMatrix CMatrix::ProduceMatrix( CMatrix B)
{
	
	if(m_ne==B.m_nv)
	{
		int n=m_nv;
		int p=B.m_ne;
		int m=m_ne;
		CMatrix C(n,p);
		for(int i=0;i<n;i++)
			for(int j=0;j<p;j++)
			{
				C[i][j]=0;
				for(int k=0;k<m;k++)
					C[i][j]+=m_val[i][k]*B[k][j];

			}
		return C;
	}
	else
		throw("no match size");
}
ifstream & operator >>(ifstream &in,CMatrix &M)
{
	int n,e;
	in>>n;
	in>>e;
	M.Resize(n,e);
	for(int i=0;i<M.m_nv;i++)
		for(int j=0;j<M.m_ne;j++)
			in>>M.m_val[i][j];
	return in;
}
ofstream & operator <<(ofstream &out,CMatrix &M)
{
	out<<M.m_nv<<" ";
	out<<M.m_ne<<" ";
	for(int i=0;i<M.m_nv;i++)
		for(int j=0;j<M.m_ne;j++)
			out<<M.m_val[i][j]<<" ";
	return out;
}
ostream & operator <<(ostream &out,CMatrix &M)
{
	for(int i=0;i<M.m_nv;i++)
		for(int j=0;j<M.m_ne;j++)
			out<<M.m_val[i][j]<<" ";
	return out;
}
void CMatrix::DeleteRow(int r)
{
	m_nv--;
	for(int i=r;i<m_nv;i++)
		m_val[i]=m_val[i+1];
}

void CMatrix::DeleteCol(int c)
{
	m_ne--;
	for(int i=0;i<m_nv;i++)
	{
		for(int j=c;j<m_ne;j++)
			m_val[i][j]=m_val[i][j+1];
	}
}

CMatrix CMatrix::operator +(const CMatrix &B)
{
	CMatrix mat(m_nv,m_ne);
	CMatrix M=B;
	if(M.m_nv==m_nv&&M.m_ne==m_ne)
	{
		for(int i=0;i<m_nv;i++)
			for(int j=0;j<m_ne;j++)
				mat[i][j]=M[i][j]+m_val[i][j];
	    return mat;		
	}
	else
		throw("errow");
}
CMatrix CMatrix::operator -(const CMatrix &B)
{
	CMatrix mat(m_nv,m_ne);
	CMatrix M=B;
	if(M.m_nv==m_nv&&M.m_ne==m_ne)
	{
		for(int i=0;i<m_nv;i++)
			for(int j=0;j<m_ne;j++)
				mat[i][j]=M[i][j]-m_val[i][j];
        return mat;
	}
	else
		throw("errow");
}

CVector CMatrix::operator *(const CVector &V)
{
	CVector vec(m_nv,false);
	for(int i=0;i<m_nv;i++)
		vec[i]=0;
	CVector B=V;
	if(B.m_p==m_ne)
	{
        for(int i=0;i<m_nv;i++)
			for(int j=0;j<m_ne;j++)
				vec[i]+=m_val[i][j]*B[j];
	    return vec;
	}
	else 
		throw("errow");
}


void CMatrix::FreeSpace()
{
	delete[] m_val;
}

CVector CMatrix::ColMax()//增广矩阵M
{
	int iRow = m_nv ;
	CVector X( iRow, false ) ;//解向量
	CMatrix C(iRow,iRow);//消去时所乘的系数
	CVector* pr=NULL;//存放r行首地址
	CVector* pk=NULL;//存放k行首地址

	for(int k=0;k<iRow;k++)
	{
		if(m_val[k][k]==0)
		{
			cout<<"主对角线元素为0"<<endl;
		}
		//选取矩阵Ａ最大值
		int r=k;//记录列主元的行数
		double temp=0;//存放最大值
		for(int i1=k;i1<iRow;i1++)
		{
			if(fabs(m_val[i1][k])>temp)
			{
				temp=fabs(m_val[i1][k]);
				r=i1;
			}
		}

		//如果当前行主对角线元素不是最大，则交换增广矩阵Ａ的r行k行
		if(r!=k)
		{
			pr=&m_val[r];
			pk=&m_val[k];
			exchange(m_val[r],m_val[k],iRow+1);
		}
		for(int i=k+1;i<iRow;i++)
		{
			C[i][k]=m_val[i][k]/m_val[k][k];
			m_val[i][iRow]=m_val[i][iRow]-C[i][k]*m_val[k][iRow];
			for(int j=k+1;j<iRow;j++)
			{
				m_val[i][j]=m_val[i][j]-C[i][k]*m_val[k][j];//化为三角型矩阵
			}
		}
	}
//	cout<<A[0][0]<<endl<<A[0][1]<<endl<<A[0][2]<<endl<<A[1][0]<<endl<<A[1][1]<<endl<<A[1][2]<<endl;
	
	//回代求解
	X[iRow-1]=m_val[iRow-1][iRow]/m_val[iRow-1][iRow-1];//最后一个x值
	//cout<<"X"<<iRow<<"="<<X[iRow-1]<<endl;
	for(int l=iRow-2;l>=0;l--)//回代求解
	{
		double c=0;
		for(int p=l+1;p<iRow;p++)
		{
			c=c+m_val[l][p]*X[p];
		}
		X[l]=(m_val[l][iRow]-c)/m_val[l][l];
		//cout<<"X"<<l+1<<"="<<X[l]<<endl;
	}
	return X;
}

void CMatrix::exchange(CVector& pr,CVector& pk,int col)         //交换增广矩阵Ａ的r行k行的函数，pr：r行首地址，pk：k行首地址
{
	double temp;
	for(int i=0;i<col;i++)
	{
		temp=pr[i];
		pr[i]=pk[i];
		pk[i]=temp;
	}
}