LU Decomposition

[JohnNichols]

Module mlu
      Implicit None
      Integer, Parameter :: SP = Kind(1d0) ! set I/O real precision
      Private
      Public luban, lusolve
Contains
      Subroutine luban (a, tol, g, h, ip, condinv, detnth)
! By Banachiewicz (1938, hereafter B38) LU decomposition method calculates such
! triangles L=G^T, and U=H  that square B=A^T=G^TH=LU. Partial pivoting
! by column permutation IP(:) is modern addition.
! Within the code a, g correspond to B38 A^T and G^T, so that a=gh holds.
!
! Normal use is for square A, however for RHS l already known
! input of (A|l)^T yields (L|y^T)^T where x in L^Tx=y is solution of Ax=l.
         Real (SP), Intent (In) :: a (:, :)! input matrix A(m,n), n<=m
         Real (SP), Intent (In) :: tol ! tolerance for near zero pivot
         Real (SP), Intent (Out) :: g (size(a,dim=1), size(a,dim=2)) ! L(m,n)
         Real (SP), Intent (Out) :: h (size(a,dim=2), size(a,dim=2)) ! U(n,n)
                                       ! note U columns are permuted
         Real (SP), Intent (Out) :: condinv ! 1/cond(A), 0 for singular A
         Real (SP), Intent (Out) :: detnth  ! sign*Abs(det(A))**(1/n)
         Integer, Intent (Out)   :: ip (size(a, dim=2)) ! columns permutation
!
         Integer :: k, n, j, l, isig
         Real (SP) :: tol0, pivmax, pivmin, piv
!
         n = size (a, dim=2)
         tol0 = Max (tol, 3._SP*epsilon(tol0))! use default for tol=0
!
! Rectangular A and G are permitted under condition:
         If (n > size(a, dim=1) .Or. n < 1) Stop 91
         Forall (k=1:n) ip (k) = k
         h = 0._SP
         g = 0._SP
         isig = 1
         detnth = 0._SP
         pivmax = Maxval (Abs (a(1, :)))
         pivmin = pivmax
!
         Do k = 1, n
! Banachiewicz (1938) Eq. (2.3)
            h(k, ip(k:)) = a(k, ip(k:)) - Matmul(g(k, :k-1), h(:k-1, ip(k:)))
!
! Find row pivot
            j = (Maxloc(Abs(h(k, ip(k:))), dim=1) + k-1)
            If (j /= k) Then ! Swap columns j and k
               isig = - isig ! Change Det(A) sign because of permutation
               l = ip (k)
               ip (k) = ip (j)
               ip (j) = l
            End If
            piv = Abs (h(k, ip(k)))
            pivmax = Max (piv, pivmax) ! Adjust condinv
            pivmin = Min (piv, pivmin)
            If (piv < tol0) Then ! singular matrix
               isig = 0
               pivmax = 1._SP
               Exit
            Else ! Account for pivot contribution to Det(A) sign and value
               If (h(k, ip(k)) < 0._SP) isig = - isig
               detnth = detnth + Log (piv)
            End If
!
! Transposed Banachiewicz (1938) Eq. (2.4)
            g (k+1:, k) = (a(k+1:, ip(k)) - &
               Matmul(g(k+1:, :k-1), h(:k-1, ip(k)))) / h (k, ip(k))
            g (k, k) = 1._SP
         End Do
!
         detnth = isig * Exp (detnth/n)
         condinv = Abs (isig) * pivmin / pivmax
! Test for square A(n,n) by uncommenting below
!         Print *, '|AQ-LU| ',Maxval (Abs(a(:,ip(:))-Matmul(g, h(:,ip(:)))))
      End Subroutine luban
      
      Subroutine lusolve(l,u,ip,x)
! Solves Ax=b system using triangle factors LU=A      
         Real (SP), Intent (In)    :: l (:, :) ! lower triangle matrix L(n,n)
         Real (SP), Intent (In)    :: u (:, :) ! upper triangle matrix U(n,n)
         Integer, Intent (In)      :: ip (:)   ! columns permutation IP(n)
         Real (SP), Intent (InOut) :: x (:, :) ! Input: m sets of RHSs B(n,m), 
                          ! Output: the corresponding sets of unknowns X(n,m)
         
         Integer :: n, m, i, j
        
         n = size(ip)
         m = size(x, dim=2)
         If (n<1.Or.m<1.Or.Any([n,n]/=shape(l)).Or.Any(shape(l)/=shape(u)).Or. &
            n/=size(x,dim=1)) Stop 91
         Do i = 1, m
            Do j = 1, n
               x(j,i) = x(j,i)-dot_product(x(:j-1,i),l(j,:j-1))
            End Do
            Do j = n, 1, -1
               x(j,i) = (x(j,i)-dot_product(x(j+1:,i),u(j,ip(j+1:)))) / &
                  u(j,ip(j))
            End Do
         End Do
      End Subroutine lusolve     
End Module mlu

I found this when I was looking for stuff on LU decomposition. There are four examples in wikipedia and this is the first.  Of course one is MATLAB.