# coding utf-8 from numpy import dot as mul from precondition import imp

 ``` 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 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134``` ```# coding=utf-8 from numpy import dot as mul from precondition import * import numpy as np from numpy.random import randn as random_matrix from random import randint from numpy import linalg import scipy """ 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 steps = 10 def pcg(A, S, b): S_t = np.transpose(S) x = [None for _ in xrange(0, steps)] r = [None for _ in xrange(0, steps)] z = [None for _ in xrange(0, steps)] betta = [None for _ in xrange(0, steps)] p = [None for _ in xrange(0, steps)] alpha = [None for _ in xrange(0, steps)] x[0] = random_matrix(n, 1) r[0] = b - mul(A, x[0]) betta[1] = 0 for k in xrange(0, steps - 1): z[k] = mul(mul(S_t, S), r[k]) if k >= 1: betta[k+1] = sum((z[k][0] * r[k])[0]) / sum((z[k-1], r[k-1])[0]) p[k+1] = z[k] + betta[k+1][0] * p[k] else: p[1] = z[0] y = mul(A, p[k+1]) alpha[k+1] = sum(z[k] * r[k])[0] / sum(mul(A, p[k+1]) * p[k+1])[0] x[k+1] = x[k] + alpha[k+1] * p[k+1] r[k+1] = r[k] - alpha[k+1] * mul(A, p[k+1]) print sum(r[k+1] * r[k+1])[0] if sum(r[k+1] * r[k+1])[0] <= epsilon: print 'k: {d}'.format(k) return x[k+1] """ # 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() ```