Skip to content

Commit

Permalink
BunchKauffman: correct permutation for rook pivoting (fixes #32080) (#…
Browse files Browse the repository at this point in the history
…32108)

(cherry picked from commit 7b70d49)
  • Loading branch information
timholy authored and KristofferC committed May 23, 2019
1 parent 1b72e97 commit ae8fc04
Show file tree
Hide file tree
Showing 2 changed files with 15 additions and 4 deletions.
13 changes: 9 additions & 4 deletions stdlib/LinearAlgebra/src/bunchkaufman.jl
Original file line number Diff line number Diff line change
Expand Up @@ -125,7 +125,7 @@ size(B::BunchKaufman, d::Integer) = size(getfield(B, :LD), d)
issymmetric(B::BunchKaufman) = B.symmetric
ishermitian(B::BunchKaufman) = !B.symmetric

function _ipiv2perm_bk(v::AbstractVector{T}, maxi::Integer, uplo::AbstractChar) where T
function _ipiv2perm_bk(v::AbstractVector{T}, maxi::Integer, uplo::AbstractChar, rook::Bool) where T
require_one_based_indexing(v)
p = T[1:maxi;]
uploL = uplo == 'L'
Expand All @@ -137,11 +137,16 @@ function _ipiv2perm_bk(v::AbstractVector{T}, maxi::Integer, uplo::AbstractChar)
p[i], p[vi] = p[vi], p[i]
i += uploL ? 1 : -1
else # the 2x2 blocks
if rook
p[i], p[-vi] = p[-vi], p[i]
end
if uploL
p[i + 1], p[-vi] = p[-vi], p[i + 1]
vp = rook ? -v[i+1] : -vi
p[i + 1], p[vp] = p[vp], p[i + 1]
i += 2
else # 'U'
p[i - 1], p[-vi] = p[-vi], p[i - 1]
vp = rook ? -v[i-1] : -vi
p[i - 1], p[vp] = p[vp], p[i - 1]
i -= 2
end
end
Expand Down Expand Up @@ -208,7 +213,7 @@ julia> F.U*F.D*F.U' - F.P*A*F.P'
function getproperty(B::BunchKaufman{T}, d::Symbol) where {T<:BlasFloat}
n = size(B, 1)
if d == :p
return _ipiv2perm_bk(getfield(B, :ipiv), n, getfield(B, :uplo))
return _ipiv2perm_bk(getfield(B, :ipiv), n, getfield(B, :uplo), B.rook)
elseif d == :P
return Matrix{T}(I, n, n)[:,invperm(B.p)]
elseif d == :L || d == :U || d == :D
Expand Down
6 changes: 6 additions & 0 deletions stdlib/LinearAlgebra/test/bunchkaufman.jl
Original file line number Diff line number Diff line change
Expand Up @@ -130,6 +130,12 @@ end
@test F\v5 == F\v6[1:5]
end

@testset "issue #32080" begin
A = Symmetric([-5 -9 9; -9 4 1; 9 1 2])
B = bunchkaufman(A, true)
@test B.U * B.D * B.U' A[B.p, B.p]
end

@test_throws DomainError logdet(bunchkaufman([-1 -1; -1 1]))
@test logabsdet(bunchkaufman([8 4; 4 2]; check = false))[1] == -Inf
@test isa(bunchkaufman(Symmetric(ones(0,0))), BunchKaufman) # 0x0 matrix
Expand Down

0 comments on commit ae8fc04

Please sign in to comment.