有限元方法中的数值技术：Crout矩阵分解

如果一个方阵    可以分解为下三角矩阵    和上三角矩阵    的乘积，这种分解称为方阵    的三角分解或    分解。

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

算法原理

对U逐行、对L逐列

• 计算      分解中的      的第一列和      的第一行元素

• 对于      计算      的第      列元素，

• 以及      的第      行元素

Fortran数值实现

subroutine crout(a,l,u,n)
! A = LU
implicitreal*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)
   enddo
   ! L的第k列元素
   l(i,k) = a(i,k) - s
enddo
do j = k+1,n 
   s = 0
   do r = 1,k-1
    s = s + l(k,r) * u(r,j)
   enddo
   ! U的第k行元素
   u(k,j) = (a(k,j)-s)/l(k,k)
enddo
  u(k,k) = 1
enddo
endsubroutine crout

数值案例

分解以下矩阵：

program CroutTest
    ! A = LU
    implicitnone
    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)
    enddo

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

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

100format(3f10.4)

endprogram 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

 

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

参考文献

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