Есть ли команда или подпрограмма для факторизации LU?

Question

В MatLab команда lu (A) дает в качестве вывода две матрицы L и U, то есть LU-факторизацию A. Мне было интересно, есть ли какая-то команда в Fortran, выполняющая точно то же самое. Я не смог найти его нигде. Я нашел много подпрограмм LAPACK, которые решают линейные системы, сначала выполнив факторизацию LU, но для моих purpouses мне нужно специально выполнить LU-факторизацию и хранить матрицы L и U.

Есть ли команда или подпрограмма, которая имеет входную квадратную матрицу A и выводит матрицы L и U факторизации LU?

Answer 1

не Там будет не встроенная команда, которая соответствует lu в Matlab, но мы можем написать простую обертку для dgetrf в LAPACK, так что интерфейс похож на 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 

Тогда мы может проверить процедуру с матрицей 3х3, показанной на странице Matlab для лу():

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 является произведением (элементарных) матриц перестановок.

Answer 2

Я добавил метод вычисления LU с использованием метода DOLITTLE. Который используется 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. Первая линия число уравнений, следующие три строки являются матричный элемент А, следующие три линии 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 

Answer 3

Вы можете попробовать «численные рецепты в fortran 77 ", была подпрограмма LU разложения.

Есть много полезных подпрограмм, linalg, stasistics и т. Д.