Changeset 475

Timestamp:

Nov 1, 2004, 3:30:45 PM (20 years ago)

Author:

ole

Message:

Modified test_cg to use sparse.py

Location:

inundation/ga/storm_surge/pyvolution

Files:

• cg_solve.py (modified) (3 diffs)
• test_cg_solve.py (modified) (6 diffs)
• test_least_squares.py (modified) (1 diff)

inundation/ga/storm_surge/pyvolution/cg_solve.py

@@ -3 +3 @@
 
 
-def conjugate_gradient(A,b,x0,imax=10000,tol=1.0e-8,iprint=0):
+def conjugate_gradient(A,b,x0=None,imax=10000,tol=1.0e-8,iprint=0):
    """
    Try to solve linear equation Ax = b using
@@ -10 +10 @@
    Input
    A: matrix or function which applies a matrix, assumed symmetric
+      A can be either dense or sparse
    b: right hand side
-   x0: inital guess
+   x0: inital guess (default the 0 vector)
    imax: max number of iterations
    tol: tolerance used for residual
@@ -19 +20 @@
    """
 
-   #if nargin<3
-   #   x0 = zeros(b);
-   #end
-
-   from Numeric import dot, asarray, Float
-   #from operator import mod
-
-   b  = asarray(b, typecode=Float)
-   x0 = asarray(x0, typecode=Float)
-   #A = asarray(A, typecode=Float)
-
-   #print "A shape",A.shape
-   #print "b shape",len(b.shape)
-   #print "x0 shape",x0.shape[0]
+   from Numeric import dot, array, Float, zeros
+   
+   b  = array(b, typecode=Float)
+   
+   if x0 is None:
+      x0 = zeros(b.shape, typecode=Float)
+   else:   
+      x0 = array(x0, typecode=Float)
+      
 
    if len(b.shape) != 1 :
       raise VectorShapeError, 'input vector should consist of only one column'
 
+   #FIXME: Should test using None
    if iprint == 0:
       iprint = imax

inundation/ga/storm_surge/pyvolution/test_cg_solve.py

@@ -3 +3 @@
 import unittest
 
-from Matrix import *
+
 from Numeric import dot, allclose, array, transpose, arange, ones, Float
 from cg_solve import *
-from scipy import sparse, mat
+from sparse import Sparse
 
 
 class TestCase(unittest.TestCase):
-
-##     def test_solve(self):
-##         """Small Matrix"""
-        
-##         A = [[2.0, -1.0, 0.0, 0.0 ],
-##              [-1.0, 2.0, -1.0, 0.0],
-##              [0.0, -1.0, 2.0, -1.0],
-##              [0.0, 0.0, -1.0, 2.0]]
-
-##         A = mat(A)
-
-##         xe = asarray([0.0, 1.0, 2.0,  3.0])
-##         b  = A*xe
-##         x =  asarray([0.0, 0.0, 0.0, 0.0])
-
-##         print A
-##         print xe
-##         print b
-##         print b.shape
-
-##         x = conjugate_gradient(A,b,x,iprint=1)
-
-##         assert allclose(x,xe)
 
     def test_sparse_solve(self):
@@ -42 +19 @@
              [0.0,0.0, -1.0, 2.0]]
         
-        A = sparse.dok_matrix(A)
+        A = Sparse(A)
 
         xe = [0.0, 1.0, 2.0, 3.0]
@@ -56 +33 @@
         """Standard 1d laplacian """
 
-        A = sparse.dok_matrix()
-
         n = 50
+        A = Sparse(n,n)
+        
         for i in arange(0,n):
             A[i,i] = 1.0
@@ -76 +53 @@
         """Standard 2d laplacian"""
         
-        A = sparse.dok_matrix()
-
         n = 20
         m = 10
+
+        A = Sparse(m*n, m*n)
+
         for i in arange(0,n):
             for j in arange(0,m):
@@ -101 +79 @@
 
 
+    def test_solve_large_2d_with_default_guess(self):
+        """Standard 2d laplacian using default first guess"""
+        
+        n = 20
+        m = 10
+
+        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)
+
+        b  = A*xe
+        x = conjugate_gradient(A,b)
+
+        assert allclose(x,xe)
+
+
     def test_vector_shape_error(self):
         """Raise VectorShapeError"""
@@ -109 +116 @@
              [0.0,0.0, -1.0, 2.0]]
         
-        A = sparse.dok_matrix(A)
+        A = Sparse(A)
 
         xe = [[0.0,2.0], [1.0,3.0], [2.0,4.0], [3.0,2.0]]
-        b  = A*xe
-        x =  [0.0, 0.0, 0.0, 0.0]
 
         try:
-            x = conjugate_gradient(A,b,b,iprint=0)
+            x = conjugate_gradient(A,xe,xe,iprint=0)
         except:
             pass

inundation/ga/storm_surge/pyvolution/test_least_squares.py

@@ -671 +671 @@
 if __name__ == "__main__":
     suite = unittest.makeSuite(TestCase,'test')
-    runner = unittest.TextTestRunner()   #(verbosity=2)
+    runner = unittest.TextTestRunner(verbosity=1)
     runner.run(suite)