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
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
# 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 = 5
epsilon = 0.1
import numpy as np
from numpy.random import randn as random_matrix
from random import randint
from numpy import linalg
import scipy
def jacobi_precondition(n=10):
A = random_matrix(n, n)
D = np.tril(np.triu(A, 0))
return A, D
steps = 10
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])
print sum(r[2] * r[2])[0]
if sum(r[2] * r[2])[0] <= epsilon:
print 'k: {d}'.format(k)
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 linalg.solve(A,b)
print pcg(A, S, b)
if __name__ == '__main__':
main()