C语言间接平差程序

教材《误差理论与测量平差基础》第二版武汉大学出版社P108页的例7-1的运行结果：源程序：#defineN5/*N是观测值个数*/#defineT3/*T是必要观测数*/#includestdio.h#includemath.hfloatNbb[T][T],Nb[T][T],W[T][1],x[T][1];main(){floatD(floata[T][N],floatb[N][N],floatc[N][T]);floatK(floata[T][N],floatb[N][N],floatc[N][1]);floatG(floata[T][T]);floatF(floatca[T-1][T-1]);floatDM(floata[1][N],floatb[N][N],floatc[N][1]);inti,j,m,n;floatB[N][T],BT[T][N],V[N][1],VT[1][N],P[N][N],C[N][1],Bx[N][1],f,g,h,x1;printf(请输入V的系数B[N][T]:\n);for(i=0;iN;i++)for(j=0;jT;j++)scanf(%8f,&B[i][j]);printf(请输入观测值的权阵P[N][N]:\n);for(i=0;iN;i++)for(j=0;jN;j++)scanf(%8f,&P[i][j]);printf(请输入常数C[N][1]:\n);for(i=0;iN;i++)for(j=0;j1;j++)scanf(%8f,&C[i][j]);for(i=0;iN;i++)for(j=0;jT;j++)BT[j][i]=B[i][j];g=D(BT,P,B);h=K(BT,P,C);f=G(Nbb);for(i=0;iT;i++)for(j=0;j1;j++){x[i][j]=Nb[i][0]*W[0][j];for(m=1;mT;m++)x[i][j]+=(Nb[i][m]*W[m][j]);}for(i=0;iT;i++)x[i][0]=x[i][0]/f;for(i=0;iN;i++)for(j=0;j1;j++){Bx[i][j]=B[i][0]*x[0][j];for(m=1;mT;m++)Bx[i][j]+=(B[i][m]*x[m][j]);}for(i=0;iN;i++)V[i][0]=(Bx[i][0]-C[i][0]);for(i=0;iN;i++)for(j=0;j1;j++)VT[j][i]=V[i][j];x1=DM(VT,P,V);x1=x1/(N-T);printf(参数x[T][1]=\n);for(i=0;iT;i++)printf(%15f,x[i][0]);printf(\n);printf(改正数V[N][1]=\n);for(i=0;iN;i++)printf(%15f,V[i][0]);printf(\n单位权中误差x1=%15f,sqrt(x1));printf(\n协因数阵Qxx[T][T]:\n);for(i=0;iT;i++){for(j=0;jT;j++)printf(%15f,Nb[i][j]/f);printf(\n);}}floatG(floata[T][T]){inti,j,m,n;floatc[T-1][T-1],y=0;for(i=0;iT;i++)for(j=0;jT;j++){for(m=0;mT;m++)for(n=0;nT;n++){if(mi&&nj)c[m][n]=a[m][n];if(mi&&nj)c[m-1][n]=a[m][n];if(mi&&nj)c[m][n-1]=a[m][n];if(mi&&nj)c[m-1][n-1]=a[m][n];}if((i+j)%2==0)Nb[j][i]=F(c);elseNb[j][i]=(-1)*F(c);}for(m=0;mT;m++)y+=(a[0][m]*Nb[m][0]);return(y);}floatF(floatca[T-1][T-1]){inti,j,m,n,s,t,k=1;floatf=1,c,x,sn;for(i=0,j=0;iT-1&&jT-1;i++,j++){if(ca[i][j]==0){for(m=i;ca[m][j]==0;m++);if(m==T-1){sn=0;return(sn);}elsefor(n=j;nT-1;n++){c=ca[i][n];ca[i][n]=ca[m][n];ca[m][n]=c;}k*=(-1);}for(s=T-2;si;s--){x=ca[s][j];for(t=j;tT-1;t++)ca[s][t]-=ca[i][t]*(x/ca[i][j]);}}for(i=0;iT-1;i++)f*=ca[i][i];sn=k*f;return(sn);}floatD(floata[T][N],floatb[N][N],floatc[N][T]){inti,j,m;floatd[T][N];for(i=0;iT;i++)for(j=0;jN;j++){d[i][j]=a[i][0]*b[0][j];for(m=1;mN;m++)d[i][j]+=(a[i][m]*b[m][j]);}for(i=0;iT;i++)for(j=0;jT;j++){Nbb[i][j]=d[i][0]*c[0][j];for(m=1;mN;m++)Nbb[i][j]+=(d[i][m]*c[m][j]);}return(Nbb[0][0]);}floatK(floata[T][N],floatb[N][N],floatc[N][1]){inti,j,m;floatd[T][N];for(i=0;iT;i++)for(j=0;jN;j++){d[i][j]=a[i][0]*b[0][j];for(m=1;mN;m++)d[i][j]+=(a[i][m]*b[m][j]);}for(i=0;iT;i++)for(j=0;j1;j++){W[i][j]=d[i][0]*c[0][j];for(m=1;mN;m++)W[i][j]+=(d[i][m]*c[m][j]);}return(W[0][0]);}floatDM(floata[1][N],floatb[N][N],floatc[N][1]){inti,j,m;floatd[1][N],x;for(i=0;i1;i++)for(j=0;jN;j++){d[i][j]=a[i][0]*b[0][j];for(m=1;mN;m++)d[i][j]+=(a[i][m]*b[m][j]);}for(i=0;i1;i++)for(j=0;j1;j++){x=d[i][0]*c[0][j];for(m=1;mN;m++)x+=(d[i][m]*c[m][j]);}return(x);}程序说明：1)用该程序前，根据具体情况输入N和T；2)该程序中的(N-T)自由度(即多余观测数)必须大于等于2，不然程序运行时会出错，原因在于求行列式的逆时，有语句for(s=R-2;si;s--)，R=1时s=-1。