← 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_debug.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_debug.py
 
1
import math
 
2
import scipy
 
3
import spmatrix, itsolvers, jdsym, precon, superlu
 
4
 
 
5
test = 1
 
6
 
 
7
if test == 1:
 
8
    # Test: compare K=eye with K=None
 
9
    #
 
10
    #   results are not the same, ritz-value in 2nd iterations differ
 
11
 
 
12
    path = '/local/home/geus/matrices'
 
13
    A = spmatrix.ll_mat_from_mtx(path + 'edge6x3x5_A.mtx')
 
14
    M = spmatrix.ll_mat_from_mtx(path + 'edge6x3x5_B.mtx')
 
15
    n = A.shape[0]
 
16
    sigma = 25.0
 
17
 
 
18
    I = spmatrix.ll_mat(n, n)
 
19
    for i in xrange(n):
 
20
        I[i,i] = 1.0
 
21
    Keye = precon.jacobi(I)
 
22
 
 
23
    k_conv, lmbd, Q, it, it_inner = jdsym.jdsym(A.to_sss(), M.to_sss(), None, 5, sigma, 1e-10, 15, itsolvers.qmrs,
 
24
                                                jmin=5, jmax=10, eps_tr=1e-4, toldecay=2.0, linitmax=200, clvl=1, strategy=1)
 
25
 
 
26
    k_conv, lmbd, Q, it, it_inner  = jdsym.jdsym(A.to_sss(), M.to_sss(), Keye, 5, sigma, 1e-10, 15, itsolvers.qmrs,
 
27
                                                 jmin=5, jmax=10, eps_tr=1e-4, toldecay=2.0, linitmax=200, clvl=1, strategy=1)
 
28
 
 
29
elif test == 2:
 
30
 
 
31
    # Test 2: K = diag(A - sigma*M), test diagonal prec using Matlab
 
32
 
 
33
    Asigma = A.copy()
 
34
    Asigma.shift(-sigma, M)
 
35
    K = precon.jacobi(Asigma.to_sss())
 
36
 
 
37
    b = numpy.ones(n, 'd')
 
38
    x = numpy.zeros(n, 'd')
 
39
    K.precon(b, x)
 
40
    print 'norm(idiag) = %.16g' % (math.sqrt(numpy.dot(x, x)), )
 
41
 
 
42
    k_conv, lmbd, Q, it, it_inner  = jdsym.jdsym(A.to_sss(), M.to_sss(), K, 5, sigma, 1e-10, 150, itsolvers.qmrs,
 
43
                                       jmin=5, jmax=10, eps_tr=1e-4, toldecay=2.0, linitmax=200, clvl=1, strategy=1)
 
44
 
 
45
elif test == 3:
 
46
 
 
47
    Asigma = A.copy()
 
48
    Asigma.shift(-sigma, M)
 
49
    K = precon.ssor(Asigma.to_sss())
 
50
    k_conv, lmbd, Q, it, it_inner  = jdsym.jdsym(A.to_sss(), M.to_sss(), K, 5, sigma, 1e-10, 150, itsolvers.qmrs,
 
51
                                       jmin=5, jmax=10, eps_tr=1e-4, toldecay=2.0, linitmax=200, clvl=1, strategy=1)
 
52
 
 
53
elif test == 4:
 
54
 
 
55
    Asigma = A.copy()
 
56
    Asigma.shift(-sigma, M)
 
57
    K = precon.jacobi(Asigma.to_sss())
 
58
    k_conv, lmbd, Q, it, it_inner  = jdsym.jdsym(A.to_sss(), M.to_sss(), K, 5, sigma, 1e-10, 150, itsolvers.cgs,
 
59
                                       jmin=5, jmax=10, eps_tr=1e-4, toldecay=2.0, linitmax=200, clvl=1,
 
60
                                       strategy=1, optype=1)
 
61
 
 
62
elif test == 5:
 
63
 
 
64
    # time jdtest -e1 -k1 -o linsolver=QMRS,optype=SYM,jmin=5,jmax=10,linitmax=1000,jdtol=1e-6,strategy=1 1.4 1 ~/matrices/cop18_el5_A.mtx ~/matrices/cop18_el5_M.mtx
 
65
 
 
66
    path = '/home/geus/matrices/'
 
67
    A = spmatrix.ll_mat_from_mtx(path + 'cop18_el5_A.mtx')
 
68
    M = spmatrix.ll_mat_from_mtx(path + 'cop18_el5_M.mtx')
 
69
    n = A.shape[0]
 
70
    sigma = 1.4
 
71
 
 
72
    Asigma = A.copy()
 
73
    Asigma.shift(-sigma, M)
 
74
    K = precon.jacobi(Asigma.to_sss())
 
75
 
 
76
    k_conv, lmbd, Q, it, it_inner  = jdsym.jdsym(A.to_sss(), M.to_sss(), K, 1, sigma, 1e-6, 150, itsolvers.qmrs,
 
77
                                       jmin=5, jmax=10, eps_tr=1e-4, toldecay=2.0, linitmax=1000, clvl=1,
 
78
                                       strategy=1)
 
79
 
 
80
    print k_conv, lmbd, it, it_inner
 
81
 
 
82
elif test == 6:
 
83
 
 
84
    class prec2lev:
 
85
        def __init__(self, A, n11):
 
86
            n = A.shape[0]
 
87
            self.n11 = n11
 
88
            self.lu11 = superlu.factorize(A[:n11,:n11].to_csr(), diag_pivot_thresh=0)
 
89
            self.K22 = precon.ssor(A[n11:,n11:].to_sss())
 
90
            self.shape = (n, n)
 
91
        def precon(self, x, y):
 
92
            self.lu11.solve(x[:n11], y[:n11])
 
93
            self.K22.precon(x[n11:], y[n11:])
 
94
 
 
95
    print 'Loading matrices...'
 
96
    path = '/home/geus/matrices/'
 
97
    A = spmatrix.ll_mat_from_mtx(path + 'cop18_el5_A.mtx')
 
98
    M = spmatrix.ll_mat_from_mtx(path + 'cop18_el5_M.mtx')
 
99
    n = A.shape[0]
 
100
    sigma = 1.4
 
101
    n11 = 4688
 
102
 
 
103
    print 'Constructing preconditioner...'
 
104
    Asigma = A.copy()
 
105
    Asigma.shift(-sigma, M)
 
106
    K = prec2lev(Asigma, n11)
 
107
 
 
108
    k_conv, lmbd, Q, it, it_inner  = jdsym.jdsym(A.to_sss(), M.to_sss(), K, 2, sigma, 1e-8, 300, itsolvers.qmrs,
 
109
                                       jmin=10, jmax=25, eps_tr=1e-3, toldecay=1.5, linitmax=3000, clvl=1,
 
110
                                       strategy=1)
 
111
 
 
112
    print k_conv, lmbd, it, it_inner