determinant Function

    pure function determinant(mat) result(det)
        !! Evaluates the determinant of a Matrix `M`.
        real(dp), intent(in) :: mat(:, :)
        real(dp):: det

        real(dp) :: work_mat(size(mat, 1), size(mat, 2))
        real(dp) :: tmp, c(size(work_mat, dim=1), size(work_mat, dim=2))
        real(dp) :: maxi
        integer:: i, k, l, m, num(size(work_mat, dim=1)), n

        work_mat(:, :) = mat(:, :)

        n = size(work_mat, dim=1)
        det = 1.0_dp
        do k = 1, n
            maxi = work_mat(k, k)
            num(k) = k
            do i = k + 1, n
                if (abs(maxi) < abs(work_mat(i, k))) then
                    maxi = work_mat(i, k)
                    num(k) = i
                end if
            end do
            if (num(k) /= k) then
                do l = k, n
                    tmp = work_mat(k, l)
                    work_mat(k, l) = work_mat(num(k), l)
                    work_mat(num(k), l) = tmp
                end do
                det = -1.0_dp * det
            end if
            do m = k + 1, n
                c(m, k) = work_mat(m, k) / work_mat(k, k)
                do l = k, n
                    work_mat(m, l) = work_mat(m, l) - c(m, k) * work_mat(k, l)
                end do
            end do !There we made matrix triangular!
        end do

        do i = 1, n
            det = det * work_mat(i, i)
        end do
    end function determinant