Computing Determinant in Fortran

It may seem that one should go to LAPACK and find a routine that does just this job. However, it turns out that LAPACK does not really have such a routine. So a little bit of detour is needed. If you want to deal with a general, complex matrix, not symmetric, not Hermitian, then the way to go is using Gauss elimination(GE). LAPACK provides routines to find the LU decomposition of a matrix, which is essentially GE and can be used to compute the determinant. Here is the code snippet:

 1 complex(kind=8) function det(N, mat)
 2     implicit none
 3     integer(kind=8), intent(in) :: N 
 4     complex(kind=8), intent(inout), dimension(:,:) :: mat
 5     integer(kind=8) :: i, info
 6     integer, allocatable :: ipiv(:)
 7     real(kind=8) :: sgn
 8 
 9     allocate(ipiv(N))
10 
11     ipiv = 0
12     
13     call zgetrf(N, N, mat, N, ipiv, info)
14     
15     det = ONE
16     do i = 1, N
17         det = det*mat(i, i)
18     end do
19     sgn = ONE
20     do i = 1, N
21         if(ipiv(i) /= i) then
22             sgn = -sgn
23         end if
24     end do
25     det = sgn*det   
26 end function det