← Back to branch summary

~ubuntu-branches/ubuntu/raring/python-scipy/raring-proposed

« back to all changes in this revision

Viewing changes to Lib/sandbox/pysparse/tests/test_superlu.py

  • Committer: Bazaar Package Importer
  • Author(s): Matthias Klose
  • Date: 2007-01-07 14:12:12 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20070107141212-mm0ebkh5b37hcpzn
http://bugs.debian.org/402783
* Remove build dependency on python-numpy-dev.
* python-scipy: Depend on python-numpy instead of python-numpy-dev.
* Package builds on other archs than i386. Closes: #402783.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/sandbox/pysparse/tests/test_superlu.py
 
1
import math, unittest
 
2
import scipy
 
3
import spmatrix, superlu
 
4
import poisson
 
5
 
 
6
def residual(A, x, b):
 
7
    n = A.shape[0]
 
8
    r = numpy.zeros(n, 'd')
 
9
    A.matvec(x, r)
 
10
    r -= b
 
11
    return math.sqrt(numpy.dot(r, r))
 
12
 
 
13
def error(x, y):
 
14
    n = len(x)
 
15
    t = x.copy(); t -= y
 
16
    return math.sqrt(numpy.dot(t, t))
 
17
 
 
18
def speye(n):
 
19
    A = spmatrix.ll_mat_sym(n, n)
 
20
    for i in xrange(n):
 
21
        A[i,i] = 1.0
 
22
    return A
 
23
 
 
24
def macheps():
 
25
    "compute machine epsilon"
 
26
    eps = 1.0
 
27
    while (1.0 + eps > 1.0):
 
28
        eps /= 2.0
 
29
    return 2.0 * eps
 
30
 
 
31
class EyeTestCase(unittest.TestCase):
 
32
    def setUp(self):
 
33
        self.n = 10000
 
34
        self.A = speye(self.n).to_csr()
 
35
        self.b = numpy.ones(self.n, 'd')
 
36
        self.x = numpy.zeros(self.n, 'd')
 
37
 
 
38
    def testTrivial(self):
 
39
        luA = superlu.factorize(self.A)
 
40
        luA.solve(self.b, self.x)
 
41
        self.failUnless(residual(self.A, self.x, self.b) == 0.0)
 
42
 
 
43
        luA = superlu.factorize(self.A, diag_pivot_thresh=0.0)
 
44
        luA.solve(self.b, self.x)
 
45
        self.failUnless(residual(self.A, self.x, self.b) == 0.0)
 
46
 
 
47
        luA = superlu.factorize(self.A, relax=20)
 
48
        luA.solve(self.b, self.x)
 
49
        self.failUnless(residual(self.A, self.x, self.b) == 0.0)
 
50
 
 
51
        luA = superlu.factorize(self.A, panel_size=1)
 
52
        luA.solve(self.b, self.x)
 
53
        self.failUnless(residual(self.A, self.x, self.b) == 0.0)
 
54
 
 
55
        for permc_spec in [0, 1, 2, 3]:
 
56
            luA = superlu.factorize(self.A, permc_spec=permc_spec)
 
57
            luA.solve(self.b, self.x)
 
58
            self.failUnless(residual(self.A, self.x, self.b) == 0.0)
 
59
 
 
60
class Poisson1dTestCase(unittest.TestCase):
 
61
    def setUp(self):
 
62
        self.n = 50000
 
63
        self.B = poisson.poisson1d(self.n).to_csr()
 
64
 
 
65
        self.b = numpy.zeros(self.n, 'd')
 
66
        self.x = numpy.zeros(self.n, 'd')
 
67
        self.x_exact = numpy.ones(self.n, 'd')
 
68
        self.x_exact /= math.sqrt(self.n)
 
69
        self.B.matvec(self.x_exact, self.b)
 
70
 
 
71
        lmbd_min = 4.0 * math.sin(math.pi/2.0/self.n) ** 2
 
72
        lmbd_max = 4.0 * math.sin((self.n - 1)*math.pi/2.0/self.n) ** 2
 
73
        cond = lmbd_max/lmbd_min
 
74
        self.tol = cond * macheps()
 
75
 
 
76
    def testPoisson1dDefault(self):
 
77
        luA = superlu.factorize(self.B)
 
78
        luA.solve(self.b, self.x)
 
79
        print error(self.x, self.x_exact), self.tol, luA.nnz
 
80
        self.failUnless(error(self.x, self.x_exact) < self.tol)
 
81
 
 
82
    def testPoisson1dNoPivot(self):
 
83
        luA = superlu.factorize(self.B, diag_pivot_thresh=0.0)
 
84
        luA.solve(self.b, self.x)
 
85
        print error(self.x, self.x_exact), self.tol, luA.nnz
 
86
        self.failUnless(error(self.x, self.x_exact) < self.tol)
 
87
 
 
88
    def testPoisson1dRelax(self):
 
89
        luA = superlu.factorize(self.B, relax=20)
 
90
        luA.solve(self.b, self.x)
 
91
        print error(self.x, self.x_exact), self.tol, luA.nnz
 
92
        self.failUnless(error(self.x, self.x_exact) < self.tol)
 
93
 
 
94
    def testPoisson1dSmallPanel(self):
 
95
        luA = superlu.factorize(self.B, panel_size=1)
 
96
        luA.solve(self.b, self.x)
 
97
        print error(self.x, self.x_exact), self.tol, luA.nnz
 
98
        self.failUnless(error(self.x, self.x_exact) < self.tol)
 
99
 
 
100
    def testPoisson1dOrigOrdering(self):
 
101
        luA = superlu.factorize(self.B, permc_spec=0)
 
102
        luA.solve(self.b, self.x)
 
103
        print error(self.x, self.x_exact), self.tol, luA.nnz
 
104
        self.failUnless(error(self.x, self.x_exact) < self.tol)
 
105
 
 
106
    def testPoisson1dMMD_AtimesA(self):
 
107
        luA = superlu.factorize(self.B, permc_spec=1)
 
108
        luA.solve(self.b, self.x)
 
109
        print error(self.x, self.x_exact), self.tol, luA.nnz
 
110
        self.failUnless(error(self.x, self.x_exact) < self.tol)
 
111
 
 
112
    def testPoisson1dMMD_AplusA(self):
 
113
        luA = superlu.factorize(self.B, permc_spec=2)
 
114
        luA.solve(self.b, self.x)
 
115
        print error(self.x, self.x_exact), self.tol, luA.nnz
 
116
        self.failUnless(error(self.x, self.x_exact) < self.tol)
 
117
 
 
118
    def testPoisson1dCOLAMD(self):
 
119
        luA = superlu.factorize(self.B, permc_spec=3)
 
120
        luA.solve(self.b, self.x)
 
121
        print error(self.x, self.x_exact), self.tol, luA.nnz
 
122
        self.failUnless(error(self.x, self.x_exact) < self.tol)
 
123
 
 
124
class Poisson2dTestCase(unittest.TestCase):
 
125
    def setUp(self):
 
126
        self.n = 200
 
127
        self.B = poisson.poisson2d(self.n).to_csr()
 
128
 
 
129
        self.b = numpy.zeros(self.n*self.n, 'd')
 
130
        self.x = numpy.zeros(self.n*self.n, 'd')
 
131
        self.x_exact = numpy.ones(self.n*self.n, 'd')
 
132
        self.x_exact /= math.sqrt(self.n*self.n)
 
133
        self.B.matvec(self.x_exact, self.b)
 
134
 
 
135
        h = 1.0 / self.n
 
136
        lmbd_min = 4.0/h/h * (math.sin(math.pi*h/2.0) ** 2 +
 
137
                              math.sin(math.pi*h/2.0) ** 2)
 
138
        lmbd_max = 4.0/h/h * (math.sin((self.n - 1)*math.pi*h/2.0) ** 2 +
 
139
                              math.sin((self.n - 1)*math.pi*h/2.0) ** 2)
 
140
        cond = lmbd_max/lmbd_min
 
141
        self.tol = cond * macheps()
 
142
 
 
143
    def testPoisson2dDefault(self):
 
144
        luA = superlu.factorize(self.B)
 
145
        luA.solve(self.b, self.x)
 
146
        print error(self.x, self.x_exact), self.tol, luA.nnz
 
147
        self.failUnless(error(self.x, self.x_exact) < self.tol)
 
148
 
 
149
    def testPoisson2dNoPivot(self):
 
150
        luA = superlu.factorize(self.B, diag_pivot_thresh=0.0)
 
151
        luA.solve(self.b, self.x)
 
152
        print error(self.x, self.x_exact), self.tol, luA.nnz
 
153
        self.failUnless(error(self.x, self.x_exact) < self.tol)
 
154
 
 
155
    def testPoisson2dRelax(self):
 
156
        luA = superlu.factorize(self.B, relax=20)
 
157
        luA.solve(self.b, self.x)
 
158
        print error(self.x, self.x_exact), self.tol, luA.nnz
 
159
        self.failUnless(error(self.x, self.x_exact) < self.tol)
 
160
 
 
161
    def testPoisson2dSmallPanel(self):
 
162
        luA = superlu.factorize(self.B, panel_size=1)
 
163
        luA.solve(self.b, self.x)
 
164
        print error(self.x, self.x_exact), self.tol, luA.nnz
 
165
        self.failUnless(error(self.x, self.x_exact) < self.tol)
 
166
 
 
167
    def testPoisson2dOrigOrdering(self):
 
168
        luA = superlu.factorize(self.B, permc_spec=0)
 
169
        luA.solve(self.b, self.x)
 
170
        print error(self.x, self.x_exact), self.tol, luA.nnz
 
171
        self.failUnless(error(self.x, self.x_exact) < self.tol)
 
172
 
 
173
    def testPoisson2dMMD_AtimesA(self):
 
174
        luA = superlu.factorize(self.B, permc_spec=1)
 
175
        luA.solve(self.b, self.x)
 
176
        print error(self.x, self.x_exact), self.tol, luA.nnz
 
177
        self.failUnless(error(self.x, self.x_exact) < self.tol)
 
178
 
 
179
    def testPoisson2dMMD_AplusA(self):
 
180
        luA = superlu.factorize(self.B, permc_spec=2)
 
181
        luA.solve(self.b, self.x)
 
182
        print error(self.x, self.x_exact), self.tol, luA.nnz
 
183
        self.failUnless(error(self.x, self.x_exact) < self.tol)
 
184
 
 
185
    def testPoisson2dCOLAMD(self):
 
186
        luA = superlu.factorize(self.B, permc_spec=3)
 
187
        luA.solve(self.b, self.x)
 
188
        print error(self.x, self.x_exact), self.tol, luA.nnz
 
189
        self.failUnless(error(self.x, self.x_exact) < self.tol)
 
190
 
 
191
if __name__ == '__main__':
 
192
    unittest.main()