source: inundation/ga/storm_surge/pyvolution/test_cg_solve.py @ 474

Last change on this file since 474 was 454, checked in by duncan, 21 years ago

bug fixes

File size: 3.0 KB
Line 
1#!/usr/bin/env python
2
3import unittest
4
5from Matrix import *
6from Numeric import dot, allclose, array, transpose, arange, ones, Float
7from cg_solve import *
8from scipy import sparse, mat
9
10
11class TestCase(unittest.TestCase):
12
13##     def test_solve(self):
14##         """Small Matrix"""
15       
16##         A = [[2.0, -1.0, 0.0, 0.0 ],
17##              [-1.0, 2.0, -1.0, 0.0],
18##              [0.0, -1.0, 2.0, -1.0],
19##              [0.0, 0.0, -1.0, 2.0]]
20
21##         A = mat(A)
22
23##         xe = asarray([0.0, 1.0, 2.0,  3.0])
24##         b  = A*xe
25##         x =  asarray([0.0, 0.0, 0.0, 0.0])
26
27##         print A
28##         print xe
29##         print b
30##         print b.shape
31
32##         x = conjugate_gradient(A,b,x,iprint=1)
33
34##         assert allclose(x,xe)
35
36    def test_sparse_solve(self):
37        """Small Sparse Matrix"""
38       
39        A = [[2.0, -1.0, 0.0, 0.0 ],
40             [-1.0, 2.0, -1.0, 0.0],
41             [0.0, -1.0, 2.0, -1.0],
42             [0.0,0.0, -1.0, 2.0]]
43       
44        A = sparse.dok_matrix(A)
45
46        xe = [0.0, 1.0, 2.0, 3.0]
47        b  = A*xe
48        x =  [0.0, 0.0, 0.0, 0.0]
49
50        x = conjugate_gradient(A,b,x,iprint=0)
51
52        assert allclose(x,xe)
53
54
55    def test_solve_large(self):
56        """Standard 1d laplacian """
57
58        A = sparse.dok_matrix()
59
60        n = 50
61        for i in arange(0,n):
62            A[i,i] = 1.0
63            if i > 0 :
64                A[i,i-1] = -0.5
65            if i < n-1 :
66                A[i,i+1] = -0.5
67               
68        xe = ones( (n,), Float)
69
70        b  = A*xe
71        x = conjugate_gradient(A,b,b,tol=1.0e-5,iprint=0)
72
73        assert allclose(x,xe)
74
75    def test_solve_large_2d(self):
76        """Standard 2d laplacian"""
77       
78        A = sparse.dok_matrix()
79
80        n = 20
81        m = 10
82        for i in arange(0,n):
83            for j in arange(0,m):
84                I = j+m*i
85                A[I,I] = 4.0
86                if i > 0  :
87                    A[I,I-m] = -1.0
88                if i < n-1 :
89                    A[I,I+m] = -1.0
90                if j > 0  :
91                    A[I,I-1] = -1.0
92                if j < m-1 :
93                    A[I,I+1] = -1.0
94               
95        xe = ones( (n*m,), Float)
96
97        b  = A*xe
98        x = conjugate_gradient(A,b,b,iprint=0)
99
100        assert allclose(x,xe)
101
102
103    def test_vector_shape_error(self):
104        """Raise VectorShapeError"""
105       
106        A = [[2.0, -1.0, 0.0, 0.0 ],
107             [-1.0, 2.0, -1.0, 0.0],
108             [0.0, -1.0, 2.0, -1.0],
109             [0.0,0.0, -1.0, 2.0]]
110       
111        A = sparse.dok_matrix(A)
112
113        xe = [[0.0,2.0], [1.0,3.0], [2.0,4.0], [3.0,2.0]]
114        b  = A*xe
115        x =  [0.0, 0.0, 0.0, 0.0]
116
117        try:
118            x = conjugate_gradient(A,b,b,iprint=0)
119        except:
120            pass
121        else:
122            msg = 'Should have raised exception'
123            raise msg
124
125       
126#-------------------------------------------------------------
127if __name__ == "__main__":
128     suite = unittest.makeSuite(TestCase,'test')
129     runner = unittest.TextTestRunner() #(verbosity=2)
130     runner.run(suite)
131
132   
133   
134
135
136
Note: See TracBrowser for help on using the repository browser.