[Scipy-svn] r3182 - in trunk/Lib/linalg: . tests

[scipy-svn at scipy.org]

Author: cdavid
Date: 2007-07-23 07:06:57 -0500 (Mon, 23 Jul 2007)
New Revision: 3182

Modified:
   trunk/Lib/linalg/decomp.py
   trunk/Lib/linalg/tests/test_decomp.py
Log:
Handle Nan eigenvalues in _make_complex_eigvecs, should close #224

Modified: trunk/Lib/linalg/decomp.py
===================================================================
--- trunk/Lib/linalg/decomp.py	2007-07-23 07:33:16 UTC (rev 3181)
+++ trunk/Lib/linalg/decomp.py	2007-07-23 12:06:57 UTC (rev 3182)
@@ -32,7 +32,9 @@
 _I = cast['F'](1j)
 def _make_complex_eigvecs(w,vin,cmplx_tcode):
     v = numpy.array(vin,dtype=cmplx_tcode)
-    ind = numpy.flatnonzero(numpy.not_equal(w.imag,0.0))
+    #ind = numpy.flatnonzero(numpy.not_equal(w.imag,0.0))
+    ind = numpy.flatnonzero(numpy.logical_and(numpy.not_equal(w.imag,0.0),
+                            numpy.isfinite(w)))
     vnew = numpy.zeros((v.shape[0],len(ind)>>1),cmplx_tcode)
     vnew.real = numpy.take(vin,ind[::2],1)
     vnew.imag = numpy.take(vin,ind[1::2],1)

Modified: trunk/Lib/linalg/tests/test_decomp.py
===================================================================
--- trunk/Lib/linalg/tests/test_decomp.py	2007-07-23 07:33:16 UTC (rev 3181)
+++ trunk/Lib/linalg/tests/test_decomp.py	2007-07-23 12:06:57 UTC (rev 3182)
@@ -109,8 +109,47 @@
             assert_array_almost_equal(dot(conjugate(transpose(a)),vl[:,i]),
                                       conjugate(w[i])*vl[:,i])
 
+    def test_singular(self):
+        """Test singular pair"""
+        # Example taken from
+        # http://www.cs.umu.se/research/nla/singular_pairs/guptri/matlab.html
+        A = array(( [22,34,31,31,17], [45,45,42,19,29], [39,47,49,26,34],
+            [27,31,26,21,15], [38,44,44,24,30]))
 
+        B = array(( [13,26,25,17,24], [31,46,40,26,37], [26,40,19,25,25],
+            [16,25,27,14,23], [24,35,18,21,22]))
 
+        w, vr = eig(A,B)
+        wt = eigvals(A,B)
+        val1 = dot(A, vr)
+        val2 = dot(B, vr) * w
+        res = val1 - val2
+        for i in range(res.shape[1]):
+            if all(isfinite(res[:, i])):
+                assert_array_almost_equal(res[:, i], 0)
+        
+        assert_array_almost_equal(w[isfinite(w)], wt[isfinite(w)])
+
+    def test_falker(self):
+        """Test matrices giving some Nan generalized eigen values."""
+        M = diag(array(([1,0,3])))
+        K = array(([2,-1,-1],[-1,2,-1],[-1,-1,2]))
+        D = array(([1,-1,0],[-1,1,0],[0,0,0]))
+        Z = zeros((3,3))
+        I = identity(3)
+        A = bmat([[I,Z],[Z,-K]])
+        B = bmat([[Z,I],[M,D]])
+        A = asarray(A)
+        B = asarray(B)
+
+        w, vr = eig(A,B)
+        val1 = dot(A, vr)
+        val2 = dot(B, vr) * w
+        res = val1 - val2
+        for i in range(res.shape[1]):
+            if all(isfinite(res[:, i])):
+                assert_array_almost_equal(res[:, i], 0)
+
 class test_eig_banded(NumpyTestCase):
 
     def __init__(self, *args):