对称正定矩阵的cholesky分解法程序,版主老大们来救命啊...怎么没有人来帮忙看下啊。。。
这个是一个cholesky分解法的程序,书上的说明:a矩阵是对称正定的系数矩阵,根据a=b*bT(b的转置矩阵)的分解方法求b矩阵的元素。且b矩阵是下三角矩阵,
#include <stdio.h>
#include <math.h>
#define m 3
void main()
{
float a[m][m],b[m][m],c[m][m],i,j,k,sum;
printf("enter a[m][m]:");
for(i=0;i<m;i++)
{ //**输入 4 -1 1
for(j=0;j<m;j++) //** -1 2 -2
{ //** 1 -2 3
scanf("%f",&a[i][j]);
}
}
printf("a[m][m]:");
for(i=0;i<m;i++)
{
for(j=0;j<m;j++)
{
printf("%3.2f ",a[i][j]);
}
printf("\n");
} //**程序执行到了这里就不再执行下去了,大家帮忙找下下面的错误,
for(i=0;i<m;i++)
{
for(j=0;j<m;j++)
{
if(j<=i) b[i][j]=1;
else b[i][j]=0;
}
}
for(i=0;i<m;i++)
{
for(j=0;j<m;j++)
{
c[i][j]=0;
}
}
b[0][0]=sqrt(a[0][0]);
for(i=1,j=0;i<m;i++)
{
b[i][j]=a[i][j]/b[0][0];
}
for(i=1;i<m;i++)
{
for(j=1;j<=i;j++)
{
sum=0;
if(j==i)
{
for(k=0;k<i;k++)
{
sum=b[i][k]*b[i][k]+sum;
}
b[i][j]=sqrt((a[i][j]-sum));
}
else if(j<i)
{
for(k=0;k<j;k++)
{
sum=b[i][k]*b[j][k]+sum;
}
b[i][j]=(a[i][j]-sum)/b[j][j];
}
}
}
for(i=0;i<m;i++)
{
for(j=0;j<=i;k++)
{
c[j][i]=b[i][j];
}
}
printf("b[m][m]:");
for(i=0;i<m;i++)
{
for(j=0;j<m;j++)
{
printf("%6.3f ",b[i][j]);
}
printf("\n");
}
printf("c[m][m]:");
for(i=0;i<m;i++)
{
for(j=0;j<m;j++)
{
printf("%6.3f ",c[i][j]);
}
printf("\n");
}
getch();
}
[ 本帖最后由 李若斌 于 2009-12-23 16:21 编辑 ]