feat: Would be nice to get diagnostic info, compute residual for solvers in debug mode

[GeometryCollective]

A common bug in geometry processing code is that you solve a linear system Ax=b that's not actually correctly constructed, but forget to check the size of the residual r := b-Ax, which might indicate that something is off. It might be nice to do this check automatically in debug mode.

Also, if there are any internal solver diagnostics (e.g., in SuiteSparse solvers), it would be nice for these to be exposed at the GC level, for debugging purposes.