BUG: sparse.linalg/lgmres: fix behavior on overflow in intermediate c…

…omputations

Also add the corresponding test for gcrotmk

scipy/sparse/linalg/isolve/lgmres.py

@@ -173,6 +173,10 @@ def lgmres(A, b, x0=None, tol=1e-5, maxiter=1000, M=None, callback=None,
                                          outer_v=outer_v,
                                          prepend_outer_v=prepend_outer_v)
             y *= inner_res_0
+            if not np.isfinite(y).all():
+                # Overflow etc. in computation. There's no way to
+                # recover from this, so we have to bail out.
+                raise LinAlgError()
         except LinAlgError:
             # Floating point over/underflow, non-finite result from
             # matmul etc. -- report failure.

scipy/sparse/linalg/isolve/tests/test_gcrotmk.py

@@ -137,6 +137,19 @@ def test_CU(self):
             assert_(count_1 < count_0/2)
             assert_allclose(x1, x0, atol=1e-14)
 
+    def test_denormals(self):
+        # Check that no warnings are emitted if the matrix contains
+        # numbers for which 1/x has no float representation, and that
+        # the solver behaves properly.
+        A = np.array([[1, 2], [3, 4]], dtype=float)
+        A *= 100 * np.nextafter(0, 1)
+
+        b = np.array([1, 1])
+
+        xp, info = gcrotmk(A, b)
+
+        if info == 0:
+            assert_allclose(A.dot(xp), b)
 
 if __name__ == "__main__":
     run_module_suite()