! permanent.f90  compute permanent of a matrix n!*n^2 flops

double precision function permanent(n, sz, A)
  implicit none
  integer, intent(in) :: n, sz
  double precision, dimension(sz,sz), intent(in) :: A
  double precision, dimension(n,n) :: B
  integer, dimension(n) :: vector, permuter
  double precision :: perm, prod
  integer :: last
  integer :: i, j, k

  interface
    double precision function udrnrt()
    end function udrnrt
    function permute(n, vector, permuter) result(last)
      integer, intent(in) :: n  ! number of items to permute
      integer, dimension(1:*), intent(inout) :: vector, permuter
      integer :: last
    end function permute 
  end interface

  do i=1,n
    vector(i) = i
    permuter(i) = n+1
  end do ! i
  last = 0
  perm = 0.0
  do while(last .eq.0)
    last = permute(n, vector, permuter)
    prod = 1.0
    do i=1,n
      j = vector(i)
      prod = prod * A(i,j)
    end do
    perm = perm + prod
  end do ! while
  permanent = perm
  return
end function permanent