Crout矩阵分解¶

约 394 个字 78 行代码预计阅读时间 2 分钟

如果一个方阵\(A\)可以分解为下三角矩阵\(L\)和上三角矩阵\(U\)的乘积，这种分解称为方阵\(A\)的三角分解或\(LU\)分解。

\[A=LU\]

\[ L = \left[ \begin{array} { c c c c } { l _ { 1 1 } } & { 0 } & { \ldots } & { 0 } \\ { l _ { 2 1 } } & { l _ { 2 2 } } & { \ldots } & { 0 } \\ { \vdots } & { \vdots } & { \ddots } & { \vdots } \\ { l _ { n 1 } } & { l _ { n 2 } } & { \ldots } & { l _ { n n } } \end{array} \right] , \ U = \left[ \begin{array} { c c c c } { 1 } & { u _ { 2 1 } } & { \ldots } & { u _ { n 1 } } \\ { 0 } & { 1 } & { \ldots } & { u _ { n 2 } } \\ { \vdots } & { \vdots } & { \ddots } & { \vdots } \\ { 0 } & { 0 } & { \ldots } & { 1 } \end{array} \right] \]

当\(U\)为单位上三角矩阵时，称为 Crout分解；如果\(L\)为单位下三角矩阵时称为 Doolittle分解；若\(U\)是\(L\)的转置（\(A\)为正定对称矩阵），则为 Cholesky分解。

算法原理¶

对U逐行、对L逐列

计算\(LU\)分解中的\(L\)的第一列和\(U\)的第一行元素

\[ \left\{ \begin{array} { l l } { \displaystyle l _ { i 1 } = a _ { i 1 } , i = 1 , 2 \cdots , n } \\ { \displaystyle u _ { 1 j } = \frac { a _ { 1 j } } { l _ { 1 1 } } , j = 2 , \cdots , n ( u _ { 1 1 } = 1 ) } \end{array} \right. \]

对于\(k=2,···,n\)计算\(L\)的第\(k\)列元素，

\[ l _ { i k } = a _ { i k } - \sum _ { r = 1 } ^ { k - 1 } l _ { i r } u _ { r k } , i = k , \cdots , n \]

以及\(U\)的第\(k\)行元素

\[ u _ { k j } = \frac { a _ { k j } - \sum _ { r = 1 } ^ { k - 1 } l _ { k r } u _ { r j } } { l _ { k k } } , j = k + 1 , \cdots , n \]

Fortran数值实现¶

subroutine crout(a,l,u,n)
    ! A = LU
    implicit real*8(a-z)
    integer :: n,i,j,k,r
    real*8 :: a(n,n),l(n,n),u(n,n)

    ! L的第一列和U的第一行
    l(:,1) = a(:,1)
    u(1,:) = a(1,:)/l(1,1)

    do k = 2,n
        do i = k,n
            s = 0
            do r = 1,k-1
                s = s + l(i,r) * u(r,k)
            end do
            ! L的第k列元素
            l(i,k) = a(i,k) - s
        end do
        do j = k+1,n 
            s = 0
            do r = 1,k-1
                s = s + l(k,r) * u(r,j)
            end do
            ! U的第k行元素
            u(k,j) = (a(k,j)-s)/l(k,k)
        end do
        u(k,k) = 1
    end do
end subroutine crout

数值案例¶

分解以下矩阵：

\[ \mathsf { A } = \left[ \begin{array} { c c c } { 2 } & { 1 } & { - 1 } \\ { 4 } & { 1 } & { 0 } \\ { 1 } & { 4 } & { - 1 } \end{array} \right] \]

program CroutTest
    ! A = LU
    implicit none
    integer, parameter :: n = 3
    real*8 :: a(n,n), l(n,n), u(n,n)
    integer :: i, j

    ! 按照行填充
    a = reshape( [ 2.0, 1.0, -1.0,   &  
                  4.0, 1.0,  0.0,   & 
                  1.0, 4.0, -1.0 ], & 
                [n, n],order=[2,1] )

    call crout(a, l, u, n)

    write(*,*) "===== original matrix A ====="
    do i = 1, n
        write(*, 100) (a(i,j), j=1,n)
    end do

    write(*,*) "===== crout decomposition ====="

    write(*,*) "lower triangular matrix L:"
    do i = 1, n
        write(*, 100) (l(i,j), j=1,n)  
    end do

    write(*,*) "upper triangular matrix U:"
    do i = 1, n
        write(*, 100) (u(i,j), j=1,n)
    end do

100 format(3f10.4)

end program CroutTest

输出：

 ===== original matrix A =====
    2.0000    1.0000   -1.0000
    4.0000    1.0000    0.0000
    1.0000    4.0000   -1.0000
 ===== crout decomposition =====
 lower triangular matrix L:
    2.0000    0.0000    0.0000
    4.0000   -1.0000    0.0000
    1.0000    3.5000    6.5000
 upper triangular matrix U:
    1.0000    0.5000   -0.5000
    0.0000    1.0000   -2.0000
    0.0000    0.0000    1.0000

Tip

子程序中默认：implicit real*8(a-z)

参考文献¶

宋叶志,茅永兴,赵秀杰.Fortran 95/2003科学计算与工程[M].清华大学出版社,2011.