Merge pull request #11826 from JuliaLang/anj/logabsdet

Fix #3070. Rename logdet2 to logabsdet + export and tests.

base/exports.jl

@@ -655,6 +655,7 @@ export
     ldltfact,
     ldltfact!,
     linreg,
+    logabsdet,
     logdet,
     lu,
     lufact!,

base/linalg.jl

@@ -83,6 +83,7 @@ export
     ldltfact!,
     ldltfact,
     linreg,
+    logabsdet,
     logdet,
     lu,
     lufact,

base/linalg/generic.jl

@@ -533,3 +533,5 @@ end
 det(x::Number) = x
 
 logdet(A::AbstractMatrix) = logdet(lufact(A))
+logabsdet(A::AbstractMatrix) = logabsdet(lufact(A))
+

base/linalg/lu.jl

@@ -160,29 +160,31 @@ function det{T,S}(A::LU{T,S})
     return prod(diag(A.factors)) * (isodd(sum(A.ipiv .!= 1:n)) ? -one(T) : one(T))
 end
 
-function logdet2{T<:Real,S}(A::LU{T,S})  # return log(abs(det)) and sign(det)
+function logabsdet{T<:Real,S}(A::LU{T,S})  # return log(abs(det)) and sign(det)
     n = chksquare(A)
     dg = diag(A.factors)
     s = (isodd(sum(A.ipiv .!= 1:n)) ? -one(T) : one(T)) * prod(sign(dg))
     sum(log(abs(dg))), s
 end
 
 function logdet{T<:Real,S}(A::LU{T,S})
-    d,s = logdet2(A)
+    d,s = logabsdet(A)
     if s < 0
         throw(DomainError())
     end
     d
 end
 
+_mod2pi(x::BigFloat) = mod(x, big(2)*π) # we don't want to export this, but we use it below
+_mod2pi(x) = mod2pi(x)
 function logdet{T<:Complex,S}(A::LU{T,S})
     n = chksquare(A)
     s = sum(log(diag(A.factors)))
     if isodd(sum(A.ipiv .!= 1:n))
         s = Complex(real(s), imag(s)+π)
     end
     r, a = reim(s)
-    a = π-mod2pi(π-a) #Take principal branch with argument (-pi,pi]
+    a = π - _mod2pi(π - a) #Take principal branch with argument (-pi,pi]
     complex(r,a)
 end
 

doc/stdlib/linalg.rst

@@ -87,6 +87,8 @@ Linear algebra functions in Julia are largely implemented by calling functions f
            ``\``            ✓                       ✓             ✓
            ``cond``         ✓                                     ✓
            ``det``          ✓                       ✓             ✓
+           ``logdet``       ✓                       ✓
+           ``logabsdet``    ✓                       ✓
            ``size``         ✓                       ✓
       ================== ====== ======================== =============
 
@@ -564,6 +566,10 @@ Linear algebra functions in Julia are largely implemented by calling functions f
 
    Log of matrix determinant. Equivalent to ``log(det(M))``, but may provide increased accuracy and/or speed.
 
+.. function:: logabsdet(M)
+
+   Log of absolute value of determinant of real matrix. Equivalent to ``(log(abs(det(M))), sign(det(M)))``, but may provide increased accuracy and/or speed.
+
 .. function:: inv(M)
 
    Matrix inverse

test/choosetests.jl

@@ -64,7 +64,7 @@ function choosetests(choices = [])
                    "linalg/lapack", "linalg/triangular", "linalg/tridiag",
                    "linalg/bidiag", "linalg/diagonal",
                    "linalg/pinv", "linalg/givens", "linalg/cholesky",
-                   "linalg/lu", "linalg/symmetric"]
+                   "linalg/lu", "linalg/symmetric", "linalg/generic"]
     if Base.USE_GPL_LIBS
         push!(linalgtests, "linalg/arnoldi")
     end

test/linalg/generic.jl

@@ -0,0 +1,29 @@
+using Base.Test
+
+debug = false
+
+srand(123)
+
+n = 5 # should be odd
+
+for elty in (Int, Rational{BigInt}, Float32, Float64, BigFloat, Complex{Float32}, Complex{Float64}, Complex{BigFloat})
+    if elty <: Int
+        A = rand(-n:n, n, n) + 10I
+    elseif elty <: Rational
+        A = Rational{BigInt}[rand(-n:n)/rand(1:n) for i = 1:n, j = 1:n] + 10I
+    elseif elty <: Real
+        A = convert(Matrix{elty}, randn(n,n)) + 10I
+    else
+        A = convert(Matrix{elty}, complex(randn(n,n), randn(n,n)))
+    end
+
+    debug && println("element type: $elty")
+
+    @test_approx_eq logdet(A) log(det(A))
+    if elty <: Real
+        @test_approx_eq logabsdet(A)[1] log(abs(det(A)))
+        @test logabsdet(A)[2] == sign(abs(det(A)))
+        @test_throws DomainError logdet(convert(Matrix{elty}, -eye(n)))
+        @test logabsdet(convert(Matrix{elty}, -eye(n)))[2] == -1
+    end
+end