#include<stdio.h>
#include<math.h>
#define MAX
10
#define MAXT 15
void cholesky(n, u, x, m)
/*
*/
float n[][MAX],u[],x[];
int m;
{
int i,j,k;
float d[MAX],y,z,w;
for(i=0;i<m;i++)
for(j=0;j<=i;j++)
{ w=n[j][i];
if (i==j)
{for(k=j-1;k>=0;k--)
{y=n[i][k];
z=y*d[k] ;
n[i][k]=z;
w=w-y*z;}
if (w<1e-7) return;
d[i]=1.0/w;}
else
for(k=j-1;k>=0;k--)
{w-=n[i][k]*n[j][k];
n[i][j]=w;}
}
for (j=0;j<m;j++)
{ y=u[j];
for(k=j-1;k>=0;k--)
y-=n[j][k]*x[k];
x[j]=y;}
for(j=m-1;j>=0;j--)
{ y=x[j]*d[j];
for(k=j+1;k<m;k++)
y-=n[k][j]*x[k];
x[j]=y;}
return;}
void
makeN( a,p, n,m,t)
float a[MAX][MAXT],p[],n[][MAX];
int m,t;
{ int i,j,k;
for(i=0;i<m;i++)
for(j=0;j<m;j++)
{ n[i][j]=0.0;
for(k=0;k<t;k++)
n[i][j]+=a[i][k]*p[k]*a[j][k];
}
}
float result(a,p,q,v,m,t)
float a[][MAXT],p[],q[],v[];
int m,t;
{
int i,j;
float s=0.0;
for(i=0;i<t;i++)
{v[i]=0.0;
for(j=0;j<m;j++)
v[i]+=a[j][i]*q[j]*p[i];
}
for (i=0;i<t;i++)
s+=v[i]*v[i]/p[i];
return sqrt(s/m);
}
main()
{ int i,r,n;
static float a[MAX][MAXT]={{1,1,1,0,0,0,0,0},{0,0,-1,0,0,-1,1,0},
{0,-1,0,1,0,0,0,-1},{0,0,0,0,1,1,0,1}}; /*a[MAX]=A*/
static float w[MAX]={28,47,-65,-60};/*w[MAX]=W]*/
static float p[MAXT]={1.81,0.94,1.42,1.76,1.35,0.99,1.38,1.4}; /*p[MAXT]=Q*/
float
k[MAX],v[MAXT],m[MAX][MAX],m0;
printf("please input r,n: " );
scanf("%d,%d",&r,&n); /* r=4,n=8 */
makeN(a,p,m,r,n) ;
cholesky(m,w,k,r);
for(i=0;i<r;i++)printf("%f," ,k[i]);
m0=result(a,p,k,v,r,n);
printf("M0=%f\n\1n",m0);
for(i=0;i<n;i++)
printf("V(%d)=%f\n",i+1,v[i]);
}
一个测量平差的序列。。。VC上的,直接复制c++显示错误。。