coding utf-8 from numpy import dot as mul from precondition import lab

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
# coding=utf-8
from numpy import dot as mul
from precondition import *
"""
2 lab
=================
PCG
полиномиальное ускорение
фиксируется число шагов
MZ = R
трёхчленная реккурентная формула(через полиномы чебышёва и лежандра)
http://www.wikiwand.com/en/Symmetric_successive_overrelaxation ssor
https://github.com/pyamg/pyamg/tree/master/pyamg/relaxation тут и Чёбышев есть
"""
n = 10
epsilon = 0.01
steps = 1000000000
def pcg(A, S, b):
S_t = np.transpose(S)
# array: [k-1, k, k+1]
x = [None, random_matrix(n, 1), None]
r = [b - mul(A, x[1]), b - mul(A, x[1]), None]
z = [mul(mul(S_t, S), r[1]), mul(mul(S_t, S), r[1]), None]
betta = [None, 0, 0]
p = [None, z[1], None]
alpha = [None, None, None]
for k in xrange(1, steps):
z[1] = mul(mul(S_t, S), r[1])
betta[2] = sum(z[1] * r[1])[0] / sum(z[0], r[0])[0]
p[2] = z[1] + betta[2][0] * p[1] # TODO error
alpha[2] = sum(z[1] * r[1])[0] / sum(mul(A, p[2]) * p[2])[0]
x[2] = x[1] + alpha[2] * p[2]
r[2] = r[1] - alpha[2] * mul(A, p[2])
if sum(r[2] * r[2])[0] <= epsilon:
return x[2]
# тут массивы сдвигаем на 1 влево:)
x[:-1], r[:-1], z[:-1], betta[:-1], p[:-1], alpha[:-1] = x[1:], r[1:], z[1:], betta[1:], p[1:], alpha[1:]
# [:-1] [1:]
"""
x_k = random_matrix(n, 1)
r_k = b - mul(A, x_k)
betta = 0
k = 0
z_k = mul(mul(S_t, S), r_k)
while sum(r_k * r_k) >= epsilon:
z_1, r_1 = z_k, r_k
z_k = mul(mul(S_t, S), r_k)
if k >= 1:
betta = sum(z_k * r_k) / sum(z_1, r_1)
p_k = z_k + betta * p_k
else:
p_k = z_1
p_k = z_k + betta * p_k
alpha = sum(z_k * r_k) / sum(mul(A, p_k), p_k)
x_k = x_k + alpha * p_k
r_k = r_k - alpha * mul(A, p_k)
k += 1
"""
def main():
#A, S = jacobi_precondition(n=n)
A, S = ssor_precondition(w=1.0, n=n)
b = random_matrix(n, 1)
print pcg(A, S, b)
if __name__ == '__main__':
main()