# 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 = 50 epsilon = 0.01 steps = 100 def pcg(A, S, b): S_t = np.transpose(S) eig_values, eig_vectros = linalg.eig(A) M = max(eig_values) + 0.001 m = min(eig_values) - 0.001 x = [None for _ in xrange(0, steps)] r = [None for _ in xrange(0, steps)] z = [None for _ in xrange(0, steps)] d = [1 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]) E = np.zeros([n, n]) E[np.diag_indices_from(E)] = [1 for i in xrange(n)] z[k+1] = (2/(M-m)) * (d[k]/d[k+1])*((2*S - (m+M)*E)*z[k] + 2*r[k])- (d[k-1])/(d[k+1])*z[k-1] 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] d[k+1] = 2*(2-m-M)/(M-m)*d[k]-d[k-1] 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]) if k == 40: 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 pcg(A, S, b) print linalg.solve(A,b) if __name__ == '__main__': main()