My problem needs to solve dense symmetric linear systems something like:

`A x = b`,
`A y = x`,
`A z = y+x`,...in sequence.

In Eigen C++, if I take advantage of the symmetry of `A` by using:

```
x=A.ldlt().solve(b);
y=A.ldlt().solve(x);
z=A.ldlt().solve(y+x);

```

The same matrix `A` has to be `Cholesky` factorized many times.

Similar issue may exists for `householderQR` and `LU` algorithms.

How can I reuse the coefficient matrix decomposition result for other systems?

It seems this is equivalent to ask whether there is a simple back-substitution solver for upper triangular linear systems.