program main

  real A(4,4)
  real A2(4,4)
  real AA(4,4)
  real Y(4)
  real X(4)
  real B(4)
  real B2(4)
  real BB(4)
  real L(4,4)
  real U(4,4)
  real Z(4)
  integer IDX(4)
  integer, parameter :: n=4
  integer i, j

  data A /1.0, 3.0, 2.0, 5.0, 4.0, 4.0, 3.0, 3.0, 5.0, 1.0, 1.0, 2.0, 1.0, 4.0, 3.0, 2.0/

  data B /3.5, 2.4, -1.2, 6.1/

  data A2/1.0, 2.0, 3.0, 4.0, 2.0, 2.0, 3.0, 4.0, 3.0, 3.0, 3.0, 4.0, 4.0, 4.0, 4.0, 4.0/

  data B2/30.0, 31.0, 34.0, 40.0/

  print *, 'LU Decomposition test, case 1'

  do i=1,n
    BB(i)=B(i)
    do j=1,n
      AA(i,j)=A(i,j)
    end do
  end do

  print *, 'A=  ', A

  call lu_decomp(n, AA, L, U, IDX)
  print *, 'AA= ', AA
  print *, 'U=  ', U
  print *, 'L=  ', L
  print *, 'IDX= ', IDX

  print *, 'B=  ', B

  call lu_solve(n, AA, IDX, BB)
  print *, 'BB= ', BB

  do i=1,n
    sum=0.0
    do j=1,n
      sum=sum+A(i,j)*BB(j)
    end do
    Z(i)=sum
  end do
  print *, 'Z= ', Z
  print *, 'Z should equal B, alias y'
  print *, ' '

  print *, 'LU Decomposition test, case 2'

  do i=1,n
    BB(i)=B2(i)
    do j=1,n
      AA(i,j)=A2(i,j)
    end do
  end do

  print *, 'A2= ', A2

  call lu_decomp(n, AA, L, U, IDX)
  print *, 'AA= ', AA
  print *, 'U=  ', U
  print *, 'L=  ', L
  print *, 'IDX= ', IDX

  print *, 'B2= ', B2

  call lu_solve(n, AA, IDX, BB)
  print *, 'BB= ', BB

  do i=1,n
    sum=0.0
    do j=1,n
      sum=sum+A(i,j)*BB(j)
    end do
    Z(i)=sum
  end do
  print *, 'Z= ', Z
  print *, 'Z should equal B2, alias y'
  print *, ' '

  print *, 'finished LU Decomposition'
end

subroutine lu_decomp(n, A, L, U, IDX)
  real A(n,n)
  real L(n,n)
  real U(n,n)
  integer IDX(n)
  real D
  real aamax, dum, sum, VV(n)
  integer imax, i, j, k

  print *, 'in LU Decomp'

  D=1.0
  do i=1,n
    aamax=0.0
    do j=1,n
      if (abs(A(i,j)) .gt. aamax) aamax=abs(A(i,j))
    end do
    if (aamax .eq. 0.0) then
      print *, 'matrix is singular'
      return
    end if
    VV(i)=1.0/aamax
  end do

  do j=1,n
    do i=1,j-1
      sum=A(i,j)
      do k=1,i-1
        sum=sum-A(i,k)*A(k,j)
      end do
      A(i,j)=sum
    end do
    aamax=0.0
    do i=j,n
      sum=A(i,j)
      do k=1,j-1
        sum=sum-A(i,k)*A(k,j)
      end do
      A(i,j)=sum
      dum=VV(i)*abs(sum)
      if (dum .ge. aamax) then
        imax=i
        aamax=dum
      end if
    end do

    if (j .ne. imax) then
      do k=1,n
        dum=A(imax,k)
        A(imax,k)=A(i,k)
        A(i,k)=dum
      end do
      D=-D
      VV(imax)=VV(j)
    end if
    IDX(j)=imax
    if (A(j,j) .eq. 0.0) then
      print *, 'matrix is singular'
      return
    end if

    if (j .ne. n) then
      dum=1.0/A(j,j)
      do i=j+1,n
        A(i,j)=A(i,j)*dum
      end do
    end if
  end do

  print *, 'copy modified A into L and U'

  do i=1,n
    do j=1,n
      L(i,j)=0.0
      U(i,j)=0.0
    end do
  end do
  do i=1,n
    L(i,i)=1.0
  end do
  do i=2,n
    do j=1,i
      L(i,j)=A(i,j)
    end do
  end do
  do i=1,n
    do j=i,n
      U(i,j)=A(i,j)
    end do
  end do

  return
end

subroutine lu_solve(n, A, IDX, B)
  real A(n,n)
  real B(n)
  integer IDX(n)
  integer n, i, j, ii, ll
  real sum

  print *, 'in LU Solve'

  ii=0
  do i=1,n
    ll=IDX(i)
    sum=B(ll)
    B(ll)=B(i)
    if (ii .ne. 0) then
      do j=ii,i-1
        sum=sum-A(i,j)*B(j)
      end do
    else if (sum .ne. 0.0) then
      ii=i
    end if
    B(i)=sum
  end do
  do i=n,1,-1
    sum=B(i)
    do j=i+1,n
      sum=sum-A(i,j)*B(j)
    end do
    B(i)=sum/A(i,i)
  end do
  return
end