← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Lib/sandbox/pysparse/examples/poisson_gmres.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/examples/poisson_gmres.py
1
 
import scipy
2
 
import math
3
 
import spmatrix
4
 
import itsolvers
5
 
import precon
6
 
import time
7
 
 
8
 
def poisson2d(n):
9
 
    L = spmatrix.ll_mat(n*n, n*n)
10
 
    for i in range(n):
11
 
        for j in range(n):
12
 
            k = i + n*j
13
 
            L[k,k] = 4
14
 
            if i > 0:
15
 
                L[k,k-1] = -1
16
 
            if i < n-1:
17
 
                L[k,k+1] = -1
18
 
            if j > 0:
19
 
                L[k,k-n] = -1
20
 
            if j < n-1:
21
 
                L[k,k+n] = -1
22
 
    return L
23
 
 
24
 
def poisson2d_sym(n):
25
 
    L = spmatrix.ll_mat_sym(n*n)
26
 
    for i in range(n):
27
 
        for j in range(n):
28
 
            k = i + n*j
29
 
            L[k,k] = 4
30
 
            if i > 0:
31
 
                L[k,k-1] = -1
32
 
            if j > 0:
33
 
                L[k,k-n] = -1
34
 
    return L
35
 
 
36
 
def poisson2d_sym_blk(n):
37
 
    L = spmatrix.ll_mat_sym(n*n)
38
 
    I = spmatrix.ll_mat_sym(n)
39
 
    P = spmatrix.ll_mat_sym(n)
40
 
    for i in range(n):
41
 
        I[i,i] = -1
42
 
    for i in range(n):
43
 
        P[i,i] = 4
44
 
        if i > 0: P[i,i-1] = -1
45
 
    for i in range(0, n*n, n):
46
 
        L[i:i+n,i:i+n] = P
47
 
        if i > 0: L[i:i+n,i-n:i] = I
48
 
    return L
49
 
 
50
 
n = 50
51
 
 
52
 
t1 = time.clock()
53
 
L = poisson2d(n)
54
 
print 'Time for constructing the matrix: %8.2f sec' % (time.clock() - t1, )
55
 
#L.export_mtx('poi2d_100.mtx')
56
 
 
57
 
A = L.to_csr()
58
 
 
59
 
 
60
 
S = L.to_sss()
61
 
print L.nnz
62
 
print S.nnz
63
 
print A.nnz
64
 
b = numpy.ones(n*n, 'd')
65
 
e = numpy.ones(n*n, 'd')
66
 
c = numpy.ones(n*n, 'd')
67
 
for loop in xrange(n*n):
68
 
    b[loop]= loop
69
 
    c[loop] = loop
70
 
y = numpy.ones(n*n, 'd')
71
 
S.matvec(b,y)
72
 
b = y
73
 
#print b
74
 
 
75
 
 
76
 
# ---------------------------------------------------------------------------------------
77
 
 
78
 
t1 = time.clock()
79
 
 
80
 
x = numpy.zeros(n*n, 'd')
81
 
info, iter, relres = itsolvers.gmres(S, b, x, 1e-12, 200, None, 100)
82
 
print 'info=%d, iter=%d, relres=%e' % (info, iter, relres)
83
 
 
84
 
print 'Time for solving the system using SSS matrix: %8.2f sec' % (time.clock() - t1, )
85
 
 
86
 
print 'norm(x) = %g' % math.sqrt(numpy.dot(x, x))
87
 
 
88
 
r = numpy.zeros(n*n, 'd')
89
 
S.matvec(x, r)
90
 
r = b - r
91
 
print 'norm(b - A*x) = %g' % math.sqrt(numpy.dot(r, r))
92
 
 
93
 
# ---------------------------------------------------------------------------------------
94
 
 
95
 
t1 = time.clock()
96
 
 
97
 
x = numpy.zeros(n*n, 'd')
98
 
info, iter, relres = itsolvers.gmres(A, b, x, 1e-12, 200)
99
 
print 'info=%d, iter=%d, relres=%e' % (info, iter, relres)
100
 
 
101
 
print 'Time for solving the system using CSR matrix: %8.2f sec' % (time.clock() - t1, )
102
 
 
103
 
print 'norm(x) = %g' % math.sqrt(numpy.dot(x, x))
104
 
 
105
 
r = numpy.zeros(n*n, 'd')
106
 
A.matvec(x, r)
107
 
r = b - r
108
 
print 'norm(b - A*x) = %g' % math.sqrt(numpy.dot(r, r))
109
 
 
110
 
# ---------------------------------------------------------------------------------------
111
 
 
112
 
t1 = time.clock()
113
 
 
114
 
x = numpy.zeros(n*n, 'd')
115
 
info, iter, relres = itsolvers.gmres(L, b, x, 1e-12, 200)
116
 
print 'info=%d, iter=%d, relres=%e' % (info, iter, relres)
117
 
 
118
 
print 'Time for solving the system using LL matrix: %8.2f sec' % (time.clock() - t1, )
119
 
 
120
 
print 'norm(x) = %g' % math.sqrt(numpy.dot(x, x))
121
 
 
122
 
r = numpy.zeros(n*n, 'd')
123
 
A.matvec(x, r)
124
 
r = b - r
125
 
print 'norm(b - A*x) = %g' % math.sqrt(numpy.dot(r, r))
126
 
# ---------------------------------------------------------------------------------------
127
 
 
128
 
K_ssor = precon.ssor(S, 1.0)
129
 
t1 = time.clock()
130
 
 
131
 
x = numpy.zeros(n*n, 'd')
132
 
info, iter, relres = itsolvers.gmres(S, b, x, 1e-12, 500, K_ssor, 20)
133
 
print 'info=%d, iter=%d, relres=%e' % (info, iter, relres)
134
 
 
135
 
print 'Time for solving the system using SSS matrix and SSOR preconditioner: %8.2f sec' % (time.clock() - t1, )
136
 
 
137
 
print 'norm(x) = %g' % math.sqrt(numpy.dot(x, x))
138
 
 
139
 
r = numpy.zeros(n*n, 'd')
140
 
S.matvec(x, r)
141
 
r = b - r
142
 
print 'norm(b - A*x) = %g' % math.sqrt(numpy.dot(r, r))
143
 
 
144
 
# ---------------------------------------------------------------------------------------
145
 
 
146
 
#import jdsym
147
 
#jdsym.jdsym(S, None, None, 5, 0.0, 1e-8, 20, itsolvers.qmrs, clvl=1)
148
 
 
149
 
x = numpy.zeros(n*n, 'd')
150
 
info, iter, relres = itsolvers.gmres(S, b, x, 1e-15, 500, K_ssor, 50)
151
 
print 'info=%d, iter=%d, relres=%e' % (info, iter, relres)
152
 
 
153
 
print 'Time for solving the system using SSS matrix and SSOR preconditioner: %8.2f sec' % (time.clock() - t1, )
154
 
 
155
 
print 'norm(x) = %g' % math.sqrt(numpy.dot(x, x))
156
 
 
157
 
r = numpy.zeros(n*n, 'd')
158
 
S.matvec(x, r)
159
 
r = b - r
160
 
print 'norm(b - A*x) = %g' % math.sqrt(numpy.dot(r, r))
161
 
print 'bye'