Index: /inundation/ga/storm_surge/pyvolution/sparse.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/sparse.py	(revision 604)
+++ /inundation/ga/storm_surge/pyvolution/sparse.py	(revision 605)
@@ -306,10 +306,7 @@
 
 
-"""
-Setup for C extensions
-"""
-
-
 def csr_mv(self, B):
+    """Python version of sparse (CSR) matrix multiplication
+    """
 
     from Numeric import zeros, Float
@@ -353,9 +350,10 @@
 
 
+
+#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__':
Index: /inundation/ga/storm_surge/pyvolution/sparse_ext.c
===================================================================
--- /inundation/ga/storm_surge/pyvolution/sparse_ext.c	(revision 604)
+++ /inundation/ga/storm_surge/pyvolution/sparse_ext.c	(revision 605)
@@ -15,4 +15,5 @@
 #include "stdio.h"
 
+//Matrix-vector routine
 int _csr_mv(int M,
 	    double* data, 
@@ -24,5 +25,4 @@
   long i, j, ckey;
 
-  //printf("test");
   for (i=0; i<M; i++ ) 
     for (ckey=row_ptr[i]; ckey<row_ptr[i+1]; ckey++) {
@@ -31,8 +31,33 @@
     }              
   
-  //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,5 +68,5 @@
   
   PyObject *csr_sparse, // input sparse matrix
-    *xin, *R;                 // output wrapped vector
+    *xin, *R;           // output wrapped vector
   
   PyArrayObject 
@@ -53,5 +78,5 @@
 
   
-  int dimensions[1], M, err;
+  int dimensions[1], M, err, columns, rows;
   
   // Convert Python arguments to C  
@@ -73,5 +98,5 @@
 /*   printf("x.data[5] = %g\n",((double*) x->data)[5]); */
 
-  
+ 
 
 
@@ -80,5 +105,4 @@
   if (!data) 
     return NULL;  
-
 
   colind = (PyArrayObject*)
@@ -91,35 +115,48 @@
   
   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		     
