Is there a command or subroutine for LU factorization? [closed]

Question

In MatLab, the command lu(A) gives as output the two matrices L and U, that is, the LU factorization of A. I was wondering whether there is some command in Fortran doing exactly the same. I have not been able to find it anywhere. I found a lot of subroutines of LAPACK which solve linear systems by first performing the LU factorization, but for my purpouses I need to specifically perform the LU factorization and store the L and U matrices.

Is there a command or subroutine which has as input a square matrix A and as outputs the matrices L and U of the LU factorization?

Answer 1

There is no built-in command that corresponds to lu in Matlab, but we can write a simple wrapper to dgetrf in LAPACK such that the interface is similar to lu, e.g.,

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

Then, we can test the routine with a 3x3 matrix shown in the Matlab page for lu():

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

which gives the expected result:

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

Here please note that the inverse of P is given by its transpose (i.e., inv(P) = P' = transpose(P)) because P is the product of (elementary) permutation matrices.

Answer 2

I have added an method to compute LU using DOLITTLE method. Which is used by MATLAB to computed LU for faster computation involving larger matrices. The algorithm is as follows. To execute the algorithm you have to provide an input file in the format given below. Since the algorithm is a subroutine, you can add it to your code and call it whenever required. Algorithm, input file, output file are as follows.

  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

Here goes the input file format for system Ax=B. First line is number of equations, next three lines are the A matrix element, next three lines are B vector ,

      3
      10 8 3
      3 20 1
      4 5 15
      18
      23
      9  

And the output is generated as,

      ********** 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