← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

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

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

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()