source: anuga_core/source/anuga/utilities/test_cg_solve.py @ 3745

Last change on this file since 3745 was 3514, checked in by duncan, 18 years ago

Hi all,
I'm doing a change in the anuga structure, moving the code to

\anuga_core\source\anuga

After you have done an svn update, the PYTHONPATH has to be changed to;
PYTHONPATH = anuga_core/source/

This is part of changes required to make installation of anuga quicker and reducing the size of our sandpits.

If any imports are broken, try fixing them. With adding anuga. to them for example. If this seems to have really broken things, email/phone me.

Cheers
Duncan

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