Cholesky XA = B solver can be improved.

[kmp5VT]

As pointed out by @wavefunction91, the POSV lapack subroutine can be used in cholesky_inverse instead of the GESV subroutine, since the Gramian is a symmetric and positive semi-definite.