Changeset 594

Timestamp:

Nov 18, 2004, 10:53:59 PM (20 years ago)

Author:

steve

Message:

Working on sparse. Sparse_CSR * vertor works, Sparse_CSR * matrix doesn't.

Location:

inundation/ga/storm_surge/pyvolution

Files:

• sparse.py (modified) (4 diffs)
• sparse_ext.c (modified) (5 diffs)
• test_cg_solve.py (modified) (2 diffs)

inundation/ga/storm_surge/pyvolution/sparse.py

@@ -302 +302 @@
             print 'FIXME: Only Numeric types implemented so far'
 
-
-##        R = csr_mv(self,B)
-        
-        #Assume numeric types from now on
+        return csr_mv(self,B) 
+
+
+
+"""
+Setup for C extensions
+"""
+
+
+def csr_mv(self, B):
+
+    from Numeric import zeros, Float
+
+
+    #Assume numeric types from now on
         
-        if len(B.shape) == 0:
-            #Scalar - use __rmul__ method
-            R = B*self
-            
-        elif len(B.shape) == 1:
-            #Vector
-            assert B.shape[0] == self.N, 'Mismatching dimensions'
-
-            R = zeros(self.M, Float) #Result
-            
-            #Multiply nonzero elements
-            for i in range(self.M):
-                for ckey in range(self.row_ptr[i],self.row_ptr[i+1]):
-                    j = self.colind[ckey]
-                    R[i] += self.data[ckey]*B[j]            
-
-        elif len(B.shape) == 2:
-        
-            R = zeros((self.M, B.shape[1]), Float) #Result matrix
-
-            #Multiply nonzero elements
-
-            for col in range(R.shape[1]):
-                #For each column
-                for i in range(self.M):
-                    for ckey in range(self.row_ptr[i],self.row_ptr[i+1]):
-                        j = self.colind[ckey]
-                        R[i, col] += self.data[ckey]*B[j,col]            
-
-        else:
-            raise ValueError, 'Dimension too high: d=%d' %len(B.shape)
-
-        return R
-
-def csr_mv_python(self, B):
-
-    from Numeric import zeros, Float
-    
-    if len(B.shape) == 1:
+    if len(B.shape) == 0:
+        #Scalar - use __rmul__ method
+        R = B*self
+        
+    elif len(B.shape) == 1:
         #Vector
         assert B.shape[0] == self.N, 'Mismatching dimensions'
-
+        
         R = zeros(self.M, Float) #Result
-            
+        
         #Multiply nonzero elements
         for i in range(self.M):
@@ -356 +333 @@
                 j = self.colind[ckey]
                 R[i] += self.data[ckey]*B[j]            
-
-        return R
+                
+    elif len(B.shape) == 2:
+        
+        R = zeros((self.M, B.shape[1]), Float) #Result matrix
+        
+        #Multiply nonzero elements
+        
+        for col in range(R.shape[1]):
+            #For each column
+            for i in range(self.M):
+                for ckey in range(self.row_ptr[i],self.row_ptr[i+1]):
+                    j = self.colind[ckey]
+                    R[i, col] += self.data[ckey]*B[j,col]            
+                    
     else:
         raise ValueError, 'Dimension too high: d=%d' %len(B.shape)
-
+    
     return R
 
@@ -367 +356 @@
 if compile.can_use_C_extension('sparse_ext.c'):
     #Replace python version with c implementations      
-    from sparse_ext import csr_mv
+    from sparse_ext import csr_mv as csr_mv
 
 
@@ -428 +417 @@
 
     print C.todense()
+
+    C = Sparse_CSR(C)
+
+    y = C*[6,14,4]
+
+    print y
+
+    y2 = C*[[6,4],[4,28],[4,8]]
+
+    print y2

inundation/ga/storm_surge/pyvolution/sparse_ext.c

@@ -16 +16 @@
 int _csr_mv(int M,
             double* data, 
-            int* colind,
-            int* row_ptr,
+            long* colind,
+            long* row_ptr,
             double* x,
             double* y) {
                 
-  int i, j, ckey;
+  long i, j, ckey;
 
   for (i=0; i<M; i++ ) 
@@ -40 +40 @@
   
   PyObject *csr_sparse, // input sparse matrix
-    *R;                 // output wrapped vector
+    *xin, *R;                 // output wrapped vector
   
   PyArrayObject 
@@ -48 +48 @@
     *x,               //Input vector array
     *y;               //Return vector array
+
+  double *xdata;
    
   int dimensions[1], M, err;
   
   // Convert Python arguments to C  
-  if (!PyArg_ParseTuple(args, "OO", &csr_sparse, &x)) 
+  if (!PyArg_ParseTuple(args, "OO", &csr_sparse, &xin)) 
     return NULL;
+
+  x = (PyArrayObject*) PyArray_ContiguousFromObject(xin,PyArray_DOUBLE,1,2);
+  if (!x)
+    return NULL;
+
+/*   printf("x.nd = %i\n",x->nd); */
+/*   printf("x.descr->type_num = %i %i\n",x->descr->type_num,PyArray_LONG); */
+/*   printf("x.dimensions[0] = %i\n",x->dimensions[0]); */
+/*   printf("x.data[0] = %g\n",((double*) x->data)[0]); */
+/*   printf("x.data[1] = %g\n",((double*) x->data)[1]); */
+/*   printf("x.data[2] = %g\n",((double*) x->data)[2]); */
+/*   printf("x.data[3] = %g\n",((double*) x->data)[3]); */
+/*   printf("x.data[4] = %g\n",((double*) x->data)[4]); */
+/*   printf("x.data[5] = %g\n",((double*) x->data)[5]); */
+
+  
 
 
@@ -61 +79 @@
     return NULL;  
 
+
   colind = (PyArrayObject*)
     PyObject_GetAttrString(csr_sparse, "colind"); 
   if (!colind) return NULL;    
-  
+
   row_ptr = (PyArrayObject*) 
     PyObject_GetAttrString(csr_sparse, "row_ptr");   
@@ -77 +96 @@
   err = _csr_mv(M,
                 (double*) data -> data, 
-                (int*)    colind -> data,
-                (int*)    row_ptr -> data,
+                (long*)    colind -> data,
+                (long*)    row_ptr -> data,
                 (double*) x -> data,
                 (double*) y -> data);     

inundation/ga/storm_surge/pyvolution/test_cg_solve.py

@@ -81 +81 @@
         """Standard 2d laplacian"""
         
-        n = 50
-        m = 50
+        n = 100
+        m = 100
 
         A = Sparse(m*n, m*n)
@@ -106 +106 @@
         #print 'finish covert'
         b = A*xe
-        x = conjugate_gradient(A,b,b,iprint=10)
+        x = conjugate_gradient(A,b,b,iprint=20)
 
         assert allclose(x,xe)