# include <math.h>
//n阶实对称方阵的cholesky分解
//a=r*r'
//只用了a的下三角
//计算出的r存在a的下三角
//成功返回0,失败返回非0值
int cholesky(float *a, int n)
{
int i,j,k;
float sum;
for(i=0;i<n;i++) //按行计算
{
for(j=0;j<i+1;j++) //只算下三角
{
sum=a[i*n+j];
for(k=0;k<j;k++)
sum-=a[i*n+k]*a[j*n+k];
if(i==j)
{
if(sum<0)
return (-1);
else
a[i*n+i]=sqrt(sum);
}
else
a[i*n+j]=sum/a[j*n+j];
}
for(j=i+1;j<n;j++) //上三角清零
a[i*n+j]=0;
}
return (0);
}
//n阶实对称方阵的cholesky分解
//a=r*r'
//只用了a的下三角
//计算出的r存在a的下三角
//成功返回0,失败返回非0值
int cholesky(float *a, int n)
{
}