fortran - 是否有 LU 分解的命令或子程序？

Question

已结束。此问题正在寻求书籍、工具、软件库等的推荐。它不满足 Stack Overflow guidelines 。目前不接受答案。

---

我们不允许提出寻求书籍、工具、软件库等推荐的问题。您可以编辑问题，以便用事实和引文来回答。

已关闭 6 年前。

Improve this question

在MatLab中，命令lu(A)给出两个矩阵L和U作为输出，即A的LU分解。我想知道Fortran中是否有一些命令做完全相同的事情。我一直没能在任何地方找到它。我发现很多 LAPACK 的子例程通过首先执行 LU 分解来求解线性系统，但出于我的目的，我需要专门执行 LU 分解并存储 L 和 U 矩阵。

是否有一个命令或子例程，其输入为方阵 A，输出为 LU 分解的矩阵 L 和 U？

Answer 1

Matlab中没有对应于lu的内置命令，但我们可以为LAPACK中的dgetrf编写一个简单的包装器，其界面类似于lu，例如，

module linalg
    implicit none
    integer, parameter :: dp = kind(0.0d0)
contains
    subroutine lufact( A, L, U, P )
        !... P * A = L * U
        !... http://www.netlib.org/lapack/explore-3.1.1-html/dgetrf.f.html
        !... (note that the definition of P is opposite to that of the above page)

        real(dp), intent(in) :: A(:,:)
        real(dp), allocatable, dimension(:,:) :: L, U, P
        integer, allocatable  :: ipiv(:)
        real(dp), allocatable :: row(:)
        integer :: i, n, info

        n = size( A, 1 )
        allocate( L( n, n ), U( n, n ), P( n, n ), ipiv( n ), row( n ) )

        L = A
        call DGETRF( n, n, L, n, ipiv, info )
        if ( info /= 0 ) stop "lufact: info /= 0"

        U = 0.0d0
        P = 0.0d0
        do i = 1, n
            U( i, i:n ) = L( i, i:n )
            L( i, i:n ) = 0.0d0
            L( i, i )   = 1.0d0
            P( i, i )   = 1.0d0
        enddo

        !... Assuming that P = P[ipiv(n),n] * ... * P[ipiv(1),1]
        !... where P[i,j] is a permutation matrix for i- and j-th rows.
        do i = 1, n
            row = P( i, : )
            P( i, : ) = P( ipiv(i), : )
            P( ipiv(i), : ) = row
        enddo
    endsubroutine
end module

然后，我们可以使用 lu() 的 Matlab 页面中显示的 3x3 矩阵来测试例程:

program test_lu
    use linalg
    implicit none
    real(dp), allocatable, dimension(:,:) :: A, L, U, P

    allocate( A( 3, 3 ) )

    A( 1, : ) = [ 1, 2, 3 ]
    A( 2, : ) = [ 4, 5, 6 ]
    A( 3, : ) = [ 7, 8, 0 ]

    call lufact( A, L, U, P )  !<--> [L,U,P] = lu( A ) in Matlab

    call show( "A =", A )
    call show( "L =", L )
    call show( "U =", U )
    call show( "P =", P )

    call show( "P * A =", matmul( P, A ) )
    call show( "L * U =", matmul( L, U ) )

    call show( "P' * L * U =", matmul( transpose(P), matmul( L, U ) ) )

contains
    subroutine show( msg, X )
        character(*) :: msg
        real(dp) :: X(:,:)
        integer i
        print "(/,a)", trim( msg )
        do i = 1, size(X,1)
            print "(*(f8.4))", X( i, : )
        enddo
    endsubroutine
end program

给出了预期的结果:

A =
  1.0000  2.0000  3.0000
  4.0000  5.0000  6.0000
  7.0000  8.0000  0.0000

L =
  1.0000  0.0000  0.0000
  0.1429  1.0000  0.0000
  0.5714  0.5000  1.0000

U =
  7.0000  8.0000  0.0000
  0.0000  0.8571  3.0000
  0.0000  0.0000  4.5000

P =
  0.0000  0.0000  1.0000
  1.0000  0.0000  0.0000
  0.0000  1.0000  0.0000

P * A =
  7.0000  8.0000  0.0000
  1.0000  2.0000  3.0000
  4.0000  5.0000  6.0000

L * U =
  7.0000  8.0000  0.0000
  1.0000  2.0000  3.0000
  4.0000  5.0000  6.0000

P' * L * U =
  1.0000  2.0000  3.0000
  4.0000  5.0000  6.0000
  7.0000  8.0000  0.0000

请注意，P 的逆是由其转置给出的(即 inv(P) = P' = transpose(P))，因为 P 是(基本)置换矩阵的乘积。