ddet Function

function ddet(A) result(x)
  ! compute the determinant of a real matrix using an LU factorization
  real(8), intent(in)  :: A(:, :)
  real(8)              :: x
  integer              :: i
  integer              :: info, n
  integer, allocatable :: ipiv(:)
  real(8), allocatable :: At(:,:)
  n = size(A(1,:))
  call assert_shape(A, [n, n], "det", "A")
  allocate(At(n,n), ipiv(n))
  At = A
  call dgetrf(n, n, At, n, ipiv, info)
  if(info /= 0) then
     print *, "dgetrf returned info =", info
     if (info < 0) then
        print *, "the", -info, "-th argument had an illegal value"
     else
        print *, "U(", info, ",", info, ") is exactly zero; The factorization"
        print *, "has been completed, but the factor U is exactly"
        print *, "singular, and division by zero will occur if it is used"
        print *, "to solve a system of equations."
     end if
     stop 'det: dgetrf error'
  end if
  ! At now contains the LU of the factorization A = PLU
  ! as L has unit diagonal entries, the determinant can be computed
  ! from the product of U's diagonal entries. Additional sign changes
  ! stemming from the permutations P have to be taken into account as well.
  x = 1d0
  do i = 1,n
     if(ipiv(i) /= i) then  ! additional sign change
        x = -x*At(i,i)
     else
        x = x*At(i,i)
     endif
  end do
end function ddet