maxadjclass.cpp

普通平差程序

#include <stdlib.h>
#include <math.h>
//***********矩阵类定义****************************************
//*************************************************************
//*************************************************************
class MAT
{        
	
    public:
        MAT();	
		MAT(int hang,int lie);
 MAT convert();
		~MAT();
		int row;
		int rank; 

		double *elem;	
       	void print();//print out the Max
		void input();
        void input1();
	private:
		char name[20]; 
};
void MAT::input()
{
	for(int i=0;i<row*rank;i++)
		cin>>elem[i];

}
void MAT::input1()
{
	if(row!=rank) for(int i=0;i<row*rank;i++)
		this->elem[i]=1;
     else for( int i=0;i<row;i++)
		 for(int j=0;j<rank;j++)
		 {  if(i==j)
			 this->elem[i+j+i*(rank-1)]=1;
		    else  this->elem[i+j+i*(rank-1)]=0;
		 }
}
MAT::MAT()
{
	
}

void MAT::print()
{
		for(int i=0;i<row;i++)
	{
		for(int j=0;j<rank;j++)
		cout<<" "<<*(elem+i+j+i*(rank-1));
		cout<<endl;}
}


MAT::MAT(int hang,int lie)
{//1.input the number of row and rank
	row=hang;
   	rank=lie;
//2.input the elements of the Mat
	elem=new double[hang*lie];

	for(int i=0;i<row*rank;i++)
		elem[i]=0;
//			
		
}
MAT::~MAT()
{
	;
}

MAT MAT::convert()
{
	MAT B;
	B.row=rank;
	B.rank=row;
	B.elem=new double[B.row*B.rank];
	  for(int i=0;i<B.rank;i++)
     for(int j=0;j<B.row;j++)
		{
		 int k=i+j+j*(B.rank-1);
		 B.elem[k]=*(elem+i+j+i*(B.row-1));
	 
	 }
 return B; 
}
//**************矩阵相乘函数定义*******************************
//*************************************************************
//*************************************************************
int tim(MAT A,MAT B,MAT C)
{   int i;  	
	if(A.rank!=B.row) return 0;
    
		 for(i=0;i<(C.row*C.rank);i++) C.elem[i]=0;
		for(i=0;i<C.row;i++)
		for(int j=0;j<C.rank;j++)
				for(int k=0;k<B.row;k++)
					C.elem[i+j+i*(C.rank-1)]+=A.elem[i+k+i*(A.rank-1)]*B.elem[k+j+k*(B.rank-1)];
		return 1;

}
//*********矩阵求逆函数定义：**********************************
//*************************************************************
//*************************************************************
	//函数说明：
	//1 MAX定义矩阵最大维数；
	//2 形参C[][]为输入矩阵，B[][]为输出矩阵
	//3 矩阵可逆,返回值为1/否则,返回值为-1；
int inverse(MAT C,MAT B)
{
	MAT A(C.row,C.row);
double b,js(0);
int m, i,j,k;
if(C.row!=C.rank) return 0;
m=C.row;
for(i=0;i<m;i++)
 for(j=0;j<m;j++)
 {
	 	A.elem[i+j+i*(A.rank-1)]=C.elem[i+j+i*(C.rank-1)]; 
     if(i==j)	B.elem[i+j+i*(B.rank-1)]=1;
	 else 	B.elem[i+j+i*(B.rank-1)]=0;
 }
	for(i=0;i<m;i++)
      { 
		  if(A.elem[i+i+i*(A.rank-1)]==0)
		    for(int l=i+1;l<m;l++)
			{if(A.elem[l+i+l*(A.rank-1)]!=0)
							{
								for(int k1=0;k1<m;k1++)
                                {A.elem[i+k1+i*(A.rank-1)]+=A.elem[l+k1+k1*(A.rank-1)];B.elem[i+k1+i*(B.rank-1)]+=B.elem[l+k1+l*(B.rank-1)];}
								goto BB;
			}}
			  	  
								  
BB:	    for(j=0;j<m;j++)
		{
			b=A.elem[j+i+j*(A.rank-1)]/A.elem[i+i+i*(A.rank-1)];
	    for(k=0;k<m;k++) 
			if(i!=j) 
			{ 
				A.elem[j+k+j*(A.rank-1)]-=b*A.elem[i+k+i*(A.rank-1)];
				B.elem[j+k+j*(B.rank-1)]-=b*B.elem[i+k+i*(B.rank-1)];
			}
		}
	  }
	  for(i=0;i<m;i++)
	  {
		  b=A.elem[i+i+i*(A.rank-1)]; if(fabs(b)<0.000000001) goto ss;
      for(j=0;j<m;j++)
	  {
            A.elem[i+j+i*(A.rank-1)]/=b; B.elem[i+j+i*(B.rank-1)]/=b;	  
	  }
	  } return 1; 
	  goto tt;
	  
ss:return 0; 
tt:;
}
//*****************矩阵求和函数定义****************************
//*************************************************************
//*************************************************************
int add(MAT A,MAT B,MAT AC)
{

	int i;
	if(A.row!=B.row || A.rank!=B.rank) { cout<<" can't do sub "<<endl;return 0;}
	for( i=0;i<A.row*A.rank;i++)
		AC.elem[i]=A.elem[i]+B.elem[i];

 return 1;
}

//*********************极大权平差类定义************************
//*************************************************************
//*************************************************************
class comadj{
public:
	MAT A;
	MAT P;
	MAT N_1;
	MAT X;
	MAT l;
	MAT v;
	comadj();
	~comadj();
	void datainput();
	void dataoutput();

private:
	int m;	
	int n;
	int k;
	double datadeal();

};
comadj::comadj()
{
	cout<<"输入观测值个数:"<<endl;cin>>m;
	cout<<"输 入 总 点 数:"<<endl;cin>>n;
    cout<<" 输入已知点数 :"<<endl;cin>>k;
    A.row=m;   A.rank=n;   A.elem=new double[m*n];
	P.row=m;   P.rank=m;   P.elem=new double[m*m];
	N_1.row=n; N_1.rank=n; N_1.elem=new double[n*n];
	X.row=n;   X.rank=1;   X.elem=new double[n*1];
	l.row=m;   l.rank=1;   l.elem=new double[m*1];
    v.row=m;   v.rank=1;   v.elem=new double[m*1];
}


void comadj::datainput()
{	
 cout<<"input elements of mat A:"<<endl;

	A.input(); 
 cout<<"input elements of P:"<<endl;
	
	P.input();
 cout<<"input elements of observators:"<<endl;
	l.input();
}
comadj::~comadj()
{
}
void comadj::dataoutput()
{   

	double u(0);
	u=datadeal();
		cout<<" the calculated result of X is: "<<endl; 
	X.print();
		cout<<" the calculated result of v is: "<<endl;
    v.print();
		cout<<" the calculated result of N-1 is: "<<endl;
   N_1.print();
   cout<<" 单位权中误差:+-"<<u<<endl;
   cout<<" 未知点的参数中误差:"<<endl;
   for(int i=0;i<N_1.row;i++)
	   cout<<"mx("<<i+1<<")=+-"<<u*sqrt(N_1.elem[i+i+i*(N_1.rank-1)])<<endl;

}

double comadj::datadeal()
{
	MAT AT(A.rank,A.row); AT=A.convert();
	MAT N(AT.row,AT.row); 
	MAT ATP(AT.row,P.rank);tim(AT,P,ATP);
	tim(ATP,A,N);int dh;
	for(int i=0;i<k;i++)
	{
	   cout<<"输入第"<<i+1<<"已知点点号:"<<endl;cin>>dh;
	     N.elem[dh+dh+dh*(N_1.rank-1)]+=1000000;
	}
	inverse(N,N_1);
	MAT ATl(ATP.row,l.rank);tim(ATP,l,ATl);
	tim(N_1,ATl,X);for(i=0;i<X.row*X.rank;i++) X.elem[i]*=-1;
	MAT AX(A.row,X.rank);tim(A,X,AX);
	add(AX,l,v);
	MAT vt(v.rank,v.row);vt=v.convert();
	MAT vtp(vt.row,P.rank);tim(vt,P,vtp);
	MAT S(vtp.row,v.rank);
	tim(vtp,v,S);cout<<"VTPV=   "<<S.elem[0]<<endl;
	double mu;
	mu=sqrt(S.elem[0]/(m-n+k));
    return mu;
}