# [Scipy-svn] r7002 - in trunk/scipy/linalg: . tests

scipy-svn@scip... scipy-svn@scip...
Sat Dec 11 15:09:41 CST 2010

```Author: ptvirtan
Date: 2010-12-11 15:09:40 -0600 (Sat, 11 Dec 2010)
New Revision: 7002

Modified:
trunk/scipy/linalg/decomp.py
trunk/scipy/linalg/tests/test_decomp.py
Log:
BUG: linalg: work around a LAPACK bug in DGGEV (#709)

Modified: trunk/scipy/linalg/decomp.py
===================================================================
--- trunk/scipy/linalg/decomp.py	2010-12-11 18:27:02 UTC (rev 7001)
+++ trunk/scipy/linalg/decomp.py	2010-12-11 21:09:40 UTC (rev 7002)
@@ -11,13 +11,13 @@
# April 2010: Functions for LU, QR, SVD, Schur and Cholesky decompositions were
# moved to their own files.  Still in this file are functions for eigenstuff
# and for the Hessenberg form.
-
+
__all__ = ['eig','eigh','eig_banded','eigvals','eigvalsh', 'eigvals_banded',
'hessenberg']

import numpy
from numpy import array, asarray_chkfinite, asarray, diag, zeros, ones, \
-        isfinite, inexact, nonzero, iscomplexobj, cast
+        isfinite, inexact, nonzero, iscomplexobj, cast, flatnonzero, conj

# Local imports
from scipy.linalg import calc_lwork
@@ -28,20 +28,17 @@

_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.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)
-    count = 0
-    conj = numpy.conjugate
-    for i in range(len(ind)//2):
-        v[:, ind[2*i]] = vnew[:, count]
-        v[:, ind[2*i+1]] = conj(vnew[:, count])
-        count += 1
+def _make_complex_eigvecs(w, vin, dtype):
+    """
+    Produce complex-valued eigenvectors from LAPACK DGGEV real-valued output
+    """
+    # - see LAPACK man page DGGEV at ALPHAI
+    v = numpy.array(vin, dtype=dtype)
+    m = (w.imag > 0)
+    m[:-1] |= (w.imag[1:] < 0) # workaround for LAPACK bug, cf. ticket #709
+    for i in flatnonzero(m):
+        v.imag[:,i] = vin[:,i+1]
+        conj(v[:,i], v[:,i+1])
return v

def _geneig(a1, b, left, right, overwrite_a, overwrite_b):

Modified: trunk/scipy/linalg/tests/test_decomp.py
===================================================================
--- trunk/scipy/linalg/tests/test_decomp.py	2010-12-11 18:27:02 UTC (rev 7001)
+++ trunk/scipy/linalg/tests/test_decomp.py	2010-12-11 21:09:40 UTC (rev 7002)
@@ -132,7 +132,7 @@
assert_array_almost_equal(w,exact_w)

-class TestEig(TestCase):
+class TestEig(object):

def test_simple(self):
a = [[1,2,3],[1,2,3],[2,5,6]]
@@ -154,6 +154,16 @@
for i in range(3):
assert_array_almost_equal(dot(transpose(a),v[:,i]),w[i]*v[:,i])

+    def test_simple_complex_eig(self):
+        a = [[1,2],[-2,1]]
+        w,vl,vr = eig(a,left=1,right=1)
+        assert_array_almost_equal(w, array([1+2j, 1-2j]))
+        for i in range(2):
+            assert_array_almost_equal(dot(a,vr[:,i]),w[i]*vr[:,i])
+        for i in range(2):
+            assert_array_almost_equal(dot(conjugate(transpose(a)),vl[:,i]),
+                                      conjugate(w[i])*vl[:,i])
+
def test_simple_complex(self):
a = [[1,2,3],[1,2,3],[2,5,6+1j]]
w,vl,vr = eig(a,left=1,right=1)
@@ -163,29 +173,32 @@
assert_array_almost_equal(dot(conjugate(transpose(a)),vl[:,i]),
conjugate(w[i])*vl[:,i])

+    def _check_gen_eig(self, A, B):
+        A, B = asarray(A), asarray(B)
+        msg = "\n%r\n%r" % (A, B)
+        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, err_msg=msg)
+
+        assert_array_almost_equal(sort(w[isfinite(w)]), sort(wt[isfinite(wt)]),
+                                  err_msg=msg)
+
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)
+        self._check_gen_eig(A, B)

-        # Disable this test, which fails now, and is not really necessary if the above
-        # succeeds ?
-        #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])))
@@ -195,17 +208,31 @@
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)
+        self._check_gen_eig(A, B)

+        # Ticket #709 (strange return values from DGGEV)
+
+        def matrices(omega):
+            c1 = -9 + omega**2
+            c2 = 2*omega
+            A = [[1, 0,  0,  0],
+                 [0, 1,  0,  0],
+                 [0, 0,  c1, 0],
+                 [0, 0,  0, c1]]
+            B = [[0, 0,  1,   0],
+                 [0, 0,  0,   1],
+                 [1, 0,  0, -c2],
+                 [0, 1, c2,   0]]
+            return A, B
+
+        # With a buggy LAPACK, this can fail for different omega on different
+        # machines -- so we need to test several values
+        for k in xrange(100):
+            A, B = matrices(omega=k*5./100)
+            self._check_gen_eig(A, B)
+
def test_not_square_error(self):
"""Check that passing a non-square array raises a ValueError."""
A = np.arange(6).reshape(3,2)

```