New GIZMO 3D matrix inversion algorithm

This MR should (finally) fix all the unit tests in the GIZMO case. testActivePair caused some headaches, since the weird setup for the test (a single pair of cells that interacts, so lots of missing neighbours for most particles) leads to very badly behaved gradient matrices in GIZMO. As a result, most of the matrices that are fed to the matrix inversion function are singular. And while the matrix inversion algorithm can cope with strictly singular matrices (i.e. the matrix determinant is exactly 0), full optimisation tends to mess things up by causing values that are close to but not exactly 0.

I have addressed this in two ways: first, I have replaced the naive matrix inversion algorithm with a safer version that uses the LU decomposition of the matrix. This algorithm can detect singular matrices early on and abort safely if they happen. The matrix inversion function now has a return value that is also used to make decisions in hydro_end_density() and complements the older condition number check. Second, I have implemented a somewhat more relaxed definition of "singular matrix" that uses a limit of 1.e-8 rather than strictly 0 in the relevant condition. This prevents divisions from blowing up and causing numerical overflow later on.

I created a separate MR to make it easier to keep track of this change; this can be merged into unit_test_fixes to actually fix the unit tests.