Index: /inundation/ga/storm_surge/pyvolution/least_squares.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/least_squares.py	(revision 440)
+++ /inundation/ga/storm_surge/pyvolution/least_squares.py	(revision 441)
@@ -20,4 +20,7 @@
 #FIXME: Current implementation uses full matrices and a general solver.
 #Later on we may consider using sparse techniques
+
+import exceptions
+class ShapeError(exceptions.Exception): pass
 
 from general_mesh import General_Mesh
@@ -178,33 +181,25 @@
         
     def build_coefficient_matrix_B(self, point_coordinates=None):
-        """Build final coefficient matrix
-        """
-        #self.alpha = 0
+        """Build final coefficient matrix"""
+        
+
         if self.alpha <> 0:
             self.build_smoothing_matrix_D()
         
-        
         if point_coordinates:
+
             self.build_interpolation_matrix_A(point_coordinates)
-            #self.At = self.At.tocsc()
-            #print "self.At ",self.At
-            #print "self.A",self.A
-            
-            #AtA = self.At.todense() * self.A.todense()
-            #print "AtA dense",AtA
-            #AtA = sparse.dok_matrix(AtA)
-            AtA = self.At * self.A 
-            #print "AtA",AtA
+
             if self.alpha <> 0:
-                self.B = AtA + self.alpha*self.D
+                self.B = self.AtA + self.alpha*self.D
             else:
-                self.B = AtA
-        #print "end of building b"
+                self.B = self.AtA
+
         
     def build_interpolation_matrix_A(self, point_coordinates):
         """Build n x m interpolation matrix, where
         n is the number of data points and
-        m is the number of basis functions phi_k (one per vertex)
-        """ 
+        m is the number of basis functions phi_k (one per vertex)"""
+        
         #Convert input to Numeric arrays
         point_coordinates = array(point_coordinates).astype(Float)
@@ -216,4 +211,5 @@
         #self.A = zeros((n,m), Float)
         self.A = sparse.dok_matrix()
+        self.AtA = sparse.dok_matrix()
 
         #Compute matrix elements
@@ -249,18 +245,27 @@
 
                     #Assign values to matrix A
-                    
+
+
                     j0 = self.mesh.triangles[k,0] #Global vertex id
-                    self.A[i, j0] = sigma0
+                    #self.A[i, j0] = sigma0
 
                     j1 = self.mesh.triangles[k,1] #Global vertex id
-                    self.A[i, j1] = sigma1
+                    #self.A[i, j1] = sigma1
 
                     j2 = self.mesh.triangles[k,2] #Global vertex id
-                    self.A[i, j2] = sigma2
-
-                    
-        #self.A = self.A.todense()        
-        #Precompute
-        self.At = self.A.transp()
+                    #self.A[i, j2] = sigma2
+
+                    sigmas = {j0:sigma0, j1:sigma1, j2:sigma2}
+                    js     = [j0,j1,j2]
+
+                    for j in js:
+                        self.A[i,j] = sigmas[j]
+                        for k in js:
+                            self.AtA[j,k] += sigmas[j]*sigmas[k]
+
+        
+        self.A   = (self.A).tocsc()
+        self.AtA = (self.AtA).tocsc()
+        self.At  =  self.A.transp()
 
     def get_A(self):
@@ -354,8 +359,9 @@
             self.D[v2,v0] += e20
             self.D[v0,v2] += e20             
-            
+
+        self.D = (self.D).tocsc()
             
     def fit(self, z):
-        """Fit a smooth surface to given data points z.
+        """Fit a smooth surface to given 1d array of data points z.
 
         The smooth surface is computed at each vertex in the underlying
@@ -366,10 +372,13 @@
           
         Inputs:
-          z: Vector or array of data at the point_coordinates.
-          If z is an array, smoothing will be done for each column
+          z: Single 1d vector or array of data at the point_coordinates.
         """
 
         #Convert input to Numeric arrays
         z = array(z).astype(Float)
+
+
+        if len(z.shape) > 1 :
+            raise VectorShapeError, 'Can only deal with 1d data vector'
         
         #Compute right hand side based on data
@@ -392,5 +401,5 @@
 
         return conjugate_gradient(self.B, Atz, Atz)
-        #FIXME: Should we store the result here for later use? (ON)
+        #FIXME: Should we store the result here for later use? (ON)        
     
     def fit_points(self, z):
@@ -432,7 +441,24 @@
           If f is an array, interpolation will be done for each column
         """
-        
-        return self.A * f
-
+
+        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
+            if m != M :
+                raise VectorShapeError, 'Mismatch between A and f dimensions'
+            
+            y = zeros( (N,n), Float)
+
+            for i in range(n):
+                y[:,i] = self.A * f[:,i]
+
+            return y
      
     #FIXME: We will need a method 'evaluate(self):' that will interpolate
Index: /inundation/ga/storm_surge/pyvolution/test_least_squares.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_least_squares.py	(revision 440)
+++ /inundation/ga/storm_surge/pyvolution/test_least_squares.py	(revision 441)
@@ -17,12 +17,13 @@
 class TestCase(unittest.TestCase):
 
-    def setUp(self):
-        pass
-        
-    def tearDown(self):
-        pass
+##     def setUp(self):
+##         pass
+        
+##     def tearDown(self):
+##         pass
 
 
     def test_datapoint_at_centroid(self):
+        
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -39,7 +40,6 @@
 
     def test_datapoints_at_vertices(self):
-        """Test that data points coinciding with vertices yield a diagonal matrix
-        """
-        print "aa"
+
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -58,9 +58,5 @@
 
     def test_datapoints_on_edge_midpoints(self):
-        """Try datapoints midway on edges -
-        each point should affect two matrix entries equally
-        """
-        print "bb"
-        
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -79,8 +75,4 @@
 
     def test_datapoints_on_edges(self):
-        """Try datapoints on edges -
-        each point should affect two matrix entries in proportion
-        """
-        print "cc"
         
         a = [0.0, 0.0]
@@ -99,6 +91,5 @@
 
     def test_arbitrary_datapoints(self):
-        """Try arbitrary datapoints
-        """
+
 
         from Numeric import sum
@@ -117,7 +108,6 @@
         
 
-    # this causes a memory error in scipy.sparce       
-    def bad_test_more_triangles(self):
-        print "dd"
+    def test_more_triangles(self):
+
         a = [-1.0, 0.0]
         b = [3.0, 4.0]
@@ -131,5 +121,5 @@
 
         #Data points
-        data_points = [ [-3., 1.9], [-2, 1], [0.0, 1], [0, 3], [2, 3], [-1.0/3,-4./3] ]
+        data_points = [ [-3., 1.9], [-2, 1], [0.0, 1], [0, 3], [2, 3], [-1.0/3,-4./3], [-1.0,-1.5 ], [1.0,-1.0]]
         interp = Interpolation(points, triangles, data_points) 
         
@@ -147,5 +137,5 @@
 
     def test_smooth_attributes_to_mesh(self):
-        print "ee"
+        
         a = [0.0, 0.0]
         b = [0.0, 5.0]
@@ -170,5 +160,5 @@
         
     def test_smooth_att_to_meshII(self):
-        print "ff"
+        
         a = [0.0, 0.0]
         b = [0.0, 5.0]
@@ -189,5 +179,5 @@
     
     def test_smooth_attributes_to_meshIII(self):
-        print "gg"
+
         a = [-1.0, 0.0]
         b = [3.0, 4.0]
@@ -221,8 +211,5 @@
 
     def test_smooth_attributes_to_meshIV(self):
-        print "hh"
-        """ Testing 2 attributes smoothed to the mesh
-        """
-        print "test_smooth_attributes_to_meshIV"
+
         a = [0.0, 0.0]
         b = [0.0, 5.0]
@@ -246,4 +233,5 @@
         
     def test_interpolate_attributes_to_points(self):
+        
         v0 = [0.0, 0.0]
         v1 = [0.0, 5.0]
@@ -268,4 +256,5 @@
        
     def test_interpolate_attributes_to_pointsII(self):
+
         a = [-1.0, 0.0]
         b = [3.0, 4.0]
@@ -301,4 +290,5 @@
     
     def test_interpolate_attributes_to_pointsIII(self):
+        
         v0 = [0.0, 0.0]
         v1 = [0.0, 5.0]
@@ -328,4 +318,5 @@
    
     def test_interpolate_attributes_to_pointsIV(self):
+
         a = [-1.0, 0.0]
         b = [3.0, 4.0]
@@ -366,8 +357,4 @@
         
     def test_smooth_attributes_to_mesh_function(self):
-        """ Testing 2 attributes smoothed to the mesh
-        """
-        """Test multiple attributes
-        """
         
         a = [0.0, 0.0]
@@ -433,4 +420,6 @@
 
     def test_smoothing_matrix_more_triangles(self):
+
+
         from Numeric import dot
         
@@ -476,4 +465,5 @@
 
     def test_fit_and_interpolation(self):
+
         from mesh import Mesh
         
@@ -663,5 +653,5 @@
 if __name__ == "__main__":
     suite = unittest.makeSuite(TestCase,'test')
-    runner = unittest.TextTestRunner()
+    runner = unittest.TextTestRunner(verbosity=2)
     runner.run(suite)
 
