@lindahua It's tempting to just change everything to N*sqrt(eps()) or something like that, but unfortunately even this is insufficient for the general case. Even the error bound on the linear solver looks something like `cond(A)*norm(b-Ax)` on the right hand side, which could still blow up every now and then when we get a nearly singular matrix. It's not clear immediately to me if we can just divide all our tests by the norm of the residual and hope that this is good enough.