# Bunch-Kaufman with rook pivoting
**URL:** https://discourse.julialang.org/t/bunch-kaufman-with-rook-pivoting/24379
**Category:** General Usage
**Tags:** linearalgebra
**Created:** [May 19, 2019, 2:39pm UTC](https://discourse.julialang.org/t/bunch-kaufman-with-rook-pivoting/24379 "2019-05-19T14:39:38Z")
**Posts on this page:** 6
**Page:** 1
### Author:  [@dpo](https://discourse.julialang.org/u/dpo)
#### Post date: [May 19, 2019, 2:39pm UTC](https://discourse.julialang.org/t/bunch-kaufman-with-rook-pivoting/24379/1 "2019-05-19T14:39:38Z")
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
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.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [May 19, 2019, 7:38pm UTC](https://discourse.julialang.org/t/bunch-kaufman-with-rook-pivoting/24379/2 "2019-05-19T19:38:13Z")
> [@dpo](#):
>
> How are we to use the permutation to recover the original matrix?
It looks like this is a bug:
>
>
> \#14389 by @timholy added support for rook-pivoted Bunch–Kaufman factorization, b…ut it seems like the permutation (\`p\` or \`P\`) property should have been updated. As reported \[on discourse\](https://discourse.julialang.org/t/bunch-kaufman-with-rook-pivoting/24379), it looks like the wrong permutation. For example:
> \`\`\`
> A = Symmetric(\[-5 -9 9; -9 4 1; 9 1 2\])
> B = bunchkaufman(A, true)
> B.U \* B.D \* B.U' ≈ A\[B.p, B.p\]
> \`\`\`
> gives \`false\`. The permutation \`B.p == \[1,3,2\]\`, whereas it seems that the correct permuation is \`\[2,1,3\]\` in this case.
>
> \#14389 did not change the \[\`\_ipiv2perm\_bk\`\](https://github.com/JuliaLang/julia/blob/ab3ecdfbcb7bceeeb516deb994ae7d00f3537106/stdlib/LinearAlgebra/src/bunchkaufman.jl#L128-L150) routine that computes the permutation from \`ipiv\`, but the documentation for \[\`sytrf\_rook\`\](https://software.intel.com/en-us/mkl-developer-reference-c-sytrf-rook#14BE4059-A295-4A4E-B033-56BA2652E409) seems to indicate that its \`ipiv\` array has a different meaning than that of \[\`sytrf\`\](https://software.intel.com/en-us/mkl-developer-reference-c-sytrf#72C2D4F7-4A77-4C93-8054-4A89F6FE98C4).
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.
---
### Author:  [@dpo](https://discourse.julialang.org/u/dpo)
#### Post date: [May 19, 2019, 8:09pm UTC](https://discourse.julialang.org/t/bunch-kaufman-with-rook-pivoting/24379/3 "2019-05-19T20:09:30Z")
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`.
---
### Author:  [@Ralph\_Smith](https://discourse.julialang.org/u/Ralph_Smith)
#### Post date: [May 20, 2019, 2:46am UTC](https://discourse.julialang.org/t/bunch-kaufman-with-rook-pivoting/24379/4 "2019-05-20T02:46:47Z")
In case you haven’t already derived the fix from the LAPACK docs, I think the following may do the job:
```julia
# 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
```julia
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.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [May 20, 2019, 11:38am UTC](https://discourse.julialang.org/t/bunch-kaufman-with-rook-pivoting/24379/5 "2019-05-20T11:38:25Z")
> [@Ralph\_Smith](#):
>
> I think the following may do the job:
A PR for Base to fix #32080 would be welcome.
---
### Author:  [@dpo](https://discourse.julialang.org/u/dpo)
#### Post date: [May 20, 2019, 2:30pm UTC](https://discourse.julialang.org/t/bunch-kaufman-with-rook-pivoting/24379/6 "2019-05-20T14:30:32Z")
Many thanks @Ralph_Smith! That does the job and appears to match the LAPACK documentation.