diff --git a/base/exports.jl b/base/exports.jl index dc449d43c84ea..f2e774b1b8b71 100644 --- a/base/exports.jl +++ b/base/exports.jl @@ -655,6 +655,7 @@ export ldltfact, ldltfact!, linreg, + logabsdet, logdet, lu, lufact!, diff --git a/base/linalg.jl b/base/linalg.jl index 8659acc51bdbb..cb8f8d36421f3 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 7b0e828a18d35..35dbf07d97a40 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 e5f64e20940ac..2bac9dfead5e8 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 49581188eed0b..1cc362630a066 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 7f4a4910a4daa..72310a5cdaa0e 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 diff --git a/test/linalg/generic.jl b/test/linalg/generic.jl new file mode 100644 index 0000000000000..e30c13448ef4e --- /dev/null +++ b/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 \ No newline at end of file