matrix.cpp

坐标转换程序,可以将空间直角坐标转化为大地坐标

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

#include "stdafx.h"
#include "Matrix.h"
#include <math.h>
#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

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

Matrix::Matrix()
{
  M=N=0;
  Data=NULL;
}

Matrix::~Matrix()
{
 if ((M!=0)||(N!=0))
 if (Data!=NULL) delete []Data ;
 Data=NULL;
}

Matrix::Matrix(UINT m)
{
	Success=FALSE;
	if ((m<1)) { Data=NULL;return;}
    Data=new double [m*1]; 
	M=m;N=1;
	if (Data!=NULL) Success=TRUE;
}


Matrix::Matrix(UINT m, UINT n)
{
	Success=FALSE;
	if ((m<1)||(n<1)) { Data=NULL;return;}
    Data=new double [m*n]; 
	M=m;N=n;
	if (Data!=NULL) Success=TRUE;
}

BOOL Matrix::SetData(UINT m/*, UINT n=0 */, double a)
{
    if (m>=M) return FALSE;
	Data[m*N+0]=a;
	return TRUE;
}

BOOL Matrix::SetData(UINT m, UINT n, double a)
{
    if ((m>=M)||(n>=N)) return FALSE;
	Data[m*N+n]=a;
	return TRUE;
}

Matrix::Matrix(Matrix &mat)
{
	M=mat.GetRow(); N=mat.GetCol();
	if((M<=0) || (N<=0)){
          Success=FALSE;return;
	}
	Data=new double[M*N];
	if(!Data){
          Success=FALSE;return;
	}
	for (i=0;i<M;i++)
		for (j=0;j<N;j++)
			SetData(i,j,mat.GetData(i,j));
    Success=TRUE;
}

UINT Matrix::GetCol()
{
 return N;
}

UINT Matrix::GetRow()
{
 return M;
}

double Matrix::GetData(UINT m, UINT n)
{
  return Data[m*N+n];
}

double Matrix::GetData(UINT m/*, UINT n=0*/)
{
  return Data[m*N];
}


BOOL Matrix::Alloc(UINT m)
{
	return Alloc(m, 1);
}

BOOL Matrix::Alloc(UINT m, UINT n)
{
	if ((m<1)||(n<1))  return FALSE;
    delete []Data;Success=FALSE;
    Data=new double [m*n]; 
	M=m;N=n;
	if (Data!=NULL) Success=TRUE;
    return TRUE;
}

void Matrix::operator =(double a)
{
	for (i=0;i<M;i++)
		for (j=0;j<N;j++)
			SetData(i,j,a);
}

Matrix Matrix::operator + (Matrix &mat)
{
	if(mat.GetRow()!=M || mat.GetCol()!=N){
		Success=FALSE;
		return mat;
	}

	Matrix m_temp(M,N);
	for (i=0;i<M;i++)
		for (j=0;j<N;j++)
		 m_temp.SetData(i,j, GetData(i,j)+mat.GetData(i,j));

	Success=TRUE;
	return m_temp;
}

Matrix Matrix::operator + (double a)
{
	Matrix m_temp(M,N);
	for (i=0;i<M;i++)
		for (j=0;j<N;j++)
		 m_temp.SetData(i,j, GetData(i,j)+a);
	Success=TRUE;
	return m_temp;
}


Matrix Matrix::operator - (Matrix &mat)
{
	if(mat.GetRow()!=M || mat.GetCol()!=N){
		Success=FALSE;
		return mat;
	}

	Matrix m_temp(M,N);
	for (i=0;i<M;i++)
		for (j=0;j<N;j++)
		 m_temp.SetData(i,j, GetData(i,j)-mat.GetData(i,j));

	Success=TRUE;
	return m_temp;
}

Matrix Matrix::operator - (double a)
{
	Matrix m_temp(M,N);
	for (i=0;i<M;i++)
		for (j=0;j<N;j++)
		 m_temp.SetData(i,j, GetData(i,j)-a);

	Success=TRUE;
	return m_temp;
}

Matrix Matrix::operator * (Matrix &mat)
{
	UINT m,n,k;
    m=mat.GetRow();
    n=mat.GetCol();
	Success=FALSE;
    if (m!=N) return mat;
	Matrix m_temp(M,n);
    double t;

	for(i=0;i<M;i++)
		for(j=0;j<n;j++){
			t=0.0;
			for(k=0;k<N;k++)
				t += GetData(i,k) * mat.GetData(k,j);
			m_temp.SetData(i,j,t);
		}

	return m_temp;
}


Matrix Matrix::operator * (double a)
{
	Matrix m_temp(M,N);
	for (i=0;i<M;i++)
		for (j=0;j<N;j++)
		 m_temp.SetData(i,j, GetData(i,j)*a);

	Success=TRUE;
	return m_temp;
}


BOOL Matrix::operator =(Matrix &mat)
{
	if(((mat.GetRow() != M) || (mat.GetCol() != N)) ){
		delete []Data;
		Data = NULL;
	}
	if(Data == NULL ){
		M = mat.GetRow();
		N = mat.GetCol();
		Data = new double[M*N];
		if(Data == NULL) return FALSE;
	}

	for (i=0;i<M;i++)
		for (j=0;j<N;j++)
			SetData(i,j,mat.GetData(i,j));
   return TRUE;
}
BOOL Matrix::operator +=(Matrix &mat)
{
	if(((mat.GetRow() != M) || (mat.GetCol() != N)) ){
       return FALSE;
	}
	for (i=0;i<M;i++)
		for (j=0;j<N;j++)
			SetData(i,j,GetData(i,j)+mat.GetData(i,j));
   return TRUE;
}

BOOL Matrix::operator +=(double a)
{
	for (i=0;i<M;i++)
		for (j=0;j<N;j++)
			SetData(i,j,GetData(i,j)+a);
   return TRUE;
}

BOOL Matrix::operator -=(Matrix &mat)
{
	if(((mat.GetRow() != M) || (mat.GetCol() != N))){
       return FALSE;
	}
	for (i=0;i<M;i++)
		for (j=0;j<N;j++)
			SetData(i,j,GetData(i,j)-mat.GetData(i,j));
   return TRUE;
}

BOOL Matrix::operator -=(double a)
{
	for (i=0;i<M;i++)
		for (j=0;j<N;j++)
			SetData(i,j,GetData(i,j)-a);
   return TRUE;
}

Matrix Matrix::MatInv()
{
	UINT k;
	double s;
	Matrix m_temp(M,N);
    Success=FALSE;
	if (M!=N) return m_temp;
    
	for (i=0;i<M;i++)
		for (j=0;j<N;j++)
			m_temp.SetData(i,j,GetData(i,j));

    for (i=0;i<M;i++)
		for (j=i;j<M;j++)
        {
         if (i!=0)
			 for (k=0;k<=i-1;k++) 
			 {
               m_temp.SetData(i,j,m_temp.GetData(i,j)-m_temp.GetData(k,i)*m_temp.GetData(k,j));
			 }

        if (j==i) m_temp.SetData(i,j,sqrt(m_temp.GetData(i,j)));
		   else   m_temp.SetData(i,j,m_temp.GetData(i,j)/m_temp.GetData(i,i));
		}
//15:
    for (i=0;i<M;i++)
	{
      m_temp.SetData(i,i,1/m_temp.GetData(i,i));
	  if (i==M-1) break;
	  for (j=i+1;j<M;j++)
	  {
        s=0;
        for (k=i;k<=j-1;k++) s=s-m_temp.GetData(i,k)*m_temp.GetData(k,j);
		m_temp.SetData(i,j,s/m_temp.GetData(j,j));
	  }
	}


    for (i=0;i<M;i++)
		for (j=i;j<M;j++)
		{
         s=0;
		 for (k=j;k<M;k++) s=s+m_temp.GetData(i,k)*m_temp.GetData(j,k);
         m_temp.SetData(i,j,s);
		}

    for (i=0;i<M;i++)
		for (j=i;j<M;j++)
           m_temp.SetData(j,i,m_temp.GetData(i,j));

    Success=TRUE;
	return m_temp;
}


Matrix Matrix::MatInv_All()
{
	int    *pivlst;
	char   *pivchk;
	UINT   i,j,k,l,lerow,lecol,l1;
	double piv,t,leval;
	double *det=new double;
	Matrix m_temp(M,N);

    Success=FALSE;
	if (M!=N) return m_temp;

	pivlst=(int *)calloc(100*2,2);
	pivchk=(char *)calloc(100,1);

    lecol=1;
//	pivlst=(int *)calloc(m*2,2);
//	pivchk=(char *)calloc(m,1);
	(*det) = 1.0;

	//copy the matrix to m_temp.
	for(i=0;i<=M-1;i++)
	{
		pivchk[i]=0;
		for(j=0;j<=N-1;j++)
			m_temp.SetData(i,j,GetData(i,j));
	}

	for(i=0;i<=M-1;++i){
		leval=0.0;
		for(j=0;j<=N-1;++j){
			if( !(pivchk[j]) ){
				for(k=0;k<=M-1;++k){
					if( !(pivchk[k]) ){
						if(fabs( m_temp.GetData(j,k) )>leval){
							lerow=j;
							lecol=k;
							leval=fabs(m_temp.GetData(j,k));
						}
					}
				}
			}
		}
		pivchk[lecol]=1;
		pivlst[i*2]=lerow;
		pivlst[i*2+1]=lecol;
		if(lerow != lecol){
			(*det)=-(*det);
			for(l=0;l<=M-1;++l){
				double x;
				x=m_temp.GetData(lerow,l);
                m_temp.SetData(lerow,l,m_temp.GetData(lecol,l));
                m_temp.SetData(lecol,l,x);
			}
		}
        
		piv=m_temp.GetData(lecol,lecol);
		(*det)=(*det)*piv;
//        if (piv==0) { piv=.1;AfxMessageBox("aaa"); };
		m_temp.SetData(lecol,lecol,1);
		for(l=0;l<=M-1;++l){
			m_temp.SetData(lecol,l,m_temp.GetData(lecol,l)/piv);
		}
		for(l1=0;l1<=M-1;++l1){
			if(l1!=lecol){
				t=m_temp.GetData(l1,lecol);
				m_temp.SetData(l1,lecol,0.0);
				for(l=0;l<=M-1;++l){
					m_temp.SetData(l1,l,m_temp.GetData(l1,l)
									  - m_temp.GetData(lecol,l)*t);
				}
			}
		}
	}
	for(i=0;i<=M-1;++i){
		l=M-i-1;
		if(pivlst[l*2] != pivlst[l*2+1]){
			lerow=pivlst[l*2];
			lecol=pivlst[l*2+1];
			for(k=0;k<=M-1;++k){
		       double x;
				x=m_temp.GetData(k,lerow);
                m_temp.SetData(k,lerow,m_temp.GetData(k,lecol));
                m_temp.SetData(k,lecol,x);
		
			}
		}
	}
	free(pivlst); free(pivchk);
	delete det;

	return m_temp; 
}
Matrix Matrix::MatTran()
{
	Matrix m_temp(N,M);
	for (i=0;i<M;i++)
		for (j=0;j<N;j++)
			m_temp.SetData(j,i,GetData(i,j));
	return m_temp;
}

BOOL Matrix::DuiCheng()
{
    if (M!=N) return FALSE;
    for (i=0;i<M;i++)
		for (j=i;j<M;j++)
           SetData(j,i,GetData(i,j));
    return TRUE;
}

double Matrix::Trace()
{
  double a;
  Success=FALSE;
  a=0;
  if (M!=N) return a;
  for (i=0;i<M;i++) a=a+GetData(i,i);
  Success=TRUE;
  return a;
}


Matrix Matrix::matChole(int *piErr)
{
	double sum;
    int    k;
	*piErr = TRUE;
	if(M != N) {
		printf("\n Matrix can't be converted!");
		exit(-4);
	}
	Matrix m_temp(M,N);

	m_temp=0;

	for(i=0;i<M;i++){
		for(j=i;j<N;j++){
			sum=GetData(i,j);
			for(k=(int)i-1;k>=0;k--)
				sum -= m_temp.GetData(i,k) * m_temp.GetData(j,k);
			if(i==j){
				if(sum<=0.0){
					*piErr = FALSE;
					m_temp=.0;
					for(k=0;k<(int)M;k++) 
						m_temp.SetData(k,k,1.0);
					return m_temp;//printf("\n Solve failed!\n"); exit(-4);
				}
				m_temp.SetData(i,i,sqrt(sum));
			}
			else m_temp.SetData(j,i, sum / m_temp.GetData(i,i));
		}
	}
	return m_temp;
}


#ifdef _DEBUG
BOOL Matrix::Print(CString FileName, int Znum,int Xnum)
{
 CString str;
 str.Format("%c%d.%dlf",'%',Znum,Xnum);
 FILE *out;
 out=fopen(FileName,"w+");
 if (out==NULL) return FALSE;
 for (i=0;i<M;i++)
 {
	 for (j=0;j<N;j++) fprintf(out,str,GetData(i,j));
     fprintf(out,"\n");
 }

 fclose(out);
 return TRUE;
}


BOOL Matrix::Print(FILE *out,CString FileName, int Znum,int Xnum)
{
 fprintf(out,"%s\n",FileName);
 CString str;
 str.Format("%c%d.%dlf",'%',Znum,Xnum);
 for (i=0;i<M;i++)
 {
	 for (j=0;j<N;j++) fprintf(out,str,GetData(i,j));
     fprintf(out,"\n");
 }
 return TRUE;
}

#endif //_DEBUG