Changeset 452 for inundation/ga/storm_surge/pyvolution/sparse.py

Timestamp:

Oct 26, 2004, 11:55:00 PM (20 years ago)

Author:

steve

Message:

Made some quick and dirty additions to sparse.

File:

• inundation/ga/storm_surge/pyvolution/sparse.py (modified) (8 diffs)

inundation/ga/storm_surge/pyvolution/sparse.py

@@ -1 +1 @@
 """Proof of concept sparse matrix code
 """
-
-
+from scipy_base import *
+from cg_solve import conjugate_gradient, VectorShapeError
 
 class Sparse:
@@ -12 +12 @@
         self.M = M
         self.N = N
+        self.shape = (M,N)
         self.A = {}
 
@@ -40 +41 @@
             return 0.0
 
+    def copy(self):
+
+        new = Sparse(self.M,self.N)
+
+        for key in self.A.keys():
+            i, j = key
+
+            new[i,j] = self.A[i,j]
+
+        return new
+
 
     def todense(self):
@@ -49 +61 @@
             for j in range(self.N):
                 if self.A.has_key( (i,j) ):                
-                    D[i, j] = A[ (i,j) ]
+                    D[i, j] = self.A[ (i,j) ]
         return D
 
@@ -67 +79 @@
 
         #Assume numeric types from now on
-        R = zeros(B.shape, Float) #Result
+        R = zeros((self.M,), Float) #Result
         
         if len(B.shape) == 1:
             #Vector
 
+##             print 'B.shape      ',B.shape
+##             print 'self.shape ',self.shape
+            
             assert B.shape[0] == self.N, 'Mismatching dimensions'
 
@@ -78 +93 @@
                 i, j = key
 
-                R[i] += A[key]*B[j]
-
-        else:
-            raise 'Numeric matrix not yet implemented'
+                R[i] += self.A[key]*B[j]
+
+        else:
+            raise ValueError, 'Numeric matrix not yet implemented'
 
         return R
-        
-
+
+    def __add__(self, other):
+        """Add this matrix onto 'other' 
+        """
+
+        from Numeric import array, zeros, Float
+        
+        new = other.copy()
+
+ #       print 'self.shape',self.shape
+ #       print 'other.shape',other.shape
+        
+        for key in self.A.keys():
+            i, j = key
+
+            new[i,j] = new[i,j] + self.A[key]
+
+        return new
+
+
+    def __rmul__(self, other):
+        """Right multiply this matrix with scalar
+        """
+
+        from Numeric import array, zeros, Float
+        
+        if isscalar(other):
+            new = self.copy()
+            #Multiply nonzero elements
+            for key in new.A.keys():
+                i, j = key
+
+                new.A[key] = other*new.A[key]
+        else:
+            raise 'only right multiple with scalar implemented'
+
+#        print 'new.shape',new.shape
+
+        return new
+
+
+    def trans_mult(self, other):
+        """Multiply the transpose of  matrix with 'other' which can be
+        a Numeric vector.
+        """
+
+        from Numeric import array, zeros, Float
+        
+        try:
+            B = array(other)
+        except:
+            print 'FIXME: Only Numeric types implemented so far'
+
+
+        #Assume numeric types from now on
+
+        
+        if len(B.shape) == 1:
+            #Vector
+
+            assert B.shape[0] == self.M, 'Mismatching dimensions'
+
+            R = zeros((self.N,), Float) #Result
+
+            #Multiply nonzero elements
+            for key in self.A.keys():
+                i, j = key
+
+                R[j] += self.A[key]*B[i]
+
+        else:
+            raise 'Can only multiply with 1d array'
+
+        return R
+
+
+    
+    
 if __name__ == '__main__':
 
@@ -122 +213 @@
     assert allclose(u, [6,14,4])
 
+    u = A.trans_mult(v)
+    print u
+    assert allclose(u, [6,6,10])
+
     #Right hand side column
     v = array([[2,4],[3,4],[4,4]])
@@ -130 +225 @@
     #u = A*v[:,1]
     #print u
+    print A.shape
+
+    B = 3*A
+    print B.todense()
+
+    B[1,0] = 2
+
+    C = A+B
+
+    print C.todense()