initial working version of SolveCholesky and MatrixInverseCholesky

Hi

Here is a working implementation of a Cholesky based linear solver and a matrix inversion based on a Cholesky factor. This allow for effectively precomputing the solutions of linear systems (with psd matrix) without the need of ever storing a matrix inverse.

Currently, there is no gradient for the Cholesky inversion. While writing tests for my implementation, I discovered a possible bug in gradients of matrix inversions, where the matrix is computed from a Cholesky, see #2066 

I included a couple of unit tests. I am already using this code in a research project of mine and it works fine. It is much faster than always re-computing the linear solve, and storing the inverted matrix is not numerically stable in my example.

Let me know what you think, and what is needed to merge this.
clever way to get the upcasted dtype.