Matrix solvers for the Jacobian

[Yurlungur]

A strenuous test in hydro solvers is to run a problem that should be symmetric, without enforcing it, and check that the solution remains symmetric with time. Running this test in riot with real materials on a strenuous problem, I found that the PTE solver can play a key role in symmetry preservation. This experience suggests that, in particular, the matrix inversion for the Jacobian can matter a great deal.

I found that the most radical solution---replacing the matrix inversion with an SVD pseudo-inverse (hacked in to the code) was effective, but much more expensive and quite brittle. In particular, Eigen's psuedo-inverse on CPU functioned well but the SVD in Kokkos-Kernels frequently failed to converge. The default algorithm provided by Kokkos, however, doesn't have any preconditioning. This suggests to me that more attention should be payed to the matrix inversion, perhaps by introducing pivoting to the QR decomposition we use in the Kokkos code-path.