[Scipy-svn] r4115 - trunk/scipy/sparse/linalg/eigen/lobpcg

[scipy-svn at scipy.org]

Author: rc
Date: 2008-04-08 09:46:52 -0500 (Tue, 08 Apr 2008)
New Revision: 4115

Modified:
   trunk/scipy/sparse/linalg/eigen/lobpcg/info.py
   trunk/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py
Log:
lobpcg fixes


Modified: trunk/scipy/sparse/linalg/eigen/lobpcg/info.py
===================================================================
--- trunk/scipy/sparse/linalg/eigen/lobpcg/info.py	2008-04-08 13:55:37 UTC (rev 4114)
+++ trunk/scipy/sparse/linalg/eigen/lobpcg/info.py	2008-04-08 14:46:52 UTC (rev 4115)
@@ -63,10 +63,10 @@
 The LOBPCG code internally solves eigenproblems of the size 3``m`` on every
 iteration by calling the "standard" dense eigensolver, so if ``m`` is not 
 small enough compared to ``n``, it does not make sense to call the LOBPCG 
-code, but rather one should use the "standard" eigensolver, e.g. symeig 
+code, but rather one should use the "standard" eigensolver, e.g. scipy or symeig 
 function in this case. If one calls the LOBPCG algorithm for 5``m``>``n``, 
-it will most likely break internally, so the code tries to call symeig 
-instead.
+it will most likely break internally, so the code tries to call the standard
+function instead.
 
 It is not that n should be large for the LOBPCG to work, but rather the
 ratio ``n``/``m`` should be large. It you call the LOBPCG code with ``m``=1

Modified: trunk/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py
===================================================================
--- trunk/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py	2008-04-08 13:55:37 UTC (rev 4114)
+++ trunk/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py	2008-04-08 14:46:52 UTC (rev 4115)
@@ -39,7 +39,10 @@
             else:
                 fun = ll.get_lapack_funcs( ['sygv'], arrays = (mtxA,) )[0]
 ##         print fun
-        w, v, info = fun( mtxA, mtxB )
+        if mtxB is None:
+            out = fun( mtxA )
+        else:
+            out = fun( mtxA, mtxB )
 ##         print w
 ##         print v
 ##         print info
@@ -51,14 +54,13 @@
         ii = nm.argsort( w )
         w = w[slice( *select )]
         if eigenvectors:
-                v = out[1][:,ii]
-                v = v[:,slice( *select )]
+            v = out[1][:,ii]
+            v = v[:,slice( *select )]
+            out = w, v, 0
+        else:
+            out = w, 0
 
-    if eigenvectors:
-        out = w, v
-    else:
-        out = w
-    return out
+    return out[:-1]
 
 def pause():
     raw_input()
@@ -549,6 +551,7 @@
     precond = spdiags( ivals, 0, n, n )
 #    precond = None
     tt = time.clock()
+#    B = None
     eigs, vecs = lobpcg( X, A, B, blockVectorY = Y,
                          M = precond,
                          residualTolerance = 1e-4, maxIterations = 40,