Index: inundation/ga/storm_surge/pyvolution/least_squares.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/least_squares.py	(revision 484)
+++ inundation/ga/storm_surge/pyvolution/least_squares.py	(revision 485)
@@ -19,4 +19,8 @@
 
 
+#FIXME (Ole): Currently datapoints outside the triangular mesh are ignored.
+#             Is there a clean way of inlcuding them?
+
+
 import exceptions
 class ShapeError(exceptions.Exception): pass
@@ -30,5 +34,4 @@
 try:
     from util import gradient
-
 except ImportError, e: 
     #FIXME reduce the dependency of modules in pyvolution
@@ -50,5 +53,4 @@
 DEFAULT_ALPHA = 0.001
     
-
 def fit_to_mesh_file(mesh_file, point_file, mesh_output_file, alpha=DEFAULT_ALPHA):
     """
@@ -67,12 +69,11 @@
     old_title_list = mesh_dict['generatedpointattributetitlelist']
 
-    
+   
     # load in the .xya file
-    
     point_dict = load_xya_file(point_file)
     point_coordinates = point_dict['pointlist']
     point_attributes = point_dict['pointattributelist']
     title_string = point_dict['title']
-    title_list = title_string.split(',') #iffy! Hard coding title delimiter  
+    title_list = title_string.split(',') #FIXME iffy! Hard coding title delimiter  
     for i in range(len(title_list)):
         title_list[i] = title_list[i].strip() 
@@ -218,9 +219,6 @@
         self.AtA = Sparse(m,m)
 
-        #Build quad tree of vertices
-        #print 'Building quad tree'
+        #Build quad tree of vertices (FIXME: Is this the right spot for that?)
         root = build_quadtree(self.mesh)
-        #print 'Quad tree built'
-        #root.show()
 
         #Compute matrix elements
@@ -235,15 +233,17 @@
             candidate_vertices = root.search(x[0], x[1])
 
-            #Find triangle containing x
+            #Find triangle containing x:
             element_found = False            
+	    
+	    #For all vertices in same cell as point x
             for v in candidate_vertices:
-                for triangle_id, vertex_id in self.mesh.vertexlist[v]:
-
-                    k = triangle_id
+	    
+                #for each triangle id (k) which has v as a vertex 
+		for k, _ in self.mesh.vertexlist[v]:
                     
                     #Get the three vertex_points of candidate triangle
                     xi0 = self.mesh.get_vertex_coordinate(k, 0)
                     xi1 = self.mesh.get_vertex_coordinate(k, 1)
-                    xi2 = self.mesh.get_vertex_coordinate(k, 2)                 
+                    xi2 = self.mesh.get_vertex_coordinate(k, 2)     
 
                     #Get the three normals
@@ -269,5 +269,7 @@
                     #Don't look for any other triangle
                     break
-
+		    
+
+	    #Update interpolation matrix A if necessary	    
             if element_found is True:        
                 #Assign values to matrix A
@@ -291,5 +293,6 @@
             else:
                 pass
-                #Ok if there is no triangle for datapoint (as in brute force version)
+                #Ok if there is no triangle for datapoint 
+		#(as in brute force version)
                 #raise 'Could not find triangle for point', x
 
@@ -301,6 +304,6 @@
         m is the number of basis functions phi_k (one per vertex)
 
-
-        This is the brute force which is too slow for large problems, but good for testing
+        This is the brute force which is too slow for large problems, 
+	but could be used for testing
         """
 
@@ -413,8 +416,6 @@
         #"element stiffness matrices:
 
-        
         m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex)
 
-        #self.D = zeros((m,m), Float)
         self.D = Sparse(m,m)
 
@@ -463,5 +464,4 @@
             self.D[v0,v2] += e20             
 
-        #self.D = (self.D).tocsc()
             
     def fit(self, z):
@@ -488,6 +488,4 @@
         Atz = self.A.trans_mult(z)
 
-
-        #print 'fit: Atz',Atz
         
         #Check sanity
@@ -500,17 +498,14 @@
 
 
-        #Solve and return
-        #return solve_linear_equations(self.get_B(), Atz)
-
-        # caused errors
-        #return sparse.solve(self.B,Atz)
-
         return conjugate_gradient(self.B, Atz, Atz,imax=2*len(Atz) )
         #FIXME: Should we store the result here for later use? (ON)        
-    
+
+	    
     def fit_points(self, z):
-        """
-        Like fit, but more robust when each point has two or more attributes
-        """
+        """Like fit, but more robust when each point has two or more attributes
+	FIXME(Ole): The name fit_points doesn't carry any meaning for me.
+	How about something like fit_multiple or fit_columns?
+        """
+	
         try:
             return self.fit(z)
@@ -525,9 +520,9 @@
             n = z.shape[1]               #Number of data points         
 
-            f = zeros((m,n), Float)
-            #f = sparse.dok_matrix() # even though it wont be sparse?
+            f = zeros((m,n), Float) #Resulting columns
             
             for i in range(z.shape[1]):
                 f[:,i] = self.fit(z[:,i])
+		
             return f
             
@@ -548,29 +543,7 @@
         """
 
-        try:
-            return self.A * f
-        except ValueError:
-            # We are here, so it probalby means that f is 2 dimensional
-            # and so will do each column separately due to problems in
-            # sparse matrix, 2d array multiplication
-            (N , M) = self.A.shape
-            f = array(f).astype(Float)
-            (m , n) = f.shape
-            #print "m",m
-            #print "M",M
-            if m != M :
-                #print "!!self.A",self.A
-                #print "!!self.A.todense() ",self.A.todense() 
-                #print "self.get_A()",self.get_A()
-                #return dot(self.get_A(),f)
-                raise VectorShapeError, 'Mismatch between A and f dimensions'
-            
-            y = zeros( (N,n), Float)
-
-            for i in range(n):
-                y[:,i] = self.A * f[:,i]
-
-            #print "!!self.A.todense() ",self.A.todense() 
-            return y
+        return self.A * f
+	
+	    
      
     #FIXME: We will need a method 'evaluate(self):' that will interpolate
@@ -597,5 +570,6 @@
     """
     import os, sys
-    usage = "usage: %s mesh_input.tsh point.xya  mesh_output.tsh alpha" % os.path.basename(sys.argv[0])
+    usage = "usage: %s mesh_input.tsh point.xya mesh_output.tsh alpha"\
+            %os.path.basename(sys.argv[0])
 
     if len(sys.argv) < 4:
Index: inundation/ga/storm_surge/pyvolution/sparse.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/sparse.py	(revision 484)
+++ inundation/ga/storm_surge/pyvolution/sparse.py	(revision 485)
@@ -116,15 +116,16 @@
             print 'FIXME: Only Numeric types implemented so far'
 
-
         #Assume numeric types from now on
-        R = zeros((self.M,), Float) #Result
-
+	
         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 key in self.A.keys():
@@ -133,5 +134,18 @@
                 R[i] += self.A[key]*B[j]
         elif len(B.shape) == 2:
-            raise ValueError, 'Numeric matrix not yet implemented'            
+	
+            
+            R = zeros((self.M, B.shape[1]), Float) #Result matrix
+
+            #Multiply nonzero elements
+	    for col in range(R.shape[1]):
+	        #For each column
+		
+                for key in self.A.keys():
+                    i, j = key
+
+                    R[i, col] += self.A[key]*B[j, col]
+	    
+	    
         else:
             raise ValueError, 'Dimension too high: d=%d' %len(B.shape)
Index: inundation/ga/storm_surge/pyvolution/test_least_squares.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/test_least_squares.py	(revision 484)
+++ inundation/ga/storm_surge/pyvolution/test_least_squares.py	(revision 485)
@@ -168,5 +168,5 @@
         #print "answer\n",answer
         
-        assert allclose(f, answer,atol=1e-7)
+        assert allclose(f, answer, atol=1e-7)
         
         
@@ -200,5 +200,5 @@
 
         vertices = [a, b, c, d,e,f]
-        triangles = [[0,1,3],[1,0,2],[0,4,5], [0,5,2]] #abd bac aef afc
+        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
 
         point_coords = [[-2.0, 2.0],
@@ -238,9 +238,9 @@
         d1 = [1.0, 1.0]
         d2 = [1.0, 3.0]
-        d3 = [3.0,1.0]
-        z1 = [2,4]
-        z2 = [4,8]
-        z3 = [4,8]
-        data_coords = [ d1, d2, d3] 
+        d3 = [3.0, 1.0]
+        z1 = [2, 4]
+        z2 = [4, 8]
+        z3 = [4, 8]
+        data_coords = [d1, d2, d3] 
         
         interp = Interpolation(points, triangles, data_coords, alpha=0.0)
@@ -260,5 +260,5 @@
         d0 = [1.0, 1.0]
         d1 = [1.0, 2.0]
-        d2 = [3.0,1.0]
+        d2 = [3.0, 1.0]
         point_coords = [ d0, d1, d2]
         
@@ -275,25 +275,25 @@
         a = [-1.0, 0.0]
         b = [3.0, 4.0]
-        c = [4.0,1.0]
+        c = [4.0, 1.0]
         d = [-3.0, 2.0] #3
-        e = [-1.0,-2.0]
+        e = [-1.0, -2.0]
         f = [1.0, -2.0] #5
 
         vertices = [a, b, c, d,e,f]
-        triangles = [[0,1,3],[1,0,2],[0,4,5], [0,5,2]] #abd bac aef afc
+        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
 
          
         point_coords = [[-2.0, 2.0],
                         [-1.0, 1.0],
-                        [0.0,2.0],
+                        [0.0, 2.0],
                         [1.0, 1.0],
                         [2.0, 1.0],
-                        [0.0,0.0],
+                        [0.0, 0.0],
                         [1.0, 0.0],
                         [0.0, -1.0],
-                        [-0.2,-0.5],
+                        [-0.2, -0.5],
                         [-0.9, -1.5],
                         [0.5, -1.9],
-                        [3.0,1.0]]
+                        [3.0, 1.0]]
         
         interp = Interpolation(vertices, triangles, point_coords)
@@ -335,25 +335,25 @@
         a = [-1.0, 0.0]
         b = [3.0, 4.0]
-        c = [4.0,1.0]
+        c = [4.0, 1.0]
         d = [-3.0, 2.0] #3
-        e = [-1.0,-2.0]
+        e = [-1.0, -2.0]
         f = [1.0, -2.0] #5
 
         vertices = [a, b, c, d,e,f]
-        triangles = [[0,1,3],[1,0,2],[0,4,5], [0,5,2]] #abd bac aef afc
+        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
 
          
         point_coords = [[-2.0, 2.0],
                         [-1.0, 1.0],
-                        [0.0,2.0],
+                        [0.0, 2.0],
                         [1.0, 1.0],
                         [2.0, 1.0],
-                        [0.0,0.0],
+                        [0.0, 0.0],
                         [1.0, 0.0],
                         [0.0, -1.0],
-                        [-0.2,-0.5],
+                        [-0.2, -0.5],
                         [-0.9, -1.5],
                         [0.5, -1.9],
-                        [3.0,1.0]]
+                        [3.0, 1.0]]
         
         interp = Interpolation(vertices, triangles, point_coords)
@@ -384,9 +384,9 @@
         d1 = [1.0, 1.0]
         d2 = [1.0, 3.0]
-        d3 = [3.0,1.0]
-        z1 = [2,4]
-        z2 = [4,8]
-        z3 = [4,8]
-        data_coords = [ d1, d2, d3]
+        d3 = [3.0, 1.0]
+        z1 = [2, 4]
+        z2 = [4, 8]
+        z3 = [4, 8]
+        data_coords = [d1, d2, d3]
         z = [z1, z2, z3]
         
Index: inundation/ga/storm_surge/pyvolution/test_sparse.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/test_sparse.py	(revision 484)
+++ inundation/ga/storm_surge/pyvolution/test_sparse.py	(revision 485)
@@ -3,5 +3,4 @@
 import unittest
 from math import sqrt
-
 
 from sparse import *
@@ -15,5 +14,4 @@
     def tearDown(self):
         pass
-
 
     def test_init1(Self):
@@ -72,5 +70,5 @@
         assert allclose(A.todense(), [[0,0,0], [0,0,0], [0,0,0]])        
 
-    def test_sparse_multiplication(self):
+    def test_sparse_multiplication_vector(self):
         A = Sparse(3,3)
         
@@ -83,9 +81,6 @@
         v = [2,3,4]
 
-        
         u = A*v
         assert allclose(u, [6,14,4])
-
-
 
         #Right hand side column
@@ -97,4 +92,20 @@
         u = A*v[:,1]
         assert allclose(u, [12,16,4])        
+
+
+    def test_sparse_multiplication_matrix(self):
+        A = Sparse(3,3)
+        
+        A[0,0] = 3
+        A[1,1] = 2
+        A[1,2] = 2
+        A[2,2] = 1
+
+        #Right hand side matrix
+        v = array([[2,4],[3,4],[4,4]])
+
+        u = A*v
+        assert allclose(u, [[6,12], [14,16], [4,4]])  
+
 
 
