Changeset 435

Timestamp:

Oct 21, 2004, 3:27:05 PM (20 years ago)

Author:

steve

Message:

Cleaned up conjugate_gradient

Location:

inundation/ga/storm_surge/pyvolution

Files:

• cg_solve.py (modified) (3 diffs)
• test_cg_solve.py (modified) (4 diffs)

inundation/ga/storm_surge/pyvolution/cg_solve.py

@@ +1 @@
+import exceptions
+class VectorShapeError(exceptions.Exception):
+   """
+   Setup Exception for routine conjugate_gradient
+   """
+   def __init__(self,args=None):
+      self.args = args
+
+
+
+
+
 def conjugate_gradient(A,b,x0,imax=2000,tol=1.0e-8,iprint=0):
    #
@@ -20 +32 @@
    #end
 
-   from Numeric import dot
+   from Numeric import dot, asarray, Float
    #from operator import mod
 
-   print "A shape",A.shape
-   print "b shape",b.shape
-   print "x0 shape",x0.shape
-   i=0
+   b  = asarray(b, typecode=Float)
+   x0 = asarray(x0, typecode=Float)
+
+   print "A shape",A.shape[0]
+   print "b shape",len(b.shape)
+   print "x0 shape",x0.shape[0]
+
+   if len(b.shape) != 1 :
+      raise VectorShapeError, 'input vector should consist of only one column'
+
+   if iprint == 0:
+      iprint = imax
+      
+   i=1
    x = x0
    r = b - A*x
@@ -47 +69 @@
        d = r + bt*d
        i = i+1
-       print "i",i
-       #if i%iprint :
-       #   #print 'i = %g rTr = %20.15e \n': i,rTr
-       #   pass
+       if i%iprint == 0 :
+          print 'i = %g rTr = %20.15e'% (i,rTr)
+
+          
    if i==imax:
      print 'max number of iterations attained'

inundation/ga/storm_surge/pyvolution/test_cg_solve.py

@@ -4 +4 @@
 #from math import sqrt
 
-from Numeric import allclose, array, transpose
+from Numeric import allclose, array, transpose, arange, ones, Float
 
 from cg_solve import *
@@ -26 +26 @@
         
         A = sparse.dok_matrix(A)
+
+        print "A"
+        print A
         
         xe = [0.0, 1.0, 2.0, 3.0]
@@ -31 +34 @@
         x =  [0.0, 0.0, 0.0, 0.0]
 
-        x = conjugate_gradient(A,b,x)
+        x = conjugate_gradient(A,b,x,iprint=1)
 
         print "x",x
@@ -38 +41 @@
 
 
-    def test_solve_more(self):
-        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]]
+    def test_solve_large(self):
+        A = sparse.dok_matrix()
+
+        n = 50
+        for i in arange(0,n):
+            A[i,i] = 1.0
+            if i > 0 :
+                A[i,i-1] = -0.5
+            if i < n-1 :
+                A[i,i+1] = -0.5
+                
+#        print "\n"
+#        print "A\n"
+#        print A
         
-        A = sparse.dok_matrix(A)
         
-        xe = [[0.0, 1.0], [1.0, 3.0], [2.0, 5.0], [3.0, 9.0]]
+        xe = ones( (n,), Float)
+
+#        print xe
         b  = A*xe
-        x =  [[0.0, 0.0], [0.0, 0.0],[0.0, 0.0],[0.0, 0.0]]
 
-        x = conjugate_gradient(A,b,x)
+        x = conjugate_gradient(A,b,b,tol=1.0e-5,iprint=10)
 
-        print "x",x
+#        print "x",x
+
+        assert allclose(x,xe)
+
+    def test_solve_large_2d(self):
+        A = sparse.dok_matrix()
+
+        n = 30
+        m = 30
+        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
+                
+#        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
 
         assert allclose(x,xe)