Index: inundation/ga/storm_surge/pyvolution/sparse.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/sparse.py	(revision 573)
+++ inundation/ga/storm_surge/pyvolution/sparse.py	(revision 586)
@@ -17,5 +17,5 @@
         """
 
-        self.A = {}
+        self.Data = {}
             
         if len(args) == 1:
@@ -32,5 +32,5 @@
                 for j in range(self.N):
                     if A[i, j] != 0.0:
-                        self.A[i, j] = A[i, j]
+                        self.Data[i, j] = A[i, j]
                 
             
@@ -45,10 +45,10 @@
 
     def __repr__(self):
-        return '%d X %d sparse matrix:\n' %(self.M, self.N) + `self.A`
+        return '%d X %d sparse matrix:\n' %(self.M, self.N) + `self.Data`
 
     def __len__(self):
         """Return number of nonzeros of A
         """
-        return len(self.A)
+        return len(self.Data)
 
     def nonzeros(self):
@@ -64,8 +64,8 @@
 
         if x != 0:
-            self.A[key] = float(x)
-        else:
-            if self.A.has_key( key ):            
-                del self.A[key]
+            self.Data[key] = float(x)
+        else:
+            if self.Data.has_key( key ):            
+                del self.Data[key]
 
     def __getitem__(self, key):
@@ -75,6 +75,6 @@
         assert 0 <= j < self.N                
 
-        if self.A.has_key( key ):
-            return self.A[ key ]
+        if self.Data.has_key( key ):
+            return self.Data[ key ]
         else:
             return 0.0
@@ -84,8 +84,8 @@
         new = Sparse(self.M,self.N)
 
-        for key in self.A.keys():
+        for key in self.Data.keys():
             i, j = key
 
-            new[i,j] = self.A[i,j]
+            new[i,j] = self.Data[i,j]
 
         return new
@@ -99,9 +99,10 @@
         for i in range(self.M):
             for j in range(self.N):
-                if self.A.has_key( (i,j) ):                
-                    D[i, j] = self.A[ (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
@@ -129,8 +130,8 @@
 	    
             #Multiply nonzero elements
-            for key in self.A.keys():
+            for key in self.Data.keys():
                 i, j = key
 
-                R[i] += self.A[key]*B[j]
+                R[i] += self.Data[key]*B[j]
         elif len(B.shape) == 2:
 	
@@ -142,8 +143,8 @@
 	        #For each column
 		
-                for key in self.A.keys():
+                for key in self.Data.keys():
                     i, j = key
 
-                    R[i, col] += self.A[key]*B[j, col]
+                    R[i, col] += self.Data[key]*B[j, col]
 	    
 	    
@@ -161,8 +162,8 @@
         
         new = other.copy()
-        for key in self.A.keys():
+        for key in self.Data.keys():
             i, j = key
 
-            new[i,j] += self.A[key]
+            new[i,j] += self.Data[key]
 
         return new
@@ -183,8 +184,8 @@
             new = self.copy()
             #Multiply nonzero elements
-            for key in new.A.keys():
+            for key in new.Data.keys():
                 i, j = key
 
-                new.A[key] = other*new.A[key]
+                new.Data[key] = other*new.Data[key]
 
         return new
@@ -213,8 +214,8 @@
 
             #Multiply nonzero elements
-            for key in self.A.keys():
+            for key in self.Data.keys():
                 i, j = key
 
-                R[j] += self.A[key]*B[i]
+                R[j] += self.Data[key]*B[i]
 
         else:
@@ -223,7 +224,119 @@
         return R
 
-
-    
-    
+class Sparse_CSR:
+
+    def __init__(self, A):
+        """Create sparse matrix in csr format.
+
+        Sparse_CSR(A) #creates csr sparse matrix from sparse matrix
+        """
+
+        from Numeric import array, Float, Int
+
+        if isinstance(A,Sparse):
+
+            from Numeric import zeros
+            keys = A.Data.keys()
+            keys.sort()
+            nnz = len(keys)
+            data    = zeros ( (nnz,), Float)
+            colind  = zeros ( (nnz,), Int)
+            row_ptr = zeros ( (A.M+1,), Int)
+            current_row = -1
+            k = 0
+            for key in keys:
+                ikey0 = int(key[0])
+                ikey1 = int(key[1])
+                if ikey0 != current_row:
+                    current_row = ikey0
+                    row_ptr[ikey0] = k
+                data[k] = A.Data[key]
+                colind[k] = ikey1
+                k += 1
+            for row in range(current_row+1, A.M+1):
+                row_ptr[row] = nnz
+            #row_ptr[-1] = nnz
+        
+            self.data    = data
+            self.colind  = colind
+            self.row_ptr = row_ptr
+            self.M       = A.M
+            self.N       = A.N
+        else:
+            raise ValueError, "Sparse_CSR(A) expects A == Sparse Matrix"
+            
+    def __repr__(self):
+        return '%d X %d sparse matrix:\n' %(self.M, self.N) + `self.data`
+
+    def __len__(self):
+        """Return number of nonzeros of A
+        """
+        return self.row_ptr[-1]
+
+    def nonzeros(self):
+        """Return number of nonzeros of A
+        """        
+        return len(self)
+
+    def todense(self):
+        from Numeric import zeros, Float
+
+        D = zeros( (self.M, self.N), Float)
+        
+        for i in range(self.M):
+            for ckey in range(self.row_ptr[i],self.row_ptr[i+1]):
+                j = self.colind[ckey]
+                D[i, j] = self.data[ckey]
+        return D
+
+    def __mul__(self, other):
+        """Multiply this matrix onto 'other' which can either be
+        a Numeric vector, a Numeric matrix or another sparse matrix.
+        """
+
+        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) == 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
+    
+
 if __name__ == '__main__':
 
Index: inundation/ga/storm_surge/pyvolution/sparse_ext.c
===================================================================
--- inundation/ga/storm_surge/pyvolution/sparse_ext.c	(revision 586)
+++ inundation/ga/storm_surge/pyvolution/sparse_ext.c	(revision 586)
@@ -0,0 +1,583 @@
+// Python - C extension for sparse module.
+//
+// To compile (Python2.3):
+//  gcc -c util_ext.c -I/usr/include/python2.3 -o util_ext.o -Wall -O 
+//  gcc -shared util_ext.o  -o util_ext.so	
+//
+// See the module sparse.py
+//
+//
+// Steve Roberts, ANU 2004
+	
+#include "Python.h"
+#include "Numeric/arrayobject.h"
+#include "math.h"
+
+//Shared code snippets
+#include "util_ext.h"
+#include "quantity_ext.h"
+
+
+
+int _compute_csr_mult(int M,
+		      double* data, 
+		      int* colind,
+		      int* row_ptr,
+		      double* R) {
+
+			
+  int i, j, ckey;
+
+  for (i=0; i<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]              
+  
+  for (k=0; k<N; k++) {
+    index3 = 3*k;
+    
+    if (number_of_boundaries[k] < 2) {
+      //Two or three true neighbours
+
+      //Get indices of neighbours (or self when used as surrogate)	    
+      //k0, k1, k2 = surrogate_neighbours[k,:]
+      
+      k0 = surrogate_neighbours[index3 + 0];
+      k1 = surrogate_neighbours[index3 + 1];
+      k2 = surrogate_neighbours[index3 + 2];           
+      if (k0 == k1 || k1 == k2) return -1;      
+
+      //Get data
+      q0 = centroid_values[k0];
+      q1 = centroid_values[k1];
+      q2 = centroid_values[k2];                    
+
+      x0 = centroids[k0*2]; y0 = centroids[k0*2+1]; 
+      x1 = centroids[k1*2]; y1 = centroids[k1*2+1]; 
+      x2 = centroids[k2*2]; y2 = centroids[k2*2+1];             
+
+      //Gradient
+      _gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2, &a[k], &b[k]);
+      
+    } else if (number_of_boundaries[k] == 2) {
+      //One true neighbour
+
+      //#Get index of the one neighbour
+      i=0; k0 = k;
+      while (i<3 && k0==k) {
+	k0 = surrogate_neighbours[index3 + i];
+	i++;
+      }
+      if (k0 == k) return -1;
+      
+      k1 = k; //self
+
+      //Get data
+      q0 = centroid_values[k0];
+      q1 = centroid_values[k1];
+	    
+      x0 = centroids[k0*2]; y0 = centroids[k0*2+1]; 
+      x1 = centroids[k1*2]; y1 = centroids[k1*2+1]; 
+
+      //Gradient
+      det = x0*y1 - x1*y0;
+      if (det != 0.0) {
+	a[k] = (y1*q0 - y0*q1)/det;
+	b[k] = (x0*q1 - x1*q0)/det;
+      }
+    } 
+    //    else:
+    //        #No true neighbours - 	
+    //        #Fall back to first order scheme
+    //        pass
+  }
+  return 0;
+}            
+
+
+int _extrapolate(int N,
+		 double* centroids, 
+		 double* centroid_values, 		 
+		 double* vertex_coordinates,
+		 double* vertex_values,
+		 double* a,
+		 double* b) {
+  
+  int k, k2, k3, k6;
+  double x, y, x0, y0, x1, y1, x2, y2;
+		     
+  for (k=0; k<N; k++) {
+    k6 = 6*k;
+    k3 = 3*k;    
+    k2 = 2*k;
+        
+    //Centroid coordinates            
+    x = centroids[k2]; y = centroids[k2+1]; 
+
+    //vertex coordinates
+    //x0, y0, x1, y1, x2, y2 = X[k,:]
+    x0 = vertex_coordinates[k6 + 0];
+    y0 = vertex_coordinates[k6 + 1];    
+    x1 = vertex_coordinates[k6 + 2];
+    y1 = vertex_coordinates[k6 + 3];        
+    x2 = vertex_coordinates[k6 + 4];
+    y2 = vertex_coordinates[k6 + 5];           
+
+    //Extrapolate
+    vertex_values[k3+0] = centroid_values[k] + a[k]*(x0-x) + b[k]*(y0-y);
+    vertex_values[k3+1] = centroid_values[k] + a[k]*(x1-x) + b[k]*(y1-y);
+    vertex_values[k3+2] = centroid_values[k] + a[k]*(x2-x) + b[k]*(y2-y);
+    
+  }
+  return 0;
+}
+
+
+
+
+int _interpolate(int N,
+		 double* vertex_values,
+		 double* edge_values) {
+  
+  int k, k3;
+  double q0, q1, q2;
+
+		     
+  for (k=0; k<N; k++) {
+    k3 = 3*k;    
+    
+    q0 = vertex_values[k3 + 0];
+    q1 = vertex_values[k3 + 1];
+    q2 = vertex_values[k3 + 2];
+            
+    //printf("%f, %f, %f\n", q0, q1, q2);
+    edge_values[k3 + 0] = 0.5*(q1+q2);
+    edge_values[k3 + 1] = 0.5*(q0+q2);  
+    edge_values[k3 + 2] = 0.5*(q0+q1);
+  }  
+  return 0;
+}
+
+int _update(int N, 
+	    double timestep,
+	    double* centroid_values,
+	    double* explicit_update,
+	    double* semi_implicit_update) {
+  //Update centroid values based on values stored in
+  //explicit_update and semi_implicit_update as well as given timestep
+  
+  
+  int k;
+  double denominator, x;
+
+  
+  //Divide semi_implicit update by conserved quantity  
+  for (k=0; k<N; k++) {  
+    x = centroid_values[k];
+    if (x == 0.0) {
+      //FIXME: Is this right
+      semi_implicit_update[k] = 0.0;
+    } else {
+      semi_implicit_update[k] /= x;
+    }    
+  }              
+  
+  
+  //Explicit updates
+  for (k=0; k<N; k++) {
+    centroid_values[k] += timestep*explicit_update[k];
+  }
+            
+  //Semi implicit updates
+  for (k=0; k<N; k++) {  
+    denominator = 1.0 - timestep*semi_implicit_update[k];  
+    
+    if (denominator == 0.0) {
+      return -1;
+    } else {
+      //Update conserved_quantities from semi implicit updates
+      centroid_values[k] /= denominator;
+    }
+  }
+  return 0;
+}
+		     
+  
+/////////////////////////////////////////////////
+// Gateways to Python
+PyObject *update(PyObject *self, PyObject *args) {
+  
+  PyObject *quantity;
+  PyArrayObject *centroid_values, *explicit_update, *semi_implicit_update; 
+  
+  double timestep;  
+  int N, err;
+  
+  // Convert Python arguments to C  
+  if (!PyArg_ParseTuple(args, "Od", &quantity, &timestep)) 
+    return NULL;
+
+  centroid_values = (PyArrayObject*) 
+    PyObject_GetAttrString(quantity, "centroid_values");           
+  if (!centroid_values) return NULL;                
+
+  explicit_update = (PyArrayObject*) 
+    PyObject_GetAttrString(quantity, "explicit_update");   
+  if (!explicit_update) return NULL;          
+  
+  semi_implicit_update = (PyArrayObject*) 
+    PyObject_GetAttrString(quantity, "semi_implicit_update");   
+  if (!semi_implicit_update) return NULL;            
+
+  N = centroid_values -> dimensions[0];  
+  
+  
+  err = _update(N, timestep,
+		(double*) centroid_values -> data,
+		(double*) explicit_update -> data,
+		(double*) semi_implicit_update -> data);		
+		
+		     
+  if (err != 0) {
+    PyErr_SetString(PyExc_RuntimeError, 
+		    "Zero division in semi implicit update - call Stephen :)");
+    return NULL;
+  }		     
+  
+  //Release and return		     
+  Py_DECREF(centroid_values);    
+  Py_DECREF(explicit_update); 
+  Py_DECREF(semi_implicit_update);   
+
+  return Py_BuildValue("");
+}
+
+
+PyObject *interpolate_from_vertices_to_edges(PyObject *self, PyObject *args) {
+  
+  PyObject *quantity;
+  PyArrayObject *vertex_values, *edge_values;
+    
+  int N, err;
+  
+  // Convert Python arguments to C  
+  if (!PyArg_ParseTuple(args, "O", &quantity)) 
+    return NULL;
+
+  vertex_values = (PyArrayObject*) 
+    PyObject_GetAttrString(quantity, "vertex_values");           
+  if (!vertex_values) return NULL;                
+
+  edge_values = (PyArrayObject*) 
+    PyObject_GetAttrString(quantity, "edge_values");   
+  if (!edge_values) return NULL;          
+
+  N = vertex_values -> dimensions[0];  
+  
+  err = _interpolate(N,
+		     (double*) vertex_values -> data,
+		     (double*) edge_values -> data);
+		     
+  if (err != 0) {
+    PyErr_SetString(PyExc_RuntimeError, "Interpolate could not be computed");
+    return NULL;
+  }		     
+  
+  //Release and return		     
+  Py_DECREF(vertex_values);    
+  Py_DECREF(edge_values); 
+
+  return Py_BuildValue("");
+}
+
+
+PyObject *compute_gradients(PyObject *self, PyObject *args) {
+  
+  PyObject *quantity, *domain, *R;
+  PyArrayObject 
+    *centroids,            //Coordinates at centroids
+    *centroid_values,      //Values at centroids    
+    *number_of_boundaries, //Number of boundaries for each triangle
+    *surrogate_neighbours, //True neighbours or - if one missing - self
+    *a, *b;                //Return values
+   
+  int dimensions[1], N, err;
+  
+  // Convert Python arguments to C  
+  if (!PyArg_ParseTuple(args, "O", &quantity)) 
+    return NULL;
+
+  domain = PyObject_GetAttrString(quantity, "domain");     
+  if (!domain) 
+    return NULL;  
+
+  //Get pertinent variables
+  centroids = (PyArrayObject*) 
+    PyObject_GetAttrString(domain, "centroid_coordinates"); 
+  if (!centroids) return NULL;    
+  
+  centroid_values = (PyArrayObject*) 
+    PyObject_GetAttrString(quantity, "centroid_values");   
+  if (!centroid_values) return NULL;        
+  
+  surrogate_neighbours = (PyArrayObject*) 
+    PyObject_GetAttrString(domain, "surrogate_neighbours");
+  if (!surrogate_neighbours) return NULL;        
+  
+  number_of_boundaries = (PyArrayObject*) 
+    PyObject_GetAttrString(domain, "number_of_boundaries");           
+  if (!number_of_boundaries) return NULL;            
+  
+  N = centroid_values -> dimensions[0];
+    
+  //Release
+  Py_DECREF(domain);      
+          
+  //Allocate space for return vectors a and b (don't DECREF) 
+  dimensions[0] = N;
+  a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
+  b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);  
+        
+
+  
+  err = _compute_gradients(N,
+			   (double*) centroids -> data, 
+			   (double*) centroid_values -> data,
+			   (int*) number_of_boundaries -> data,
+			   (int*) surrogate_neighbours -> data,
+			   (double*) a -> data,
+			   (double*) b -> data);     
+			   
+  if (err != 0) {
+    PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed");
+    return NULL;
+  }
+  		     
+  //Release		     
+  Py_DECREF(centroids);    
+  Py_DECREF(centroid_values);    
+  Py_DECREF(number_of_boundaries); 
+  Py_DECREF(surrogate_neighbours);           
+         
+  //Build result, release and return
+  R = Py_BuildValue("OO", PyArray_Return(a), PyArray_Return(b)); 
+  Py_DECREF(a);      
+  Py_DECREF(b);        
+  return R;
+}
+
+
+
+PyObject *extrapolate_second_order(PyObject *self, PyObject *args) {
+  
+  PyObject *quantity, *domain;
+  PyArrayObject 
+    *centroids,            //Coordinates at centroids
+    *centroid_values,      //Values at centroids    
+    *vertex_coordinates,   //Coordinates at vertices            
+    *vertex_values,        //Values at vertices        
+    *number_of_boundaries, //Number of boundaries for each triangle
+    *surrogate_neighbours, //True neighbours or - if one missing - self
+    *a, *b;                //Gradients
+   
+  //int N, err;
+  int dimensions[1], N, err;  
+  //double *a, *b;  //Gradients
+  
+  // Convert Python arguments to C  
+  if (!PyArg_ParseTuple(args, "O", &quantity)) 
+    return NULL;
+
+  domain = PyObject_GetAttrString(quantity, "domain");     
+  if (!domain) 
+    return NULL;  
+
+  //Get pertinent variables
+  centroids = (PyArrayObject*) 
+    PyObject_GetAttrString(domain, "centroid_coordinates"); 
+  if (!centroids) return NULL;    
+  
+  centroid_values = (PyArrayObject*) 
+    PyObject_GetAttrString(quantity, "centroid_values");   
+  if (!centroid_values) return NULL;        
+  
+  surrogate_neighbours = (PyArrayObject*) 
+    PyObject_GetAttrString(domain, "surrogate_neighbours");
+  if (!surrogate_neighbours) return NULL;        
+  
+  number_of_boundaries = (PyArrayObject*) 
+    PyObject_GetAttrString(domain, "number_of_boundaries");           
+  if (!number_of_boundaries) return NULL;            
+  
+  vertex_coordinates = (PyArrayObject*) 
+    PyObject_GetAttrString(domain, "vertex_coordinates");           
+  if (!vertex_coordinates) return NULL;              
+  
+  vertex_values = (PyArrayObject*) 
+    PyObject_GetAttrString(quantity, "vertex_values");           
+  if (!vertex_values) return NULL;                
+
+  
+  /*
+  printf("In extrapolate C routine\n");  
+  printf("d0=%d, d1=%d\n", 
+	 vertex_values -> dimensions[0],
+	 vertex_values -> dimensions[1]);
+  */
+  
+  vertex_values = (PyArrayObject*) 
+    PyObject_GetAttrString(quantity, "vertex_values");           
+  if (!vertex_values) return NULL;                
+  
+  N = centroid_values -> dimensions[0];
+    
+  //Release
+  Py_DECREF(domain);      
+          
+  //Allocate space for return vectors a and b (don't DECREF) 
+  dimensions[0] = N;
+  a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
+  b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);  
+
+  //FIXME: Odd that I couldn't use normal arrays	  
+  //Allocate space for return vectors a and b (don't DECREF) 
+  //a = (double*) malloc(N * sizeof(double));
+  //if (!a) return NULL;  
+  //b = (double*) malloc(N * sizeof(double));  
+  //if (!b) return NULL;
+
+  
+  err = _compute_gradients(N,
+			   (double*) centroids -> data, 
+			   (double*) centroid_values -> data,
+			   (int*) number_of_boundaries -> data,
+			   (int*) surrogate_neighbours -> data,
+			   (double*) a -> data,
+			   (double*) b -> data);
+			   
+  if (err != 0) {
+    PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed");
+    return NULL;
+  }
+  
+  err = _extrapolate(N, 
+		     (double*) centroids -> data, 
+		     (double*) centroid_values -> data,		     
+		     (double*) vertex_coordinates -> data,
+		     (double*) vertex_values -> data,
+		     (double*) a -> data,
+		     (double*) b -> data);     		     
+  //a, b);
+ 		     
+		     
+  if (err != 0) {
+    PyErr_SetString(PyExc_RuntimeError, 
+		    "Internal function _extrapolate failed");
+    return NULL;
+  }		     
+		
+  
+
+  //Release  
+  Py_DECREF(centroids);    
+  Py_DECREF(centroid_values);    
+  Py_DECREF(number_of_boundaries); 
+  Py_DECREF(surrogate_neighbours);
+  Py_DECREF(vertex_coordinates);                   
+  Py_DECREF(vertex_values);                 
+  Py_DECREF(a);                   
+  Py_DECREF(b);                     
+  //free(a);      
+  //free(b);        
+  
+  return Py_BuildValue("");    
+}    
+
+
+
+PyObject *limit(PyObject *self, PyObject *args) {
+  
+  PyObject *quantity, *domain, *Tmp;
+  PyArrayObject 
+    *qv, //Conserved quantities at vertices 
+    *qc, //Conserved quantities at centroids
+    *neighbours;
+    
+  int k, i, n, N, k3;
+  double beta; //Safety factor
+  double *qmin, *qmax, qn;
+  
+  // Convert Python arguments to C  
+  if (!PyArg_ParseTuple(args, "O", &quantity)) 
+    return NULL;
+
+  domain = PyObject_GetAttrString(quantity, "domain");     
+  if (!domain) 
+    return NULL;  
+    
+  neighbours = (PyArrayObject*) PyObject_GetAttrString(domain, "neighbours");
+
+  //Get safety factor beta
+  Tmp = PyObject_GetAttrString(domain, "beta"); 
+  if (!Tmp) 
+    return NULL;
+      
+  beta = PyFloat_AsDouble(Tmp);  
+  
+  Py_DECREF(Tmp);    
+  Py_DECREF(domain);      
+
+  qc = (PyArrayObject*) PyObject_GetAttrString(quantity, "centroid_values"); 
+  qv = (PyArrayObject*) PyObject_GetAttrString(quantity, "vertex_values");   
+  N = qc -> dimensions[0];
+    
+  //Find min and max of this and neighbour's centroid values
+  qmin = malloc(N * sizeof(double));
+  qmax = malloc(N * sizeof(double));  
+  for (k=0; k<N; k++) {
+    k3=k*3;
+    
+    qmin[k] = ((double*) qc -> data)[k];
+    qmax[k] = qmin[k];
+    
+    for (i=0; i<3; i++) {
+      n = ((int*) neighbours -> data)[k3+i];
+      if (n >= 0) {
+	qn = ((double*) qc -> data)[n]; //Neighbour's centroid value
+	  
+	qmin[k] = min(qmin[k], qn);
+	qmax[k] = max(qmax[k], qn);
+      }
+    }
+  }
+  
+  // Call underlying routine
+  _limit(N, beta, (double*) qc -> data, (double*) qv -> data, qmin, qmax);
+	 
+  free(qmin);
+  free(qmax);  
+  return Py_BuildValue("");  
+}
+
+
+
+// Method table for python module
+static struct PyMethodDef MethodTable[] = {
+  {"limit", limit, METH_VARARGS, "Print out"},    
+  {"update", update, METH_VARARGS, "Print out"},      
+  {"compute_gradients", compute_gradients, METH_VARARGS, "Print out"},      
+  {"extrapolate_second_order", extrapolate_second_order, 
+   METH_VARARGS, "Print out"},        
+  {"interpolate_from_vertices_to_edges", 
+   interpolate_from_vertices_to_edges,
+   METH_VARARGS, "Print out"},           
+  {NULL, NULL, 0, NULL}   /* sentinel */
+};
+
+// Module initialisation   
+void initquantity_ext(void){
+  Py_InitModule("quantity_ext", MethodTable);
+  
+  import_array();     //Necessary for handling of NumPY structures  
+}
+
Index: inundation/ga/storm_surge/pyvolution/test_cg_solve.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/test_cg_solve.py	(revision 573)
+++ inundation/ga/storm_surge/pyvolution/test_cg_solve.py	(revision 586)
@@ -6,5 +6,5 @@
 from Numeric import dot, allclose, array, transpose, arange, ones, Float
 from cg_solve import *
-from sparse import Sparse
+from sparse import Sparse, Sparse_CSR
 
 
@@ -78,4 +78,36 @@
         assert allclose(x,xe)
 
+    def test_solve_large_2d_csr_matrix(self):
+        """Standard 2d laplacian"""
+        
+        n = 50
+        m = 50
+
+        A = Sparse(m*n, m*n)
+
+        for i in arange(0,n):
+            for j in arange(0,m):
+                I = j+m*i
+                A[I,I] = 4.0
+                if i > 0  :
+                    A[I,I-m] = -1.0
+                if i < n-1 :
+                    A[I,I+m] = -1.0
+                if j > 0  :
+                    A[I,I-1] = -1.0
+                if j < m-1 :
+                    A[I,I+1] = -1.0
+                
+        xe = ones( (n*m,), Float)
+
+        # Convert to csr format
+        #print 'start covert'
+        A = Sparse_CSR(A)
+        #print 'finish covert'
+        b = A*xe
+        x = conjugate_gradient(A,b,b,iprint=10)
+
+        assert allclose(x,xe)
+
 
     def test_solve_large_2d_with_default_guess(self):
Index: inundation/ga/storm_surge/pyvolution/test_sparse.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/test_sparse.py	(revision 573)
+++ inundation/ga/storm_surge/pyvolution/test_sparse.py	(revision 586)
@@ -164,5 +164,33 @@
         C = A+B
 
-        assert allclose(C.todense(), [[12,0,0], [2,8,8], [0,0,4]])        
+        assert allclose(C.todense(), [[12,0,0], [2,8,8], [0,0,4]])
+
+    def test_sparse_tocsr(self):
+
+        A = Sparse(4,3)
+
+        A[0,0] = 3
+        A[1,1] = 2
+        A[1,2] = 2
+        A[2,2] = 1
+        A[0,2] = 4
+        A[2,0] = 5
+
+        print ' ' 
+        print A.todense()
+
+        B = Sparse_CSR(A)
+
+        print B.todense()
+
+        C = [1, 2, 3]
+
+        assert allclose(B*C, [15.0, 10.0 ,8.0, 0.0])
+
+        C2 = [[1,2],[2,4],[3,6]]
+
+        print B*C2
+
+        assert allclose(B*C2, [[15.0, 30.0],[10.0, 20.0],[8.0, 16.0],[0.0, 0.0]])
         
 
