source: inundation-numpy-branch/utilities/test_cg_solve.py @ 2545

Last change on this file since 2545 was 2545, checked in by ole, 18 years ago

Got utilities to compile and pass all tests using numpy

File size: 4.3 KB
Line 
1#!/usr/bin/env python
2
3import unittest
4
5
6from numpy import dot, allclose, array, transpose, arange, ones, Float
7from cg_solve import *
8from sparse import Sparse, Sparse_CSR
9
10
11class Test_CG_Solve(unittest.TestCase):
12
13    def test_sparse_solve(self):
14        """Small Sparse 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 = Sparse(A)
22
23        xe = [0.0, 1.0, 2.0, 3.0]
24        b  = A*xe
25        x =  [0.0, 0.0, 0.0, 0.0]
26
27        x = conjugate_gradient(A,b,x,iprint=0)
28
29        assert allclose(x,xe)
30
31    def test_max_iter(self):
32        """Small Sparse Matrix"""
33
34        A = [[2.0, -1.0, 0.0, 0.0 ],
35             [-1.0, 2.0, -1.0, 0.0],
36             [0.0, -1.0, 2.0, -1.0],
37             [0.0,0.0, -1.0, 2.0]]
38
39        A = Sparse(A)
40
41        xe = [0.0, 1.0, 2.0, 3.0]
42        b  = A*xe
43        x =  [0.0, 0.0, 0.0, 0.0]
44
45        x = conjugate_gradient(A,b,x,iprint=0,imax=2)
46
47
48    def test_solve_large(self):
49        """Standard 1d laplacian """
50
51        n = 50
52        A = Sparse(n,n)
53
54        for i in arange(0,n):
55            A[i,i] = 1.0
56            if i > 0 :
57                A[i,i-1] = -0.5
58            if i < n-1 :
59                A[i,i+1] = -0.5
60
61        xe = ones( (n,), Float)
62
63        b  = A*xe
64        x = conjugate_gradient(A,b,b,tol=1.0e-5,iprint=0)
65
66        assert allclose(x,xe)
67
68    def test_solve_large_2d(self):
69        """Standard 2d laplacian"""
70
71        n = 20
72        m = 10
73
74        A = Sparse(m*n, m*n)
75
76        for i in arange(0,n):
77            for j in arange(0,m):
78                I = j+m*i
79                A[I,I] = 4.0
80                if i > 0  :
81                    A[I,I-m] = -1.0
82                if i < n-1 :
83                    A[I,I+m] = -1.0
84                if j > 0  :
85                    A[I,I-1] = -1.0
86                if j < m-1 :
87                    A[I,I+1] = -1.0
88
89        xe = ones( (n*m,), Float)
90
91        b  = A*xe
92        x = conjugate_gradient(A,b,b,iprint=0)
93
94        assert allclose(x,xe)
95
96    def test_solve_large_2d_csr_matrix(self):
97        """Standard 2d laplacian with csr format
98        """
99
100        n = 100
101        m = 100
102
103        A = Sparse(m*n, m*n)
104
105        for i in arange(0,n):
106            for j in arange(0,m):
107                I = j+m*i
108                A[I,I] = 4.0
109                if i > 0  :
110                    A[I,I-m] = -1.0
111                if i < n-1 :
112                    A[I,I+m] = -1.0
113                if j > 0  :
114                    A[I,I-1] = -1.0
115                if j < m-1 :
116                    A[I,I+1] = -1.0
117
118        xe = ones( (n*m,), Float)
119
120        # Convert to csr format
121        #print 'start covert'
122        A = Sparse_CSR(A)
123        #print 'finish covert'
124        b = A*xe
125        x = conjugate_gradient(A,b,b,iprint=20)
126
127        assert allclose(x,xe)
128
129
130    def test_solve_large_2d_with_default_guess(self):
131        """Standard 2d laplacian using default first guess"""
132
133        n = 20
134        m = 10
135
136        A = Sparse(m*n, m*n)
137
138        for i in arange(0,n):
139            for j in arange(0,m):
140                I = j+m*i
141                A[I,I] = 4.0
142                if i > 0  :
143                    A[I,I-m] = -1.0
144                if i < n-1 :
145                    A[I,I+m] = -1.0
146                if j > 0  :
147                    A[I,I-1] = -1.0
148                if j < m-1 :
149                    A[I,I+1] = -1.0
150
151        xe = ones( (n*m,), Float)
152
153        b  = A*xe
154        x = conjugate_gradient(A,b)
155
156        assert allclose(x,xe)
157
158
159    def test_vector_shape_error(self):
160        """Raise VectorShapeError"""
161
162        A = [[2.0, -1.0, 0.0, 0.0 ],
163             [-1.0, 2.0, -1.0, 0.0],
164             [0.0, -1.0, 2.0, -1.0],
165             [0.0,0.0, -1.0, 2.0]]
166
167        A = Sparse(A)
168
169        xe = [[0.0,2.0], [1.0,3.0], [2.0,4.0], [3.0,2.0]]
170
171        try:
172            x = conjugate_gradient(A,xe,xe,iprint=0)
173        except:
174            pass
175        else:
176            msg = 'Should have raised exception'
177            raise msg
178
179
180#-------------------------------------------------------------
181if __name__ == "__main__":
182     suite = unittest.makeSuite(Test_CG_Solve,'test')
183     runner = unittest.TextTestRunner() #(verbosity=2)
184     runner.run(suite)
185
186
187
Note: See TracBrowser for help on using the repository browser.