Changeset 438 for inundation/ga/storm_surge/pyvolution

Timestamp:

Oct 21, 2004, 11:54:49 PM (20 years ago)

Author:

steve

Message:

CG tests

Location:

inundation/ga/storm_surge/pyvolution

Files:

• cg_solve.py (modified) (2 diffs)
• test_cg_solve.py (modified) (5 diffs)

inundation/ga/storm_surge/pyvolution/cg_solve.py

@@ -1 +1 @@
 import exceptions
-class VectorShapeError(exceptions.Exception):
-   """
-   Setup Exception for routine conjugate_gradient
-   """
-   def __init__(self,args=None):
-      self.args = args
-
-
-
+class VectorShapeError(exceptions.Exception): pass
 
 
 def conjugate_gradient(A,b,x0,imax=2000,tol=1.0e-8,iprint=0):
-   #
-   # Try to solve linear equation Ax = b using
-   # conjugate gradient method
-   #
-   # Input
-   # A: matrix or function which applies a matrix, assumed symmetric
-   # b: right hand side
-   # x0: inital guess
-   # imax: max number of iterations
-   # tol: tolerance used for residual
-   #
-   # Output
-   # x: approximate solution
-   # i: number of iterations
-   #
+   """
+   Try to solve linear equation Ax = b using
+   conjugate gradient method
+   
+   Input
+   A: matrix or function which applies a matrix, assumed symmetric
+   b: right hand side
+   x0: inital guess
+   imax: max number of iterations
+   tol: tolerance used for residual
+   
+   Output
+   x: approximate solution
+   """
 
    #if nargin<3
@@ -37 +28 @@
    b  = asarray(b, typecode=Float)
    x0 = asarray(x0, typecode=Float)
+   #A = asarray(A, typecode=Float)
 
-   print "A shape",A.shape[0]
-   print "b shape",len(b.shape)
-   print "x0 shape",x0.shape[0]
+   #print "A shape",A.shape
+   #print "b shape",len(b.shape)
+   #print "x0 shape",x0.shape[0]
 
    if len(b.shape) != 1 :

inundation/ga/storm_surge/pyvolution/test_cg_solve.py

@@ -2 +2 @@
 
 import unittest
-#from math import sqrt
 
-from Numeric import allclose, array, transpose, arange, ones, Float
+from Matrix import *
+from Numeric import dot, allclose, array, transpose, arange, ones, Float
+from cg_solve import *
+from scipy import sparse, mat
 
-from cg_solve import *
 
-from scipy import sparse
-        
 class TestCase(unittest.TestCase):
 
-    def setUp(self):
-        pass
+##     def test_solve(self):
+##         """Small Matrix"""
         
-    def tearDown(self):
-        pass
+##         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)
 
-    def test_solve(self):
+##         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):
+        """Small Sparse Matrix"""
+        
         A = [[2.0, -1.0, 0.0, 0.0 ],
              [-1.0, 2.0, -1.0, 0.0],
@@ -27 +44 @@
         A = sparse.dok_matrix(A)
 
-        print "A"
-        print A
-        
         xe = [0.0, 1.0, 2.0, 3.0]
         b  = A*xe
         x =  [0.0, 0.0, 0.0, 0.0]
 
-        x = conjugate_gradient(A,b,x,iprint=1)
-
-        print "x",x
+        x = conjugate_gradient(A,b,x,iprint=0)
 
         assert allclose(x,xe)
@@ -42 +54 @@
 
     def test_solve_large(self):
+        """Standard 1d laplacian """
+
         A = sparse.dok_matrix()
 
@@ -52 +66 @@
                 A[i,i+1] = -0.5
                 
-#        print "\n"
-#        print "A\n"
-#        print A
-        
-        
         xe = ones( (n,), Float)
 
-#        print xe
         b  = A*xe
-
-        x = conjugate_gradient(A,b,b,tol=1.0e-5,iprint=10)
-
-#        print "x",x
+        x = conjugate_gradient(A,b,b,tol=1.0e-5,iprint=0)
 
         assert allclose(x,xe)
 
     def test_solve_large_2d(self):
+        """Standard 2d laplacian"""
+        
         A = sparse.dok_matrix()
 
-        n = 30
-        m = 30
+        n = 20
+        m = 10
         for i in arange(0,n):
             for j in arange(0,m):
@@ -86 +93 @@
                     A[I,I+1] = -1.0
                 
-#        print "\n"
-#        print "A\n"
-#        print A.todense()
-        
-        
         xe = ones( (n*m,), Float)
 
-#        print xe
         b  = A*xe
-
-        x = conjugate_gradient(A,b,b,iprint=10)
-
-#        print "x",x
+        x = conjugate_gradient(A,b,b,iprint=0)
 
         assert allclose(x,xe)
+
+
+    def test_vector_shape_error(self):
+        """Raise VectorShapeError"""
+        
+        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 = sparse.dok_matrix(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)
+        except:
+            pass
+        else:
+            msg = 'Should have raised exception'
+            raise msg
 
         
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    suite = unittest.makeSuite(TestCase,'test')
-    runner = unittest.TextTestRunner()
-    runner.run(suite)
+     suite = unittest.makeSuite(TestCase,'test')
+     runner = unittest.TextTestRunner(verbosity=2)
+     runner.run(suite)