[Scipy-svn] r7005 - in trunk/scipy/sparse/linalg: . isolve

[scipy-svn at scipy.org]

Author: ptvirtan
Date: 2010-12-11 20:21:18 -0600 (Sat, 11 Dec 2010)
New Revision: 7005

Modified:
   trunk/scipy/sparse/linalg/interface.py
   trunk/scipy/sparse/linalg/isolve/utils.py
Log:
ENH: sparse.linalg: cut unnecessary LinearOperator interface overhead for matrices and arrays (#1036)

Modified: trunk/scipy/sparse/linalg/interface.py
===================================================================
--- trunk/scipy/sparse/linalg/interface.py	2010-12-12 01:19:38 UTC (rev 7004)
+++ trunk/scipy/sparse/linalg/interface.py	2010-12-12 02:21:18 UTC (rev 7005)
@@ -200,6 +200,38 @@
 
         return '<%dx%d LinearOperator with %s>' % (M,N,dt)
 
+class MatrixLinearOperator(LinearOperator):
+    def __init__(self, A):
+        LinearOperator.__init__(self, shape=A.shape, dtype=A.dtype,
+                                matvec=None, rmatvec=self.rmatvec)
+        self.matvec = A.dot
+        self.matmat = A.dot
+        self.__mul__ = A.dot
+        self.A = A
+        self.A_conj = None
+
+    def rmatvec(self, x):
+        if self.A_conj is None:
+            self.A_conj = self.A.T.conj()
+        return self.A_conj.dot(x)
+
+class IdentityOperator(LinearOperator):
+    def __init__(self, shape, dtype):
+        LinearOperator.__init__(self, shape=shape, dtype=dtype, matvec=None,
+                                rmatvec=self.rmatvec)
+
+    def matvec(self, x):
+        return x
+
+    def rmatvec(self, x):
+        return x
+
+    def matmat(self, x):
+        return x
+
+    def __mul__(self, x):
+        return x
+
 def aslinearoperator(A):
     """Return A as a LinearOperator.
 
@@ -226,27 +258,11 @@
     elif isinstance(A, np.ndarray) or isinstance(A, np.matrix):
         if A.ndim > 2:
             raise ValueError('array must have rank <= 2')
-
         A = np.atleast_2d(np.asarray(A))
+        return MatrixLinearOperator(A)
 
-        def matvec(v):
-            return np.dot(A, v)
-        def rmatvec(v):
-            return np.dot(A.conj().transpose(), v)
-        def matmat(V):
-            return np.dot(A, V)
-        return LinearOperator(A.shape, matvec, rmatvec=rmatvec,
-                              matmat=matmat, dtype=A.dtype)
-
     elif isspmatrix(A):
-        def matvec(v):
-            return A * v
-        def rmatvec(v):
-            return A.conj().transpose() * v
-        def matmat(V):
-            return A * V
-        return LinearOperator(A.shape, matvec, rmatvec=rmatvec,
-                              matmat=matmat, dtype=A.dtype)
+        return MatrixLinearOperator(A)
 
     else:
         if hasattr(A, 'shape') and hasattr(A, 'matvec'):

Modified: trunk/scipy/sparse/linalg/isolve/utils.py
===================================================================
--- trunk/scipy/sparse/linalg/isolve/utils.py	2010-12-12 01:19:38 UTC (rev 7004)
+++ trunk/scipy/sparse/linalg/isolve/utils.py	2010-12-12 02:21:18 UTC (rev 7005)
@@ -6,7 +6,8 @@
 
 from numpy import asanyarray, asarray, asmatrix, array, matrix, zeros
 
-from scipy.sparse.linalg.interface import aslinearoperator, LinearOperator
+from scipy.sparse.linalg.interface import aslinearoperator, LinearOperator, \
+     IdentityOperator
 
 _coerce_rules = {('f','f'):'f', ('f','d'):'d', ('f','F'):'F',
                  ('f','D'):'D', ('d','f'):'d', ('d','d'):'d',
@@ -116,7 +117,11 @@
             rpsolve = A_.rpsolve
         else:
             rpsolve = id
-        M = LinearOperator(A.shape, matvec=psolve, rmatvec=rpsolve, dtype=A.dtype)
+        if psolve is id and rpsolve is id:
+            M = IdentityOperator(shape=A.shape, dtype=A.dtype)
+        else:
+            M = LinearOperator(A.shape, matvec=psolve, rmatvec=rpsolve,
+                               dtype=A.dtype)
     else:
         M = aslinearoperator(M)
         if A.shape != M.shape: