← 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/examples/jdsym_test.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/examples/jdsym_test.py
 
1
import spmatrix, jdsym, itsolvers
 
2
from numpy import zeros, dot, allclose, multiply
 
3
from math import sqrt
 
4
import RandomArray
 
5
 
 
6
class diagPrecShifted:
 
7
    def __init__(self, A, M, sigma):
 
8
        self.shape = A.shape
 
9
        n = self.shape[0]
 
10
        self.dinv = zeros(n, 'd')
 
11
        for i in xrange(n):
 
12
            self.dinv[i] = 1.0 / (A[i,i] - sigma*M[i,i])
 
13
    def precon(self, x, y):
 
14
        multiply(x, self.dinv, y)
 
15
 
 
16
def computeResiduals(A, M, lmbd, Q):
 
17
    kconv = lmbd.shape[0]
 
18
    residuals = zeros((kconv, ), 'd')
 
19
    r = zeros((n, ), 'd')
 
20
    u = zeros((n, ), 'd')
 
21
    t = zeros((n, ), 'd')
 
22
    for k in xrange(kconv):
 
23
        u = Q[:,k].copy()
 
24
        A.matvec(u, r)
 
25
        if M <> None:
 
26
            M.matvec(u, t)
 
27
        else:
 
28
            t = u
 
29
        r = r - lmbd[k]*t
 
30
        residuals[k] = sqrt(dot(r,r))
 
31
    return residuals
 
32
 
 
33
n = 1000; ncv = 5; tol = 1e-6
 
34
 
 
35
A = spmatrix.ll_mat_sym(n)
 
36
for i in xrange(n):
 
37
    A[i,i] = i+1.0
 
38
As = A.to_sss()
 
39
 
 
40
M = spmatrix.ll_mat_sym(n)
 
41
for i in xrange(n):
 
42
    M[i,i] = float(n/2) + i
 
43
Ms = M.to_sss()
 
44
normM = M[n-1,n-1]
 
45
 
 
46
K = diagPrecShifted(A, M, 0.006)
 
47
 
 
48
#-------------------------------------------------------------------------------
 
49
# Test 1: M = K = None
 
50
 
 
51
print 'Test 1'
 
52
 
 
53
lmbd_exact = zeros(ncv, 'd')
 
54
for k in xrange(ncv):
 
55
    lmbd_exact[k] =  A[k,k]
 
56
 
 
57
kconv, lmbd, Q, it, it_inner = jdsym.jdsym(As, None, None, ncv, 0.0, tol, 150, itsolvers.qmrs,
 
58
                                           jmin=5, jmax=10, eps_tr=1e-4, clvl=1)
 
59
 
 
60
assert ncv == kconv
 
61
assert allclose(computeResiduals(As, None, lmbd, Q), zeros(kconv), 0.0, tol)
 
62
assert allclose(lmbd, lmbd_exact, tol*tol, 0.0)
 
63
 
 
64
print 'OK'
 
65
 
 
66
#-------------------------------------------------------------------------------
 
67
# Test 2: K = None
 
68
 
 
69
print 'Test 2',
 
70
 
 
71
lmbd_exact = zeros(ncv, 'd')
 
72
for k in xrange(ncv):
 
73
    lmbd_exact[k] =  A[k,k]/M[k,k]
 
74
 
 
75
 
 
76
X0 = RandomArray.random((n,ncv))
 
77
 
 
78
kconv, lmbd, Q, it, it_inner = jdsym.jdsym(As, Ms, None, ncv, 0.0, tol, 150, itsolvers.qmrs,
 
79
                                           jmin=5, jmax=10, eps_tr=1e-4, clvl=1)
 
80
 
 
81
assert ncv == kconv
 
82
assert allclose(computeResiduals(As, Ms, lmbd, Q), zeros(kconv), 0.0, normM*tol)
 
83
assert allclose(lmbd, lmbd_exact, normM*tol*tol, 0.0)
 
84
 
 
85
print 'OK'
 
86
 
 
87
#-------------------------------------------------------------------------------
 
88
# Test 3: general case
 
89
 
 
90
print 'Test 3',
 
91
 
 
92
lmbd_exact = zeros(ncv, 'd')
 
93
for k in xrange(ncv):
 
94
    lmbd_exact[k] =  A[k,k]/M[k,k]
 
95
 
 
96
kconv, lmbd, Q, it, it_inner = jdsym.jdsym(As, Ms, K, ncv, 0.0, tol, 150, itsolvers.qmrs,
 
97
                                           jmin=5, jmax=10, eps_tr=1e-4, clvl=1)
 
98
assert ncv == kconv
 
99
assert allclose(computeResiduals(As, Ms, lmbd, Q), zeros(kconv), 0.0, normM*tol)
 
100
assert allclose(lmbd, lmbd_exact, normM*tol*tol, 0.0)
 
101
 
 
102
print 'OK'
 
103
 
 
104
#-------------------------------------------------------------------------------
 
105
# Test 4: K = None, with X0
 
106
 
 
107
print 'Test 4',
 
108
 
 
109
lmbd_exact = zeros(ncv, 'd')
 
110
for k in xrange(ncv):
 
111
    lmbd_exact[k] =  A[k,k]/M[k,k]
 
112
 
 
113
# Fixme: RandomArray.random is broken AMD64
 
114
# X0 = RandomArray.random((n,ncv))
 
115
X0 = zeros((n,ncv), 'd')
 
116
for k in xrange(ncv):
 
117
    X0[k,k] = 10000
 
118
 
 
119
kconv, lmbd, Q, it, it_inner = jdsym.jdsym(As, Ms, None, ncv, 0.0, tol, 150, itsolvers.qmrs,
 
120
                                           jmin=5, jmax=10, eps_tr=1e-4, clvl=1, V0=X0)
 
121
 
 
122
assert ncv == kconv
 
123
assert allclose(computeResiduals(As, Ms, lmbd, Q), zeros(kconv), 0.0, normM*tol)
 
124
assert allclose(lmbd, lmbd_exact, normM*tol*tol, 0.0)
 
125
 
 
126
print 'OK'