program LUfactor

      implicit none
      !declare variables
      integer, parameter :: dkind = selected_real_kind(15,307)
      integer :: i,j,n,k
      real (kind=dkind), allocatable, dimension(:,:) :: A, L, U
      real (kind=dkind), allocatable, dimension(:) :: x, b, y

      !Prompt user for the order of the matrix
      print*, "Enter order of matrix"
      read*, n

      !allocate matrices
      allocate(A(n,n))
      allocate(b(n))
      allocate(x(n))
      allocate(L(n,n))
      allocate(U(n,n))
      allocate(y(n))

      !Enter in entries of the given A matrix
      do i = 1,n
        do j = 1, n
        print*,"enter elements of A ", i,j
        read*,A(i,j)
        enddo
      enddo

      !Enter in entries of the column vector b
      do i = 1,n
        print*, "enter column vector b elements ", i
        read*,b(i)
      enddo

      !stop program if one of the diagonal elements = 0
      do k=1,n
        if (A(k,k) == 0.0) then
        print*,'A has no LU factorization'
        stop
        end if
      end do

      !Factor A into L*U
      do k=1,n-1
        A(k+1:n,k) = A(k+1:n,k)/A(k,k)
        forall (i=k+1:n,j=k+1:n)
            A(i,j) = A(i,j) - A(i,k)*A(k,j)
        end forall
      end do


      !Extract L and U from the new A matrix
      L=0
      U=0
      do i=1,n
        L(i,i) = 1
      end do
      do j=1,n
        do i=1,n
            if (i .gt. j) then
            L(i,j) = A(i,j)
            else
            U(i,j) = A(i,j)
            end if
        end do
      end do

      print *, 'L:'
      do i = 1,n
        print *, (L(i,j), j=1,n)
      enddo
      print *, 'U:'
      do i = 1,n
        print *, (U(i,j), j=1,n)
      enddo


      !Forward substitution Ly=b
      do k=1,n
        do i=1,k-1
            b(k) = b(k)-L(k,i)*b(i)
        end do
        y(k) = b(k)/L(k,k)
      end do
      print *, 'y:', (y(k),k=1,n)

      !Back substitution Ux=y
      do k=n,1,-1
        do i=k+1,n
          y(k) = y(k) - U(k,i)*x(i)
        end do
        x(k) = y(k)/U(k,k)
      end do
      print *, 'x:', (x(k),k=1,n)

      deallocate(A)
      deallocate(b)
      deallocate(x)
      deallocate(L)
      deallocate(U)
      deallocate(y)

      end program LUfactor