The bkfact routine calls the sytrf routine in LAPACK (http://www.netlib.org/lapack/lapack-3.1.1/html/dsytrf.f.html). As described in the LAPACK documentation, its Bunch-Kaufman factorization does not actually compute the factors explicitly, but rather expresses them implicitly as a product of sparse factors.
Given this implicit representation, however, you can still solve linear systems, which is exposed in Julia via \. That is, you can do: B = bkfact(A) x = B \ b to solve Ax = b. You can also do inv(B) to get the inverse of A, and det(B) to get the determinant. If there are other operations that would be useful to perform with the Bunch-Kaufman factorization, patches to Julia would be welcome. If there is a good use-case for explicitly computing the individual factors, this could be exposed as well by implementing a computation of the factors from their implicit representation.
