! test_lup_decomp.f90

! test program for solving simultaneous equations and matrix inverse by the
! L U Decomposition method, using the general solution that requires the
! permutations P

! 11/3/96 JSS modified from Fortran 77 version

program test_lup_decomp
  use lup_decomposition          ! bring in the module
  implicit none
  integer, parameter   :: n = 4  ! keep it small, must fit data
  real, dimension(n,n) :: A = &
      reshape( (/4.0, 4.0, 3.0, 3.0,   5.0, 1.0, 1.0, 2.0, &
                 1.0, 4.0, 3.0, 2.0,   1.0, 2.0, 1.0, 4.0/), shape(A))
  real, dimension(n,n) :: A2 = &
      reshape( (/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/), shape(A2))
  real, dimension(n,n) :: A3 = &
      reshape( (/4.0, 4.0, 4.0, 4.0,   4.0, 3.0, 3.0, 3.0, &
                 4.0, 3.0, 2.0, 2.0,   4.0, 3.0, 2.0, 1.0/), shape(A3))

  real, dimension(n,n) :: AA  ! temporary
  real, dimension(n)   :: X   ! computed by solve
  real, dimension(n), parameter :: Y =  (/ 3.5, 2.4, -1.2, 6.1 /)
  real, dimension(n), parameter :: Y2 = (/30.0, 31.0, 34.0, 40.0/)
  real, dimension(n), parameter :: Y3 = (/16.0, 13.0, 11.0, 10.0/)
  real, dimension(n,n) :: L   ! lower triangular matrix computed
  real, dimension(n,n) :: U   ! upper triangular matrix computed
  real, dimension(n)   :: Z   ! scratch vector for checking solve
  real, dimension(n,n) :: P   ! permutation matrix, used for checking
  integer, dimension(n) :: IDX  ! permutation indices, computed

  print *, "LUP Decomposition test, case 1"

  print *, 'A=  ', A
  call lup_decomp(A, L, U, IDX)
  print *, 'U=  ', U
  print *, 'L=  ', L
  print *, 'IDX= ', IDX
  print *, ''

  print *, 'Y=  ', Y
  call lup_solve(L, U, IDX, X, Y)
  print *, 'X=  ', X
  Z = matmul(A, X)
  print *, "Z=  ", Z
  print *, "Z should equal Y"
  call make_perm(P, IDX)  ! IDX vector becomes P matrix
  print *, 'P= ', P
  AA = matmul(P, A) - matmul(L, U)
  print *, 'P*A-L*U=0 ', AA
  print *, ' '


  print *, 'LUP Decomposition test, case 2'

  print *, 'A2= ', A2
  call lup_decomp(A2, L, U, IDX)
  print *, 'U=  ', U
  print *, 'L=  ', L
  print *, 'IDX= ', IDX
  print *, ''

  print *, 'Y2= ', Y2
  call lup_solve(L, U, IDX, X, Y2)
  print *, 'X=  ', X

  Z = matmul(A2, X)
  print *, 'Z= ', Z
  print *, 'Z should equal Y2'
  call make_perm(P, IDX)
  print *, 'P= ', P
  AA = matmul(P, A2) - matmul(L, U)
  print *, 'P*A2-L*U=0 ', AA
  print *, ' '


  print *, 'LUP Decomposition test, case 3'

  print *, 'A3= ', A3
  call lup_decomp(A3, L, U, IDX)
  print *, 'U=  ', U
  print *, 'L=  ', L
  print *, 'IDX= ', IDX
  print *, ''

  print *, 'Y3= ', Y3
  call lup_solve(L, U, IDX, X, Y3)
  print *, 'X=  ', X
  Z = matmul(A3, X)
  print *, 'Z= ', Z
  print *, 'Z should equal Y3'
  call make_perm(P, IDX)
  print *, 'P= ', P
  AA = matmul(P,A3) - matmul(L,U)
  print *, 'P*A3-L*U=0 ', AA
  print *, ' '

  print *, 'finished LUP Decomposition'
end program test_lup_decomp