Changeset 8127

Timestamp:

Mar 7, 2011, 5:32:19 PM (14 years ago)

Author:

steve

Message:

Moved kinematic code into a new directory called "operators"

Location:

trunk/anuga_core/source/anuga

Files:

• operators (added)
• operators/kinematic_viscosity.py (added)
• operators/kinematic_viscosity_ext.c (added)
• operators/test_all.py (added)
• operators/test_kinematic_viscosity.py (added)
• utilities/sparse.py (modified) (20 diffs)
• utilities/test_sparse.py (modified) (1 diff)

trunk/anuga_core/source/anuga/utilities/sparse.py

@@ -1 @@
-"""Proof of concept sparse matrix code
+"""
+Proof of concept sparse matrix code
 """
 
@@ -12 +13 @@
         Usage:
 
-        Sparse(A)     #Creates sparse matrix from dense matrix A 
+        Sparse(A)     #Creates sparse matrix from dense matrix A
         Sparse(M, N)  #Creates empty MxN sparse matrix
         """
 
         self.Data = {}
-            
+
         if len(args) == 1:
             try:
@@ -25 +26 @@
 
             assert len(A.shape) == 2, 'Input must be a 2d matrix'
-            
+
             self.M, self.N = A.shape
             for i in range(self.M):
@@ -31 +32 @@
                     if A[i, j] != 0.0:
                         self.Data[i, j] = A[i, j]
-                
-            
+
+
         elif len(args) == 2:
             self.M = args[0]
@@ -38 +39 @@
         else:
             raise Exception('Invalid construction')
-            
-        self.shape = (self.M, self.N) 
+
+        self.shape = (self.M, self.N)
 
 
@@ -52 +53 @@
     def nonzeros(self):
         """Return number of nonzeros of A
-        """        
+        """
         return len(self)
-    
+
     def __setitem__(self, key, x):
 
@@ -60 +61 @@
         # removing these asserts will not speed things up
         assert 0 <= i < self.M
-        assert 0 <= j < self.N        
+        assert 0 <= j < self.N
 
         if x != 0:
             self.Data[key] = float(x)
         else:
-            if self.Data.has_key( key ):            
+            if self.Data.has_key( key ):
                 del self.Data[key]
 
     def __getitem__(self, key):
-        
+
         i,j = key
         # removing these asserts will not speed things up
         assert 0 <= i < self.M
-        assert 0 <= j < self.N                
+        assert 0 <= j < self.N
 
         if self.Data.has_key( key ):
@@ -94 +95 @@
     def todense(self):
         D = num.zeros( (self.M, self.N), num.float)
-        
+
         for i in range(self.M):
             for j in range(self.N):
-                if self.Data.has_key( (i,j) ):                
+                if self.Data.has_key( (i,j) ):
                     D[i, j] = self.Data[ (i,j) ]
         return D
 
 
-    
+
     def __mul__(self, other):
         """Multiply this matrix onto 'other' which can either be
@@ -113 +114 @@
             msg = 'FIXME: Only numeric types implemented so far'
             raise Exception(msg)
-            
+
 
         # Assume numeric types from now on
-        
+
         if len(B.shape) == 0:
             # Scalar - use __rmul__ method
             R = B*self
-            
+
         elif len(B.shape) == 1:
             # Vector
@@ -128 +129 @@
 
             R = num.zeros(self.M, num.float) #Result
-            
+
             # Multiply nonzero elements
             for key in self.Data.keys():
@@ -135 +136 @@
                 R[i] += self.Data[key]*B[j]
         elif len(B.shape) == 2:
-        
-            
+
+
             R = num.zeros((self.M, B.shape[1]), num.float) #Result matrix
 
@@ -142 +143 @@
             for col in range(R.shape[1]):
                 # For each column
-                
+
                 for key in self.Data.keys():
                     i, j = key
 
                     R[i, col] += self.Data[key]*B[j, col]
-            
-            
+
+
         else:
             raise ValueError('Dimension too high: d=%d' %len(B.shape))
 
         return R
-    
+
 
     def __add__(self, other):
-        """Add this matrix onto 'other' 
+        """Add this matrix onto 'other'
         """
 
@@ -220 +221 @@
 class Sparse_CSR:
 
-    def __init__(self, A):
+    def __init__(self, A=None, data=None, Colind=None, rowptr=None, m=None, n=None):
         """Create sparse matrix in csr format.
 
@@ -239 +240 @@
                  Regard it as a pointer into the colind array, for the ith row.
 
-                 
+
         """
 
@@ -264 +265 @@
                 row_ptr[row] = nnz
             #row_ptr[-1] = nnz
-        
+
             self.data    = data
             self.colind  = colind
@@ -270 +271 @@
             self.M       = A.M
             self.N       = A.N
-        else:
-            raise ValueError("Sparse_CSR(A) expects A == Sparse Matrix")
-            
+        elif isinstance(data,num.ndarray) and isinstance(Colind,num.ndarray) and isinstance(rowptr,num.ndarray) and isinstance(m,int) and isinstance(n,int):
+            msg = "Sparse_CSR: data is array of wrong dimensions"
+            #assert len(data.shape) == 1, msg
+            nnz = data.size
+
+            msg = "Sparse_CSR: Colind is array of wrong dimensions"
+            assert Colind.shape == (nnz,), msg
+
+            msg = "Sparse_CSR: rowptr is array of wrong dimensions"
+            assert rowptr.shape == (m+1,), msg
+
+            self.data = data
+            self.colind = Colind
+            self.row_ptr = rowptr
+            self.M = m
+            self.N = n
+        else:
+            raise ValueError('Sparse_CSR(A) expects A == Sparse Matrix *or* data==array,colind==array,rowptr==array,m==int,n==int')
+
     def __repr__(self):
         return '%d X %d sparse matrix:\n' %(self.M, self.N) + `self.data`
@@ -283 +300 @@
     def nonzeros(self):
         """Return number of nonzeros of A
-        """        
+        """
         return len(self)
 
     def todense(self):
         D = num.zeros( (self.M, self.N), num.float)
-        
+
         for i in range(self.M):
             for ckey in range(self.row_ptr[i],self.row_ptr[i+1]):
@@ -305 +322 @@
             print 'FIXME: Only numeric types implemented so far'
 
-        return csr_mv(self,B) 
+        return csr_mv(self,B)
 
 
@@ -317 +334 @@
 if __name__ == '__main__':
     # A little selftest
-    
+
     A = Sparse(3,3)
 
@@ -329 +346 @@
 
     print A
-    print A.todense()    
+    print A.todense()
 
     A[1,2] = 0

trunk/anuga_core/source/anuga/utilities/test_sparse.py

@@ -199 +199 @@
         assert num.allclose(B*C2, [[15.0, 30.0],[10.0, 20.0],[8.0, 16.0],[0.0, 0.0]])
 
+    def test_sparse_csr_init(self):
+        A = num.array([[1.0,0.0,-1.0,0.0],[0.0,2.0,0.0,0.0],[0.0,0.0,0.0,-3.0],[0.0,0.0,4.0,0.0]])
+        data = num.array([1.0,-1.0,2.0,-3.0,4.0])
+        Colind = num.array([0,2,1,3,2])
+        rowptr = num.array([0,2,3,4,5]) #the 5 does not correspond to any row, it is just there so we know how many nonzero entries there are!
+        A_CSR = Sparse_CSR(None,data,Colind,rowptr,4,4)
+        A_dense = A_CSR.todense()
+        assert num.allclose(A, A_dense)
+
 ################################################################################