Changeset 435


Ignore:
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:
2 edited

Legend:

Unmodified
Added
Removed
  • inundation/ga/storm_surge/pyvolution/cg_solve.py

    r432 r435  
     1import exceptions
     2class VectorShapeError(exceptions.Exception):
     3   """
     4   Setup Exception for routine conjugate_gradient
     5   """
     6   def __init__(self,args=None):
     7      self.args = args
     8
     9
     10
     11
     12
    113def conjugate_gradient(A,b,x0,imax=2000,tol=1.0e-8,iprint=0):
    214   #
     
    2032   #end
    2133
    22    from Numeric import dot
     34   from Numeric import dot, asarray, Float
    2335   #from operator import mod
    2436
    25    print "A shape",A.shape
    26    print "b shape",b.shape
    27    print "x0 shape",x0.shape
    28    i=0
     37   b  = asarray(b, typecode=Float)
     38   x0 = asarray(x0, typecode=Float)
     39
     40   print "A shape",A.shape[0]
     41   print "b shape",len(b.shape)
     42   print "x0 shape",x0.shape[0]
     43
     44   if len(b.shape) != 1 :
     45      raise VectorShapeError, 'input vector should consist of only one column'
     46
     47   if iprint == 0:
     48      iprint = imax
     49     
     50   i=1
    2951   x = x0
    3052   r = b - A*x
     
    4769       d = r + bt*d
    4870       i = i+1
    49        print "i",i
    50        #if i%iprint :
    51        #   #print 'i = %g rTr = %20.15e \n': i,rTr
    52        #   pass
     71       if i%iprint == 0 :
     72          print 'i = %g rTr = %20.15e'% (i,rTr)
     73
     74          
    5375   if i==imax:
    5476     print 'max number of iterations attained'
  • inundation/ga/storm_surge/pyvolution/test_cg_solve.py

    r432 r435  
    44#from math import sqrt
    55
    6 from Numeric import allclose, array, transpose
     6from Numeric import allclose, array, transpose, arange, ones, Float
    77
    88from cg_solve import *
     
    2626       
    2727        A = sparse.dok_matrix(A)
     28
     29        print "A"
     30        print A
    2831       
    2932        xe = [0.0, 1.0, 2.0, 3.0]
     
    3134        x =  [0.0, 0.0, 0.0, 0.0]
    3235
    33         x = conjugate_gradient(A,b,x)
     36        x = conjugate_gradient(A,b,x,iprint=1)
    3437
    3538        print "x",x
     
    3841
    3942
    40     def test_solve_more(self):
    41         A = [[2.0, -1.0, 0.0, 0.0 ],
    42              [-1.0, 2.0, -1.0, 0.0],
    43              [0.0, -1.0, 2.0, -1.0],
    44              [0.0,0.0, -1.0, 2.0]]
     43    def test_solve_large(self):
     44        A = sparse.dok_matrix()
     45
     46        n = 50
     47        for i in arange(0,n):
     48            A[i,i] = 1.0
     49            if i > 0 :
     50                A[i,i-1] = -0.5
     51            if i < n-1 :
     52                A[i,i+1] = -0.5
     53               
     54#        print "\n"
     55#        print "A\n"
     56#        print A
    4557       
    46         A = sparse.dok_matrix(A)
    4758       
    48         xe = [[0.0, 1.0], [1.0, 3.0], [2.0, 5.0], [3.0, 9.0]]
     59        xe = ones( (n,), Float)
     60
     61#        print xe
    4962        b  = A*xe
    50         x =  [[0.0, 0.0], [0.0, 0.0],[0.0, 0.0],[0.0, 0.0]]
    5163
    52         x = conjugate_gradient(A,b,x)
     64        x = conjugate_gradient(A,b,b,tol=1.0e-5,iprint=10)
    5365
    54         print "x",x
     66#        print "x",x
     67
     68        assert allclose(x,xe)
     69
     70    def test_solve_large_2d(self):
     71        A = sparse.dok_matrix()
     72
     73        n = 30
     74        m = 30
     75        for i in arange(0,n):
     76            for j in arange(0,m):
     77                I = j+m*i
     78                A[I,I] = 4.0
     79                if i > 0  :
     80                    A[I,I-m] = -1.0
     81                if i < n-1 :
     82                    A[I,I+m] = -1.0
     83                if j > 0  :
     84                    A[I,I-1] = -1.0
     85                if j < m-1 :
     86                    A[I,I+1] = -1.0
     87               
     88#        print "\n"
     89#        print "A\n"
     90#        print A.todense()
     91       
     92       
     93        xe = ones( (n*m,), Float)
     94
     95#        print xe
     96        b  = A*xe
     97
     98        x = conjugate_gradient(A,b,b,iprint=10)
     99
     100#        print "x",x
    55101
    56102        assert allclose(x,xe)
Note: See TracChangeset for help on using the changeset viewer.