[Scipy-svn] r7023 - in branches/0.9.x/scipy/linalg: . tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Tue Jan 11 07:47:39 EST 2011
Author: rgommers
Date: 2011-01-11 06:47:38 -0600 (Tue, 11 Jan 2011)
New Revision: 7023
Modified:
branches/0.9.x/scipy/linalg/decomp.py
branches/0.9.x/scipy/linalg/tests/test_decomp.py
Log:
BUG: linalg: work around a LAPACK bug in DGGEV (#709)
(backported from r7002)
Modified: branches/0.9.x/scipy/linalg/decomp.py
===================================================================
--- branches/0.9.x/scipy/linalg/decomp.py 2011-01-10 14:45:25 UTC (rev 7022)
+++ branches/0.9.x/scipy/linalg/decomp.py 2011-01-11 12:47:38 UTC (rev 7023)
@@ -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: branches/0.9.x/scipy/linalg/tests/test_decomp.py
===================================================================
--- branches/0.9.x/scipy/linalg/tests/test_decomp.py 2011-01-10 14:45:25 UTC (rev 7022)
+++ branches/0.9.x/scipy/linalg/tests/test_decomp.py 2011-01-11 12:47:38 UTC (rev 7023)
@@ -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)
+ def test_bad_geneig(self):
+ # 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)
More information about the Scipy-svn
mailing list