# Bunch-Kaufman with rook pivoting

I asked this question on Slack but didn’t get an answer. Could someone explain how to use the Bunch-Kaufman factorization with rook pivoting?

``````julia> A = [0.271242 0.657867 1.14586; 0.657867 1.20111 1.13588; 1.14586 1.13588 0.434028]
3×3 Array{Float64,2}:
0.271242  0.657867  1.14586
0.657867  1.20111   1.13588
1.14586   1.13588   0.434028

julia> SA = Symmetric(A, :L);

julia> F = bunchkaufman(SA, true, check=true);  # rook pivoting = true

julia> F.L * F.D * F.L' ≈ A[F.p, F.p]  # this part appears to succeed with rook = false
false

julia> F.L * F.D * F.L' ≈ A
false
``````

How are we to use the permutation to recover the original matrix? Thanks.

It looks like this is a bug:

and the rook-pivoted Bunch–Kaufman factorization is returning the wrong permutation.

Note, however, than most of the time you use the factorization object simply by `x = F \ b` to solve `Ax=b`, and the `\` solver for rook-pivoted Bunch–Kaufman appears to be correct.

1 Like

I noticed that `\` appears to return the correct value (presumably because it calls LAPACK directly instead of relying on the permutation). However, in my use case I specifically need access to `L` and `D`.

In case you haven’t already derived the fix from the LAPACK docs, I think the following may do the job:

``````# adapted from LinearAlgebra stdlib
function _ipiv2p_rook(v::AbstractVector{T}, maxi::Integer, uplo::AbstractChar) where T
p = T[1:maxi;]
uploL = uplo == 'L'
i = uploL ? 1 : maxi
# if uplo == 'U' we construct the permutation backwards
@inbounds while 1 <= i <= length(v)
vi = v[i]
if vi > 0 # the 1x1 blocks
p[i], p[vi] = p[vi], p[i]
i += uploL ? 1 : -1
else # the 2x2 blocks
if uploL
p[i], p[-vi] = p[-vi], p[i]
vp = -v[i+1]
p[i + 1], p[vp] = p[vp], p[i + 1]
i += 2
else # 'U'
p[i], p[-vi] = p[-vi], p[i]
vp = -v[i-1]
p[i - 1], p[vp] = p[vp], p[i - 1]
i -= 2
end
end
end
return p
end
function permvec(B::BunchKaufman{T}) where {T}
B.rook || return B.p
n = size(B,1)
return _ipiv2p_rook(getfield(B, :ipiv), n, getfield(B, :uplo))
end
``````

Then

``````F = bunchkaufman(A,rook)
p = permvec(F)
F.L * F.D * tfunc(F.L) ≈ A[p,p]
``````

where `tfunc` is `adjoint` or `transpose` as appropriate.

5 Likes

A PR for Base to fix #32080 would be welcome.

1 Like

Many thanks @Ralph_Smith! That does the job and appears to match the LAPACK documentation.