ENH: Improve "Kalman smoother" in utils.linalg

The function utils.linalg.block_banded_solve computes the Kalman
smoother solution. It is now slightly improved:
  * one less Cholesky solution for a matrix is required because the
    result is stored in the forward recursion
  * plates are handled better as the plates in y does not affect the
    plates of the diagonal and super-diagonal blocks. Thus, it is
    efficient to compute a large number of solutions for the same system
    matrix.