! test_allocate.f90

! read data, allocate arrays, do computation, deallocate arrays
! test program for solving simultaneous equations by the
! L U Decomposition method, using the general solution that requires the
! permutations P

! 11/12/96 JSS 

program test_allocate
  use lup_decomposition          ! bring in the module
  implicit none
  integer :: n 
  real, dimension(:,:), allocatable  :: A
  real, dimension(:), allocatable    :: X   ! computed by solve
  real, dimension(:), allocatable    :: Y
  real, dimension(:,:), allocatable  :: L   ! lower triangular matrix computed
  real, dimension(:,:), allocatable  :: U   ! upper triangular matrix computed
  real, dimension(:), allocatable    :: Z   ! scratch vector for checking solve
  integer, dimension(:), allocatable :: IDX  ! permutation indices, computed
  integer :: iostat, stat
  integer :: i, j

  print *, "test_allocate"
  open(unit=11, file="matrix.dat", action="read", iostat=iostat)
  if(iostat /= 0) then
    print *, "can not open file:  matrix.dat"
    stop
  end if

  do
    read(unit=11, fmt="(i1)", iostat=iostat) n
    if(iostat < 0) exit  ! end of file
    if(iostat > 0) cycle ! bad data
    print *, "allocating working arrays of size ", n
    allocate(A(n,n), stat=stat)
    if(stat > 0) print *, "out of memory"
    allocate(L(n,n))
    allocate(U(n,n))
    allocate(X(n))
    allocate(Y(n))
    allocate(Z(n))
    allocate(IDX(n))

    do i=1,n
      read(unit=11, fmt="(9g9.0)", iostat=iostat) (A(i,j),j=1,n)
      if(iostat /= 0) print *, "may have bad data"
    end do

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

    do i = 1,n
      Y(i) = real(i)
    end do
    print *, 'Y=  ', Y
    call lup_solve(L, U, IDX, X, Y)
    print *, 'X=  ', X
    Z = matmul(A, X)
    print *, "Z=  ", Z, "check, should be Y"

    deallocate(A)
    deallocate(L)
    deallocate(U)
    deallocate(X)
    deallocate(Y)
    deallocate(Z)
    deallocate(IDX)
  end do
  print *, 'finished test_allocate'
end program test_allocate