教材《误差理论与测量平差基础》第二版武汉大学出版社
P108页的例7-1的运行结果：
源程序：
#define N 5 /*N是观测值个数*/
#define T 3 /*T是必要观测数*/
#include<stdio.h>
#include<math.h>
float Nbb[T][T],Nb[T][T],W[T][1],x[T][1];
main()
{
float D(float a[T][N],float b[N][N],float c[N][T]);
float K(float a[T][N],float b[N][N],float c[N][1]);
float G(float a[T][T]);
float F(float ca[T-1][T-1]);
float DM(float a[1][N],float b[N][N] ,float c[N][1]);
int i,j,m,n;
float B[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;i<N;i++)
for(j=0;j<T;j++)
scanf("%8f",&B[i][j]);
printf("请输入观测值的权阵P[N][N]:\n");
for(i=0;i<N;i++)
for(j=0;j<N;j++)
scanf("%8f",&P[i][j]);
printf("请输入常数C[N][1]:\n");
for(i=0;i<N;i++)
for(j=0;j<1;j++)
scanf("%8f",&C[i][j]);
for(i=0;i<N;i++)
for(j=0;j<T;j++)
BT[j][i]=B[i][j];
g=D(BT, P, B);
h=K(BT, P, C);
f=G(Nbb);
for(i=0;i<T;i++)
for(j=0;j<1;j++)
{
x[i][j]=Nb[i][0]*W[0][j];
for(m=1;m<T;m++)
x[i][j]+=(Nb[i][m]*W[m][j]);
}
for(i=0;i<T;i++)
x[i][0]=x[i][0]/f;
for(i=0;i<N;i++)
for(j=0;j<1;j++)
{
Bx[i][j]=B[i][0]*x[0][j];
for(m=1;m<T;m++)
Bx[i][j]+=(B[i][m]*x[m][j]);
}
for(i=0;i<N;i++)
V[i][0]=(Bx[i][0]-C[i][0]);
for(i=0;i<N;i++)
for(j=0;j<1;j++)
VT[j][i]=V[i][j];
x1=DM(VT,P,V);
x1=x1/(N-T);
printf("参数x[T][1]=\n");
for(i=0;i<T;i++)
printf("%15f",x[i][0]);
printf("\n");
printf("改正数V[N][1]=\n");
for(i=0;i<N;i++)
printf("%15f",V[i][0]);
printf("\n单位权中误差x1=%15f",sqrt(x1));
printf("\n协因数阵Qxx[T][T]:\n");
for(i=0;i<T;i++)
{
for(j=0;j<T;j++)
printf("%15f",Nb[i][j]/f);
printf("\n");
}
}
float G(float a[T][T])
{
int i,j,m,n;
float c[T-1][T-1],y=0;
for(i=0;i<T;i++)
for(j=0;j<T;j++)
{
for(m=0;m<T;m++)
for(n=0;n<T;n++)
{
if(m<i&&n<j)
c[m][n]=a[m][n];
if(m>i&&n<j)
c[m-1][n]=a[m][n];
if(m<i&&n>j)
c[m][n-1]=a[m][n];
if(m>i&&n>j)
c[m-1][n-1]=a[m][n];
}
if((i+j)%2==0)
Nb[j][i]=F(c);
else
Nb[j][i]=(-1)*F(c);
}
for(m=0;m<T;m++)
y+=(a[0][m]*Nb[m][0]);
return (y);
}
float F(float ca[T-1][T-1])
{
int i,j,m,n,s,t,k=1;
float f=1,c,x,sn;
for (i=0,j=0;i<T-1&&j<T-1;i++,j++) {
if (ca[i][j]==0)
{
for (m=i;ca[m][j]==0;m++);
if (m==T-1)
{
sn=0;
return (sn);
}
else
for (n=j;n<T-1;n++)
{
c=ca[i][n];
ca[i][n]=ca[m][n];
ca[m][n]=c;
}
k*=(-1);
}
for (s=T-2;s>i;s--)
{
x=ca[s][j];
for (t=j;t<T-1;t++)
ca[s][t]-=ca[i][t]*(x/ca[i][j]);
}
}
for (i=0;i<T-1;i++)
f*=ca[i][i];
sn=k*f;
return (sn);
}
float D(float a[T][N],float b[N][N] ,float c[N][T]) {
int i,j,m;
float d[T][N];
for(i=0;i<T;i++)
for(j=0;j<N;j++)
{
d[i][j]=a[i][0]*b[0][j];
for(m=1;m<N;m++)
d[i][j]+=(a[i][m]*b[m][j]);
}
for(i=0;i<T;i++)
for(j=0;j<T;j++)
{
Nbb[i][j]=d[i][0]*c[0][j];
for(m=1;m<N;m++)
Nbb[i][j]+=(d[i][m]*c[m][j]);
}
return (Nbb[0][0]);
}
float K(float a[T][N],float b[N][N],float c[N][1]) {
int i,j,m;
float d[T][N];
for(i=0;i<T;i++)
for(j=0;j<N;j++)
{
d[i][j]=a[i][0]*b[0][j];
for(m=1;m<N;m++)
d[i][j]+=(a[i][m]*b[m][j]);
}
for(i=0;i<T;i++)
for(j=0;j<1;j++)
{
W[i][j]=d[i][0]*c[0][j];
for(m=1;m<N;m++)
W[i][j]+=(d[i][m]*c[m][j]);
}
return (W[0][0]);
}
float DM(float a[1][N],float b[N][N] ,float c[N][1]) {
int i,j,m;
float d[1][N],x;
for(i=0;i<1;i++)
for(j=0;j<N;j++)
{
d[i][j]=a[i][0]*b[0][j];
for(m=1;m<N;m++)
d[i][j]+=(a[i][m]*b[m][j]);
}
for(i=0;i<1;i++)
for(j=0;j<1;j++)
{
x=d[i][0]*c[0][j];
for(m=1;m<N;m++)
x+=(d[i][m]*c[m][j]);
}
return (x);
}
程序说明：
1) 用该程序前，根据具体情况输入N和T；
2) 该程序中的(N-T)自由度(即多余观测数)必须大于等于2，不然程序运行时会出错，原因在于求行列式的逆时，有语句for (s=R-2;s>i;s--)，R=1时s=-1。