Cholesky分解

# 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);
}