[Scipy-svn] r6273 - trunk/scipy/sparse/linalg/eigen/arpack

[scipy-svn at scipy.org]

Author: cdavid
Date: 2010-03-26 00:35:35 -0500 (Fri, 26 Mar 2010)
New Revision: 6273

Modified:
   trunk/scipy/sparse/linalg/eigen/arpack/arpack.py
Log:
ENH: use implicit A'A computation for sparse SVD.

Modified: trunk/scipy/sparse/linalg/eigen/arpack/arpack.py
===================================================================
--- trunk/scipy/sparse/linalg/eigen/arpack/arpack.py	2010-03-26 05:35:25 UTC (rev 6272)
+++ trunk/scipy/sparse/linalg/eigen/arpack/arpack.py	2010-03-26 05:35:35 UTC (rev 6273)
@@ -519,22 +519,35 @@
     else:
         op = lambda x: x.T
 
-    def _left(x):
+    tp = A.dtype.char
+    linear_at = aslinearoperator(op(A))
+    linear_a = aslinearoperator(A)
+
+    def _left(x, sz):
         x = csc_matrix(x)
-        m = op(x) * x
 
-        eigvals, eigvec = eigen_symmetric(m, k)
+        matvec = lambda x: linear_at.matvec(linear_a.matvec(x))
+        params = _SymmetricArpackParams(sz, k, tp, matvec)
+
+        while not params.converged:
+            params.iterate()
+        eigvals, eigvec = params.extract(True)
         s = np.sqrt(eigvals)
 
         v = eigvec
         u = (x * v) / s
         return u, s, op(v)
 
-    def _right(x):
+    def _right(x, sz):
         x = csr_matrix(x)
-        m = x * op(x)
 
-        eigvals, eigvec = eigen_symmetric(m, k)
+        matvec = lambda x: linear_a.matvec(linear_at.matvec(x))
+        params = _SymmetricArpackParams(sz, k, tp, matvec)
+
+        while not params.converged:
+            params.iterate()
+        eigvals, eigvec = params.extract(True)
+
         s = np.sqrt(eigvals)
 
         u = eigvec
@@ -542,6 +555,6 @@
         return u, s, vh
 
     if n > m:
-        return _left(A)
+        return _left(A, m)
     else:
-        return _right(A)
+        return _right(A, n)