diff --git a/base/exports.jl b/base/exports.jl index 79f69ce9fb4657..7b7fe53c0d3724 100644 --- a/base/exports.jl +++ b/base/exports.jl @@ -653,6 +653,7 @@ export ldltfact, ldltfact!, linreg, + logabsdet, logdet, lu, lufact!, diff --git a/base/linalg.jl b/base/linalg.jl index 8659acc51bdbbc..cb8f8d36421f34 100644 --- a/base/linalg.jl +++ b/base/linalg.jl @@ -83,6 +83,7 @@ export ldltfact!, ldltfact, linreg, + logabsdet, logdet, lu, lufact, diff --git a/base/linalg/generic.jl b/base/linalg/generic.jl index 7b0e828a18d358..35dbf07d97a40b 100644 --- a/base/linalg/generic.jl +++ b/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)) + diff --git a/base/linalg/lu.jl b/base/linalg/lu.jl index e5f64e20940ac6..2bac9dfead5e8a 100644 --- a/base/linalg/lu.jl +++ b/base/linalg/lu.jl @@ -160,7 +160,7 @@ 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)) @@ -168,13 +168,15 @@ function logdet2{T<:Real,S}(A::LU{T,S}) # return log(abs(det)) and sign(det) 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))) @@ -182,7 +184,7 @@ function logdet{T<:Complex,S}(A::LU{T,S}) 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 diff --git a/doc/stdlib/linalg.rst b/doc/stdlib/linalg.rst index 49581188eed0b8..1cc362630a0667 100644 --- a/doc/stdlib/linalg.rst +++ b/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 diff --git a/test/choosetests.jl b/test/choosetests.jl index 7f4a4910a4daaf..72310a5cdaa0e1 100644 --- a/test/choosetests.jl +++ b/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