是否有用于 LU 分解的命令或子例程？
Is there a command or subroutine for LU factorization?
在 MatLab 中，命令 lu(A) 给出两个矩阵 L 和 U 作为输出，即 A 的 LU 分解。我想知道 Fortran 中是否有一些命令完全相同。我无法在任何地方找到它。我发现了很多 LAPACK 的子程序，它们通过首先执行 LU 分解来求解线性系统，但出于我的目的，我需要专门执行 LU 分解并存储 L 和 U 矩阵。
是否有命令或子程序以方阵 A 作为输入，并以 LU 分解的矩阵 L 和 U 作为输出？
在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是（初等）置换矩阵的乘积。
我添加了一种使用 DOLITTLE 方法计算 LU 的方法。 MATLAB 使用它来计算 LU，以便更快地计算涉及更大的矩阵。算法如下。要执行该算法，您必须提供以下格式的输入文件。由于该算法是一个子例程，您可以将它添加到您的代码中并在需要时调用它。算法、输入文件、输出文件如下
PROGRAM DOLITTLE
IMPLICIT NONE
INTEGER :: n
!**********************************************************
! READS THE NUMBER OF EQUATIONS TO BE SOLVED.
OPEN(UNIT=1,FILE='input.dat',ACTION='READ',STATUS='OLD')
READ(1,*) n
CLOSE(1)
!**********************************************************
CALL LU(n)
END PROGRAM
!==========================================================
! SUBROUTINES TO MAIN PROGRAM
SUBROUTINE LU(n)
IMPLICIT NONE
INTEGER :: i,j,k,p,n,z,ii,itr = 500000
REAL(KIND=8) :: temporary,s1,s2
REAL(KIND=8),DIMENSION(1:n) :: x,b,y
REAL(KIND=8),DIMENSION(1:n,1:n) :: A,U,L,TEMP
REAL(KIND=8),DIMENSION(1:n,1:n+1) :: AB
! READING THE SYSTEM OF EQUATIONS
OPEN(UNIT=2,FILE='input.dat',ACTION='READ',STATUS='OLD')
READ(2,*)
DO I=1,N
READ(2,*) A(I,:)
END DO
DO I=1,N
READ(2,*) B(I)
END DO
CLOSE(2)
DO z=1,itr
U(:,:) = 0
L(:,:) = 0
DO j = 1,n
L(j,j) = 1.0d0
END DO
DO j = 1,n
U(1,j) = A(1,j)
END DO
DO i=2,n
DO j=1,n
DO k=1,i1
s1=0
if (k==1)then
s1=0
else
DO p=1,k1
s1=s1+L(i,p)*U(p,k)
end DO
endif
L(i,k)=(A(i,k)-s1)/U(k,k)
END DO
DO k=i,n
s2=0
DO p=1,i-1
s2=s2+l(i,p)*u(p,k)
END DO
U(i,k)=A(i,k)*s2
END DO
END DO
END DO
IF(z.eq.1)THEN
OPEN(UNIT=3,FILE='output.dat',ACTION='write')
WRITE(3,*) ''
WRITE(3,*) '********** SOLUTIONS *********************'
WRITE(3,*) ''
WRITE(3,*) ''
WRITE(3,*) 'UPPER TRIANGULAR MATRIX'
DO I=1,N
WRITE(3,*) U(I,:)
END DO
WRITE(3,*) ''
WRITE(3,*) ''
WRITE(3,*) 'LOWER TRIANGULAR MATRIX'
DO I=1,N
WRITE(3,*) L(I,:)
END DO
END SUBROUTINE
这是系统 Ax=B 的输入文件格式。第一行是方程数，接下来三行是A矩阵元素，接下来三行是B向量，
3
10 8 3
3 20 1
4 5 15
18
23
9
输出生成为，
********** SOLUTIONS *********************
UPPER TRIANGULAR MATRIX
10.000000000000000 8.0000000000000000 3.0000000000000000
0.0000000000000000 17.600000000000001 0.1000000000000009
0.0000000000000000 0.0000000000000000 13.789772727272727
LOWER TRIANGULAR MATRIX
1.0000000000000000 0.0000000000000000 0.0000000000000000
0.2999999999999999 1.0000000000000000 0.0000000000000000
0.4000000000000002 0.1022727272727273 1.0000000000000000
你可以试试"numerical recipes in fortran 77",
有LU分解子程序。
有很多有用的子程序，linalg，statistics等
在 MatLab 中，命令 lu(A) 给出两个矩阵 L 和 U 作为输出，即 A 的 LU 分解。我想知道 Fortran 中是否有一些命令完全相同。我无法在任何地方找到它。我发现了很多 LAPACK 的子程序，它们通过首先执行 LU 分解来求解线性系统，但出于我的目的，我需要专门执行 LU 分解并存储 L 和 U 矩阵。
是否有命令或子程序以方阵 A 作为输入，并以 LU 分解的矩阵 L 和 U 作为输出？
在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是（初等）置换矩阵的乘积。
我添加了一种使用 DOLITTLE 方法计算 LU 的方法。 MATLAB 使用它来计算 LU，以便更快地计算涉及更大的矩阵。算法如下。要执行该算法，您必须提供以下格式的输入文件。由于该算法是一个子例程，您可以将它添加到您的代码中并在需要时调用它。算法、输入文件、输出文件如下
PROGRAM DOLITTLE
IMPLICIT NONE
INTEGER :: n
!**********************************************************
! READS THE NUMBER OF EQUATIONS TO BE SOLVED.
OPEN(UNIT=1,FILE='input.dat',ACTION='READ',STATUS='OLD')
READ(1,*) n
CLOSE(1)
!**********************************************************
CALL LU(n)
END PROGRAM
!==========================================================
! SUBROUTINES TO MAIN PROGRAM
SUBROUTINE LU(n)
IMPLICIT NONE
INTEGER :: i,j,k,p,n,z,ii,itr = 500000
REAL(KIND=8) :: temporary,s1,s2
REAL(KIND=8),DIMENSION(1:n) :: x,b,y
REAL(KIND=8),DIMENSION(1:n,1:n) :: A,U,L,TEMP
REAL(KIND=8),DIMENSION(1:n,1:n+1) :: AB
! READING THE SYSTEM OF EQUATIONS
OPEN(UNIT=2,FILE='input.dat',ACTION='READ',STATUS='OLD')
READ(2,*)
DO I=1,N
READ(2,*) A(I,:)
END DO
DO I=1,N
READ(2,*) B(I)
END DO
CLOSE(2)
DO z=1,itr
U(:,:) = 0
L(:,:) = 0
DO j = 1,n
L(j,j) = 1.0d0
END DO
DO j = 1,n
U(1,j) = A(1,j)
END DO
DO i=2,n
DO j=1,n
DO k=1,i1
s1=0
if (k==1)then
s1=0
else
DO p=1,k1
s1=s1+L(i,p)*U(p,k)
end DO
endif
L(i,k)=(A(i,k)-s1)/U(k,k)
END DO
DO k=i,n
s2=0
DO p=1,i-1
s2=s2+l(i,p)*u(p,k)
END DO
U(i,k)=A(i,k)*s2
END DO
END DO
END DO
IF(z.eq.1)THEN
OPEN(UNIT=3,FILE='output.dat',ACTION='write')
WRITE(3,*) ''
WRITE(3,*) '********** SOLUTIONS *********************'
WRITE(3,*) ''
WRITE(3,*) ''
WRITE(3,*) 'UPPER TRIANGULAR MATRIX'
DO I=1,N
WRITE(3,*) U(I,:)
END DO
WRITE(3,*) ''
WRITE(3,*) ''
WRITE(3,*) 'LOWER TRIANGULAR MATRIX'
DO I=1,N
WRITE(3,*) L(I,:)
END DO
END SUBROUTINE
这是系统 Ax=B 的输入文件格式。第一行是方程数，接下来三行是A矩阵元素，接下来三行是B向量，
3
10 8 3
3 20 1
4 5 15
18
23
9
输出生成为，
********** SOLUTIONS *********************
UPPER TRIANGULAR MATRIX
10.000000000000000 8.0000000000000000 3.0000000000000000
0.0000000000000000 17.600000000000001 0.1000000000000009
0.0000000000000000 0.0000000000000000 13.789772727272727
LOWER TRIANGULAR MATRIX
1.0000000000000000 0.0000000000000000 0.0000000000000000
0.2999999999999999 1.0000000000000000 0.0000000000000000
0.4000000000000002 0.1022727272727273 1.0000000000000000
你可以试试"numerical recipes in fortran 77", 有LU分解子程序。
有很多有用的子程序，linalg，statistics等