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

Dollitle分解法是建立在Gauss消元法之上的，限制条件少于Gauss消元，即只需要A的各阶顺序主子式不为0即可，又称直接三角分解法

关于L、U的计算方法为先算U的第一行，再算L的第一列，再算U的第二行，再算L的第二列，以此类推

Dollitle分解法的时间复杂度约为三分之一n的立方，与Gauss消元法相当。直接三角分解的优点在于当需要求解具有同系数矩阵的一系列方程组时

（即A相同，而b不相同），可以节省大量计算量，此时在完成三角分解后，并贮存了矩阵L和U之后，右端项每改变一次仅需增加n的平方次运算

其实Dollitle分解法就是记忆了消元后矩阵信息的Gauss消元

并且由于Dollitle分解每一个元素时，只处理了该元素左上子矩阵的元素，故在处理大型带状矩阵时，其有Gaus s消元法无可媲美的优势，可以大量节省存贮空间，建立在其上的Cholesky、Crout分解也具有同样的特点

#include <stdio.h>
#define N 1000
double a[N][N];
double x[N];
double y[N];
int n;
int min(int a,int b)
{
    if(b<a)
        return b;
    return a;
}
double he(int i,int j)
{
    int k;
    double sum=0;
    if(i==-1)
    {
        for(k=0;k<j;++k)
            sum+=y[k]*a[j][k];
        return sum;
    }
    else if(i==-2)
    {
        for(k=j+1;k<n;++k)
            sum+=x[k]*a[j][k];
        return sum;
    }
    else
    {
        if(i==0 || j==0)
            return 0;
        for(k=0;k<min(i,j);++k)
            sum+=a[i][k]*a[k][j];
        return sum;
    }
}
int main()
{
    int i,j,k;
    double tmp;
    int lie;
    double liez;
    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)
        {
            lie=-1;
            liez=-1;
            for(j=i;j<n;++j)
            {
                tmp=a[j][i]-he(j,i);
                if(lie==-1 || tmp>liez)
                {
                    lie=j;
                    liez=tmp;
                }
            }
            if(liez==0 || a[lie][i]==0)
            {
                flag=0;
                break;
            }
            if(lie!=i)
            {
                for(j=0;j<=n;++j)
                {
                    tmp=a[lie][j];
                    a[lie][j]=a[i][j];
                    a[i][j]=tmp;
                }
            }
            for(j=i;j<n;++j)
                a[i][j]=a[i][j]-he(i,j);
            for(j=i+1;j<n;++j)
                a[j][i]=(a[j][i]-he(j,i))/a[i][i];
        }
        if(flag)
        {
            y[0]=a[0][n];
            for(i=1;i<n;++i)
                y[i]=a[i][n]-he((-1),i);
            x[n-1]=y[n-1]/a[n-1][n-1];
            for(i=n-1;i>=0;--i)
                x[i]=(y[i]-he((-2),i))/a[i][i];
            for(i=0;i<n;++i)
                printf("%lf ",x[i]);
        }
        else
            printf("no unique solution");
        printf("\n");
    }
    return 0;
}