Changeset 435
- Timestamp:
- Oct 21, 2004, 3:27:05 PM (20 years ago)
- Location:
- inundation/ga/storm_surge/pyvolution
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
inundation/ga/storm_surge/pyvolution/cg_solve.py
r432 r435 1 import exceptions 2 class 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 1 13 def conjugate_gradient(A,b,x0,imax=2000,tol=1.0e-8,iprint=0): 2 14 # … … 20 32 #end 21 33 22 from Numeric import dot 34 from Numeric import dot, asarray, Float 23 35 #from operator import mod 24 36 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 29 51 x = x0 30 52 r = b - A*x … … 47 69 d = r + bt*d 48 70 i = i+1 49 print "i",i50 #if i%iprint :51 # #print 'i = %g rTr = %20.15e \n': i,rTr 52 # pass71 if i%iprint == 0 : 72 print 'i = %g rTr = %20.15e'% (i,rTr) 73 74 53 75 if i==imax: 54 76 print 'max number of iterations attained' -
inundation/ga/storm_surge/pyvolution/test_cg_solve.py
r432 r435 4 4 #from math import sqrt 5 5 6 from Numeric import allclose, array, transpose 6 from Numeric import allclose, array, transpose, arange, ones, Float 7 7 8 8 from cg_solve import * … … 26 26 27 27 A = sparse.dok_matrix(A) 28 29 print "A" 30 print A 28 31 29 32 xe = [0.0, 1.0, 2.0, 3.0] … … 31 34 x = [0.0, 0.0, 0.0, 0.0] 32 35 33 x = conjugate_gradient(A,b,x )36 x = conjugate_gradient(A,b,x,iprint=1) 34 37 35 38 print "x",x … … 38 41 39 42 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 45 57 46 A = sparse.dok_matrix(A)47 58 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 49 62 b = A*xe 50 x = [[0.0, 0.0], [0.0, 0.0],[0.0, 0.0],[0.0, 0.0]]51 63 52 x = conjugate_gradient(A,b, x)64 x = conjugate_gradient(A,b,b,tol=1.0e-5,iprint=10) 53 65 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 55 101 56 102 assert allclose(x,xe)
Note: See TracChangeset
for help on using the changeset viewer.