Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Store info in Cholesky type #21976

Merged
merged 3 commits into from
May 30, 2017
Merged

Conversation

fredrikekre
Copy link
Member

Implements JuliaLang/LinearAlgebra.jl#422

The first commit is #21595 (bump on that, its pure code deletion), so when that is merged I will rebase.

Tests pass locally, but there are still things that needs to be changed, for instance

julia> A = [1. 2.; 2. 1.];

julia> det(A)
-3.0

julia> det(cholfact(A))
9.0

@ararslan ararslan added the linear algebra Linear algebra label May 19, 2017
and delay throwing for non-positive definiteness
@fredrikekre fredrikekre changed the title [WIP] Store info in Cholesky type Store info in Cholesky type May 25, 2017
@fredrikekre
Copy link
Member Author

I think I will deal with the isposdef part of JuliaLang/LinearAlgebra.jl#422 in a separate PR (it requires more changes than I first thought). Is this news worthy?

@@ -18,6 +18,11 @@
# through the Hermitian and Symmetric views or exact symmetric or Hermitian elements which
# is checked for and an error is thrown if the check fails.

# The internal structure is as follows
# - _chol! return the factor and info without checking positive definiteness
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

returns ?

@@ -99,30 +113,33 @@ function _chol!(A::AbstractMatrix, ::Type{LowerTriangular})
end
end
end
return LowerTriangular(A)
return LowerTriangular(A), convert(BlasInt, 0) # TODO: If we get here, do we know A is pos. def?
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'd think so. What could be wrong?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Okay thanks. Just wanted to make sure, will remove the comment.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I was paranoid so confirmed with experiments:

for N in (10, 100, 1000), i in 1:1000
    A = rand(N, N); B = A'A
    for M in (A, B)
        C, info = invoke(Base.LinAlg._chol!, Tuple{AbstractMatrix, Type{LowerTriangular}}, copy(M), LowerTriangular)
        infoblas = cholfact(Hermitian(M, :L)).info
        info == infoblas == 0 || (info > 0 && infoblas > 0) || error()
    end
    for M in (A, B)
        C, info = invoke(Base.LinAlg._chol!, Tuple{AbstractMatrix, Type{UpperTriangular}}, copy(M), UpperTriangular)
        infoblas = cholfact(Hermitian(M, :U)).info
        info == infoblas == 0 || (info > 0 && infoblas > 0) || error()
    end
end

@KristofferC KristofferC merged commit 9b51302 into JuliaLang:master May 30, 2017
for T in (Float32, Float64, Complex64, Complex128)
A = T[1 2; 2 1]; B = T[1, 1]
C = cholfact(A)
@show typeof(A), typeof(B), typeof(C.factors)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

debugging output left in?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yes, sry

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
linear algebra Linear algebra
Projects
None yet
Development

Successfully merging this pull request may close these issues.

5 participants