数值分析与多项式研究-Cholesky分解法

Cholesky分解法，即对称正定矩阵的平方根分解法，是由直接三角分解法演化而来，是由于上三角和下三角互为转制而产生的特殊解法

Cholesky分解法同Dollitle分解一样，按列计算，只算下三角（这点与Dollitle不同）

Cholesky分解具有数值稳定性，故不需要选取列主元

#include <stdio.h>
#include <math.h>
#define N 1000
double a[N][N];
double x[N];
double y[N];
int n;
double he(int i,int j)
{
    double sum=0;
    int k;
    if(i==-1)
    {
        for(k=0;k<j;++k)
            sum+=a[j][k]*a[j][k];
        return sum;
    }
    else if(i==-2)
    {
        for(k=0;k<j;++k)
            sum+=a[j][k]*y[k];
        return sum;
    }
    else if(i==-3)
    {
        for(k=j+1;k<n;++k)
            sum+=a[k][j]*x[k];
        return sum;
    }
    else
    {
        for(k=0;k<j;++k)
            sum+=a[i][k]*a[j][k];
        return sum;
    }
}
int main()
{
    int i,j,k;
    //int p,q;
    int flag;
    while(scanf("%d",&n),n)
    {
        flag=1;
        for(i=0;i<n;++i)
            for(j=0;j<n;++j)
                scanf("%lf",&a[i][j]);
        for(k=0;k<n;++k)
            scanf("%lf",&a[k][n]);
        for(i=0;i<n;++i)
        {
            a[i][i]=sqrt(a[i][i]-he((-1),i));
            if(a[i][i]==0)
            {
                flag=0;
                break;
            }
            for(j=i+1;j<n;++j)
                a[j][i]=(a[j][i]-he(j,i))/a[i][i];
            y[i]=(a[i][n]-he((-2),i))/a[i][i];
        }
        if(flag)
        {
            /*for(p=0;p<n;++p)
            {
                for(q=0;q<=p;++q)
                    printf("%lf ",a[p][q]);
                printf("\n");
            }*/
            x[n-1]=y[n-1]/a[n-1][n-1];
            for(i=n-2;i>=0;--i)
                x[i]=(y[i]-he((-3),i))/a[i][i];
            for(i=0;i<n;++i)
                printf("%lf ",x[i]);
        }
        else
            printf("no unique solution");
        printf("\n");
    }
    return 0;
}