! lup_decomp.f90  permuting rows for maximum accuracy
! solve for x,  A x = b, using L U = A decomposition
!
! | L11 0   0   0   |   | 1 U12 U13 U14 |   | A11 A12 A13 A14 |
! | L21 L22 0   0   | * | 0 1   U23 U24 | = | A21 A22 A23 A24 |
! | L31 L32 L33 0   |   | 0 0   1   U34 |   | A31 A32 A33 A34 |
! | L41 L42 L43 L44 |   | 0 0   0   1   |   | A41 A42 A43 A44 |
!
! do m=1,n
!   compute m th column of L
!   compute m th row of U
! end do 
! compute bb using b and L     solve L bb = b  get bb   easy
! compute x  using bb and U    solve U x = bb  get x    easy
!
! the matrix U and L can over write A (then no checking possible)
! the vector bb can over write b (then no checking possible)

program lup_decomp
  implicit none
  double precision, dimension(4,4) :: A
  double precision, dimension(4,4) :: U, L
  double precision, dimension(4) :: x, b, bb
  double precision :: t, pivot
  integer :: i, j, k, n, m, im, it
  
  A(1,1) = 2.0
  A(1,2) = 3.0
  A(1,3) = 5.0
  A(1,4) = 6.0
  A(2,1) = 3.0
  A(2,2) = 5.0
  A(2,3) = 6.0
  A(2,4) = 8.0
  A(3,1) = 1.0
  A(3,2) = 4.0
  A(3,3) = 5.0
  A(3,4) = 7.0
  A(4,1) = 3.0
  A(4,2) = 4.0
  A(4,3) = 7.0
  A(4,4) = 8.0
  b(1) = 0.5
  b(2) = 1.5
  b(3) = 4.5
  b(4) = 6.5
  n = 4
  print *, "lup_decomp.f90  n=", n
  
  ! zero U and L, diagonal of U=1
  do i=1,n
    do j= 1,n
      U(i,j) = 0.0
      L(i,j) = 0.0
    end do ! on j
    U(i,i) = 1.0
  end do ! on i

  print *, 'A='
  do i=1,n
    print *,(A(i,j),j=1,n)
  end do ! on i
  print *, 'initial U='
  do i=1,n
    print *,(U(i,j),j=1,n)
  end do ! on i
  print *, 'initial L='
  do i=1,n
    print *,(L(i,j),j=1,n)
  end do ! on i

  ! initialize A row interchange
  do m=1,n-1
    pivot = abs(A(m,m))
    im = m
    do i=m+1,n
      if(abs(A(i,m))>pivot) then
        im = i
        pivot = abs(A(i,m))
      end if
    end do ! i
    print *, 'interchange m=', m, ', im=', im    
    do j=1,n
      t = A(im,j)
      A(im,j) = A(m,j)
      A(m,j) = t
    end do ! on j
    t = b(im)
    b(im) = b(m)
    b(m) = t
  end do ! on i

  ! main loop on m, compute column of L and row of U
  do m=1,n

    ! compute column of L
    do i=m,n
      L(i,m) = A(i,m)
      do k=1,m-1
        L(i,m) = L(i,m) - L(i,k)*U(k,m)
        print *, 'm=', m, ', i=', i, ', k=', k, ', L(i,m)=', L(i,m)
      end do ! on k
      print *, 'L(', i, ',', m, '(=', L(i,m)
    end do ! on i
    ! compute row of U
    do j=m+1,n
      U(m,j) = A(m,j)
      do k=1,m-1
        U(m,j) = U(m,j) - L(m,k)*U(k,j)
        print *, 'm=', m, ', j=', j, ', k=', k, ', U(m,j)=', U(m,j)
      end do ! on k
      U(m,j) = U(m,j)/L(m,m)
      print *, 'U(', m, ',', j, '(=', U(m,j)
    end do ! on j
  end do ! on m

  print *, 'computed U='
  do i=1,n
    print *,(U(i,j),j=1,n)
  end do ! on i
  print *, 'computed L='
  do i=1,n
    print *,(L(i,j),j=1,n)
  end do ! on i

  print *, 'list errors in L*U-A = 0'
  do i=1,n
    do j= 1,n
      t = 0.0
      do k=1,n
        t = t + L(i,k)*U(k,j)
      end do ! on k
      if(t .ne. A(i,j)) then
        print *, 'error A(', i, ',', j, ')-L*U=', A(i,j)-t, &
                 'A(i,j)=', A(i,j), 't=', t
      end if
    end do ! on j
  end do ! on i

  ! solve for x given b
  print *, 'b=', b
  do i=1,n
    bb(i) = b(i)
    do k=1,i-1
      bb(i) = bb(i) - L(i,k)*bb(k)
    end do ! on k
    bb(i) = bb(i)/L(i,i)
  end do ! on i
  print *, 'bb=', (bb(i), i=1,n)

  do j=n,1,-1
    x(j) = bb(j)
    do k=j+1,n
      x(j) = x(j) - U(j,k)*x(k)
    end do ! on k
  end do ! on j
  print *, 'x=', x
  
  print *, 'list errors in A*x-b = 0'
  do i=1,n
    t = 0.0
    do j= 1,n
      t = t + A(i,j)*x(j)
    end do ! on j
    if(t .ne. b(i)) then
      print *, 'error Ax-b(', i, ')=', b(i)-t, 'b(i)=', b(i), 't=', t
    end if
  end do ! on i
  
end program lup_decomp