Changeset 605

Timestamp:

Nov 21, 2004, 12:28:38 PM (20 years ago)

Author:

ole

Message:

Implemented matrix-matrix mult in c-extension using CSR format - all tests work again.

Location:

inundation/ga/storm_surge/pyvolution

Files:

• sparse.py (modified) (2 diffs)
• sparse_ext.c (modified) (8 diffs)

inundation/ga/storm_surge/pyvolution/sparse.py

@@ -306 +306 @@
 
 
-"""
-Setup for C extensions
-"""
-
-
 def csr_mv(self, B):
+    """Python version of sparse (CSR) matrix multiplication
+    """
 
     from Numeric import zeros, Float
@@ -353 +350 @@
 
 
+
+#Setup for C extensions
 import compile
 if compile.can_use_C_extension('sparse_ext.c'):
-    #Replace python version with c implementations      
-    from sparse_ext import csr_mv as csr_mv
-
+    #Replace python version with c implementation
+    from sparse_ext import csr_mv
 
 if __name__ == '__main__':

inundation/ga/storm_surge/pyvolution/sparse_ext.c

@@ -15 +15 @@
 #include "stdio.h"
 
+//Matrix-vector routine
 int _csr_mv(int M,
             double* data, 
@@ -24 +25 @@
   long i, j, ckey;
 
-  //printf("test");
   for (i=0; i<M; i++ ) 
     for (ckey=row_ptr[i]; ckey<row_ptr[i+1]; ckey++) {
@@ -31 +31 @@
     }              
   
-  //printf("end");
   return 0;
 }            
 
+//Matrix-matrix routine
+int _csr_mm(int M,
+            int columns, 
+            double* data, 
+            long* colind,
+            long* row_ptr,
+            double* x,
+            double* y) {
+                
+  long i, j, ckey, c, rowind_i, rowind_j;
+
+  for (i=0; i<M; i++ ) {
+    rowind_i = i*columns;
+    
+    for (ckey=row_ptr[i]; ckey<row_ptr[i+1]; ckey++) {
+      j = colind[ckey];
+      rowind_j = j*columns;
+          
+      for (c=0; c<columns; c++) {
+        y[rowind_i+c] += data[ckey]*x[rowind_j+c];
+      }             
+    } 
+  }
+  
+  return 0;
+}            
 
                      
@@ -43 +68 @@
   
   PyObject *csr_sparse, // input sparse matrix
-    *xin, *R;                 // output wrapped vector
+    *xin, *R;           // output wrapped vector
   
   PyArrayObject 
@@ -53 +78 @@
 
   
-  int dimensions[1], M, err;
+  int dimensions[1], M, err, columns, rows;
   
   // Convert Python arguments to C  
@@ -73 +98 @@
 /*   printf("x.data[5] = %g\n",((double*) x->data)[5]); */
 
-  
+ 
 
 
@@ -80 +105 @@
   if (!data) 
     return NULL;  
-
 
   colind = (PyArrayObject*)
@@ -91 +115 @@
   
   M = (row_ptr -> dimensions[0])-1;
+    
+  if (x -> nd == 1) {
+    // Multiplicant is a vector
+  
+    //Allocate space for return vectors y (don't DECREF) 
+    dimensions[0] = M;
+    y = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
+  
+    err = _csr_mv(M,
+                  (double*) data -> data, 
+                  (long*)   colind -> data,
+                  (long*)   row_ptr -> data,
+                  (double*) x -> data,
+                  (double*) y -> data); 
 
+                           
+    if (err != 0) {
+      PyErr_SetString(PyExc_RuntimeError, "matrix vector mult could not be calculated");
+      return NULL;
+    }
+  } else if(x -> nd == 2) {
+  
 
-  //printf(" x.dim[0] = %i  x.dim[1] = %i \n", x -> dimensions[0], x -> dimensions[1]);
+    rows = x -> dimensions[0];     //Number of rows in x        
+    columns = x -> dimensions[1];  //Number of columns in x        
     
-  //Allocate space for return vectors y (don't DECREF) 
-  dimensions[0] = M;
-  y = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
+    //Allocate space for return matrix y (don't DECREF) 
+    dimensions[0] = M;                   //Number of rows in sparse matrix  
+    dimensions[1] = columns;
+    y = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
+    
+    err = _csr_mm(M, columns,
+                  (double*) data -> data, 
+                  (long*)   colind -> data,
+                  (long*)   row_ptr -> data,
+                  (double*) x -> data,
+                  (double*) y -> data); 
+    
+  } else {
+    PyErr_SetString(PyExc_RuntimeError, 
+                    "Allowed dimensions in sparse_ext.c restricted to 1 or 2");
+    return NULL;  
+  }
   
-/*   for (i=0;i<25;i++){ */
-/*     printf(" data[%i]=%i\n",i,(data->data)[i]); */
-/*   } */
-
-/*   for (i=0;i<25;i++){ */
-/*     printf(" colind[%i]=%i\n",i,(colind->data)[i]); */
-/*   } */
-
-  err = _csr_mv(M,
-                (double*) data -> data, 
-                (long*)    colind -> data,
-                (long*)    row_ptr -> data,
-                (double*) x -> data,
-                (double*) y -> data); 
-
-/*   for (i=0;i<25;i++){ */
-/*     printf(" data[%i]=%i\n",i,data->data[i]); */
-/*   } */
-                           
-  if (err != 0) {
-    PyErr_SetString(PyExc_RuntimeError, "matrix vector mult could not be calculated");
-    return NULL;
-  }
                      
   //Release