Changeset 435 for inundation/ga/storm_surge/pyvolution/cg_solve.py

Timestamp:

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

Author:

steve

Message:

Cleaned up conjugate_gradient

File:

• inundation/ga/storm_surge/pyvolution/cg_solve.py (modified) (3 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'