线性方程组的LU分解法的FORTRAN程序

设A是一个方块矩阵。A的LU分解是将它分解成如下形式：

其中L和U分别是下三角矩阵和上三角矩阵。

例如对于一个的矩阵，就有

。

一个LDU分解是一个如下形式的分解：

A = LDU

其中D是对角矩阵，L和U是单位三角矩阵（对角线上全是1的三角矩阵）。

LU分解在本质上是高斯消元法的一种表达形式。实质上是将A通过初等行变换变成一个上三角矩阵，其变换矩阵就是一个单位下三角矩阵。这正是所谓的杜尔里特算法（Doolittle algorithm）：从下至上地对矩阵A做初等行变换，将对角线左下方的元素变成零，然后再证明这些行变换的效果等同于左乘一系列单位下三角矩阵，这一系列单位下三角矩阵的乘积的逆就是L矩阵，它也是一个单位下三角矩阵。

这类算法的复杂度一般在左右，对充分消元的分解则不然。

[编辑]杜尔里特算法

对给定的N × N矩阵

A = (a_n,n)

有

A⁽⁰⁾: = A

然后定义对于n = 1,...,N-1的情况如下：

在第n步，消去矩阵A^(n-1)的第n列主对角线下的元素：将A^(n-1)的第n行乘以之后加到第i行上去。其中。

这相当于在A^(n-1)的左边乘上一个单位下三角矩阵：

于是，定义为：设

A⁽ⁿ⁾: = L_nA^{(n − 1)}.

经过N-1轮操作后，所有在主对角线下的系数都为0了，于是我们得到了一个上三角矩阵：A^(N-1)，这时就有：

这时，矩阵A^(N-1) 就是U，。于是我们得到分解：A = LU。

显然，要是算法成立，在每步操作时必须有。如果这一条件不成立，就要将第n行和另一行交换，由此就会出现一个置换矩阵P。这就是为什么一般来说LU分解里会带有一个置换矩阵的原因。

!Author:RainMan
!2009-12-20
!LU分解法解线性方程组AX=B
!该程序从文件输入，格式为
!矩阵的维数是 3
!A矩阵
!4 2 1
!8 7 4
!16 17 15
!B矩阵
!11 34 95
PROGRAM solve_equation
    REAL,DIMENSION(:,:),ALLOCATABLE::A
    REAL,DIMENSION(:),ALLOCATABLE::B
    CHARACTER*(80):: inputfile = 'dengjin.txt'
    INTEGER DIM
    CALL GET_DIM(inputfile,DIM)
    ALLOCATE(A(DIM,DIM),B(DIM))
    CALL INPUT_MATRIX(A,B,DIM,inputfile)
    CALL LU(A,DIM)
    CALL solve(A,B,B,DIM)
    PRINT*,B
END

SUBROUTINE GET_DIM(inputfile,n)
    CHARACTER*(*) inputfile
    CHARACTER*(80) str
    INTEGER n
    OPEN(10,FILE = inputfile)
    READ(10,*)str,n
    CLOSE(10)
END

SUBROUTINE INPUT_MATRIX(A,B,DIM,inputfile)
    CHARACTER*(*) inputfile
    INTEGER DIM
    INTEGER A(DIM,DIM),B(DIM)
    
    OPEN(20,FILE=inputfile)
    READ(20,*)
    READ(20,*)
    READ(20,*) ((A(I,J),J=1,DIM),I=1,DIM)
    READ(20,*)
    READ(20,*) (B(I),I=1,DIM)
    CLOSE(20)
END

SUBROUTINE LU(A,DIM)
    INTEGER DIM
    REAL A(DIM,DIM)
    
    DO I=2,DIM
        DO J=I,DIM
            A(J,I-1) = A(J,I-1)*1.0/A(I-1,I-1)
            m = A(J,I-1)
            DO K=I,DIM
                A(J,K) = A(J,K)-m*A(I-1,K)
            ENDDO
        ENDDO
    ENDDO
END

SUBROUTINE solve(A,B,X,DIM)
    INTEGER:: DIM
    REAL A(DIM,DIM),B(DIM),X(DIM)
!----------LY = B
    DO I=2,DIM
        DO J=1,I-1
            B(I) = B(I) - A(I,J)*B(J)
        ENDDO
    ENDDO

    B(DIM) = B(DIM)/A(DIM,DIM)
!---------UX=Y Y stores in B,same as X
    DO I=DIM-1,1,-1
        DO J=I+1,DIM
            B(I) = B(I) - A(I,J)*B(J)
        ENDDO
        B(I) = B(I)/A(I,I)
    ENDDO
END