source: trunk/anuga_core/source/anuga/utilities/test_cg_solve.py @ 7848

Last change on this file since 7848 was 7848, checked in by steve, 13 years ago

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