平差程序.txt

此程序为测量中的平差程序

#include<iostream>
#include<iomanip> 
using namespace std;
 
int main()
{

    //矩阵转置开始
	const int a=5,b=3;
	double m[a][b]={
		{1,0,0},
		{-1,1,0},
		{0,1,0},
		{0,1,-1},
		{0,0,1}};
	double mm[b][a];
	for(int ii=0;ii<a;ii++)
		for(int jj=0;jj<b;jj++)
		{
			mm[jj][ii]=m[ii][jj];
		}
	//矩阵转置结束







    //输入权阵开始
		double pp[a][a]={
			{2.9,0,0,0,0},
			{0,3.7,0,0,0},
			{0,0,2.5,0,0},
			{0,0,0,3.3,0},
			{0,0,0,0,4.0}};
	 
    //输入权阵结束


	//输入L开始
			double ml[a][1]={0,-23,0,14,0};
	//输入L结束








	//第一次矩阵相乘开始

	
	//B^TP开始
	double cc[b][a]={0.0};
	for(int mn=0;mn<b;mn++)
	{
		for(int nn=0;nn<a;nn++)
		{
			for(int nnn=0;nnn<a;nnn++)
			{
				cc[mn][nn]+=mm[mn][nnn]*pp[nnn][nn];//mm[][]=B^T
			}
		}
	}
	//B^TP结束

    //(B^TP)B开始
	double c[b][b]={0.0};
	for(int mnn=0;mnn<b;mnn++)
	{
		for(int nnm=0;nnm<b;nnm++)
		{
			for(int nnn1=0;nnn1<a;nnn1++)
			{
				c[mnn][nnm]+=cc[mnn][nnn1]*m[nnn1][nnm];
			}
		}
	}
	//(B^TP)B结束


	//第一次矩阵相乘结束,下面对C求逆



	 
	
	//矩阵求逆程序片断开始
	const int n=b; 
    int i,j,h,k; 
    double p; 
    double q[n][12];
    for(i=0;i<n;i++)//构造高斯矩阵 
        for(j=0;j<n;j++) 
            q[i][j]=c[i][j]; 
    for(i=0;i<n;i++) 
        for(j=n;j<12;j++) 
		{
			if(i+6==j) 
            q[i][j]=1; 
            else 
            q[i][j]=0;
		}	 
    for(h=k=0;k<n-1;k++,h++)//消去对角线以下的数据 
        for(i=k+1;i<n;i++) 
		{
		    for(int ii=i;ii<n;ii++)
			{
	        	if((q[k][k]==0)&&(q[ii][k]!=0))
				{
			         for(int r=k;r<12;r++)
					 {
				          double rr;
			              rr=q[k][r];
				          q[k][r]=q[ii][r];
				          q[ii][r]=rr;
						  					 
					 }				  
			     break;    
				}
			}
			 
			if(q[i][h]==0) 
            continue; 
            p=q[k][h]/q[i][h]; 
            for(j=0;j<12;j++) 
			{
				q[i][j]*=p; 
                q[i][j]-=q[k][j]; 
			} 
		} 
    for(h=k=n-1;k>0;k--,h--) // 消去对角线以上的数据 
        for(i=k-1;i>=0;i--) 
		{
			if(q[i][h]==0) 
            continue; 
            p=q[k][h]/q[i][h]; 
            for(j=0;j<12;j++) 
			{
				q[i][j]*=p; 
                q[i][j]-=q[k][j];
			}
		} 
     for(i=0;i<n;i++)//将对角线上数据化为1 
	 { 
		 p=1.0/q[i][i]; 
         for(j=0;j<12;j++) 
         q[i][j]*=p;
	 } 
     for(i=0;i<n;i++) //提取逆矩阵 
          for(j=0;j<n;j++) 
              c[i][j]=q[i][j+6]; 
     
	 //矩阵求逆程序片断结束
	 //((B^TP)B)-1结束,此时C=((B^TP)B)-1


     //第二次矩阵相乘开始

	 //CB^T开始
	 double x[b][a]={0.0};
	 for(int y=0;y<b;y++)
	 {
		for(int yy=0;yy<a;yy++)
		{
			for(int yyy=0;yyy<b;yyy++)
			{
				x[y][yy]+=c[y][yyy]*mm[yyy][yy];
			}
		}
	 }
	//CB^T结束


    //(CB^T)P开始
	 double xx[b][a]={0.0};
	for(int y1=0;y1<b;y1++)
	{
		for(int yy1=0;yy1<a;yy1++)
		{
			for(int yyy1=0;yyy1<a;yyy1++)
			{
				xx[y1][yy1]+=x[y1][yyy1]*pp[yyy1][yy1];
			}
		}
	}
	//(CB^T)P结束

    //(CB^T)PL开始
    double t[b][1]={0.0};
	for(int y2=0;y2<b;y2++)
	{
		for(int yy2=0;yy2<1;yy2++)
		{
			for(int yyy2=0;yyy2<a;yyy2++)
			{
				t[y2][yy2]+=xx[y2][yyy2]*ml[yyy2][yy2];
			}
		}
	}
	//(CB^T)PL结束


	//第二次矩阵相乘结束
	for(int w=0;w<b;w++)
	cout<<"\n\nx"<<w+1<<"="<<fixed<<setprecision(8)<<setw(12)<<t[w][0]<<"\n";
	cout<<endl;

    double v[a][1]={0.0};
	for(int w1=0;w1<a;w1++)
	{
		for(int ww1=0;ww1<1;ww1++)
		{
			for(int www1=0;www1<b;www1++)
			{
				v[w1][ww1]+=m[w1][www1]*t[www1][ww1];
			}
            v[w1][ww1]-=ml[w1][ww1];
		}
	}
	for(int ww=0;ww<a;ww++)
	cout<<"\n\nv"<<ww+1<<"="<<fixed<<setprecision(8)<<setw(12)<<v[ww][0]<<"\n";
	cout<<endl;

	 


     return 0;
}