[Scipy-svn] r3904 - trunk/scipy/splinalg/eigen/arpack

[scipy-svn@scip...]

Author: hagberg
Date: 2008-02-08 08:14:23 -0600 (Fri, 08 Feb 2008)
New Revision: 3904

Modified:
   trunk/scipy/splinalg/eigen/arpack/arpack.py
Log:
More careful handling on which eigenvalues and eigenvectors returned for arpack
nonsymmetric case.



Modified: trunk/scipy/splinalg/eigen/arpack/arpack.py
===================================================================
--- trunk/scipy/splinalg/eigen/arpack/arpack.py	2008-02-07 06:09:12 UTC (rev 3903)
+++ trunk/scipy/splinalg/eigen/arpack/arpack.py	2008-02-08 14:14:23 UTC (rev 3904)
@@ -264,26 +264,30 @@
                 i+=1
             i+=1
                       
-        # Now we have k+1 eigenvalues and eigenvectors
+        # Now we have k+1 possible eigenvalues and eigenvectors
         # Return the ones specified by the keyword "which"
+        nreturned=iparam[4] # number of good eigenvalues returned
+        if nreturned==k:    # we got exactly how many eigenvalues we wanted
+            d=d[:k]
+            z=z[:,:k]
+        else:   # we got one extra eigenvalue (likely a cc pair, but which?)
+            # cut at approx precision for sorting
+            rd=np.round(d,decimals=_ndigits[typ])
+            if which in ['LR','SR']:
+                ind=np.argsort(rd.real) 
+            elif which in ['LI','SI']:
+                # for LI,SI ARPACK returns largest,smallest abs(imaginary) why?
+                ind=np.argsort(abs(rd.imag)) 
+            else:
+                ind=np.argsort(abs(rd))
+            if which in ['LR','LM','LI']:
+                d=d[ind[-k:]]
+                z=z[:,ind[-k:]]
+            if which in ['SR','SM','SI']:
+                d=d[ind[:k]]
+                z=z[:,ind[:k]]
 
-        # cut at approx precision for sorting
-        rd=np.round(d,decimals=_ndigits[typ])
-        if which in ['LR','SR']:
-            ind=np.argsort(rd.real) 
-        elif which in ['LI','SI']:
-            # for LI,SI ARPACK returns largest,smallest abs(imaginary) why?
-            ind=np.argsort(abs(rd.imag)) 
-        else:
-            ind=np.argsort(abs(rd))
-        if which in ['LR','LM','LI']:
-            d=d[ind[-k:]]
-            z=z[:,ind[-k:]]
-        if which in ['SR','SM','SI']:
-            d=d[ind[:k]]
-            z=z[:,ind[:k]]
 
-
     else:
         # complex is so much simpler...
         d,z,info =\