Bunch-kaufmann

 ``` 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``` ```import numpy as np from operator import itemgetter #-------------------------------------------------------- from math import factorial def binomial(n, k): if k < 0 or k > n: return 0 if k == 0 or k == n: return 1 return factorial(n) // (factorial(k) * factorial(n-k)) def hilb(n, m=0): if n < 1 or m < 0: raise ValueError("Matrix size must be one or greater") elif n == 1 and (m == 0 or m == 1): return np.array([[1]]) elif m == 0: m = n return 1. / (np.arange(1, n + 1) + np.arange(0, m)[:, np.newaxis]) def invhilb(n): H = np.empty((n, n)) for i in range(n): for j in range(n): H[i, j] = ((-1)**(i + j)) * (i + j + 1) * binomial(n + i, n - j - 1) * \ binomial(n + j, n - i - 1) * binomial(i + j, i) ** 2 return H #-------------------------------------------------------- def trans_matrix(size, idx1, idx2): permutation = np.array([[0 for i in xrange(size)] for j in xrange(size)]) for i in xrange(size): permutation[i][i] = 1 if idx1 != idx2: permutation[idx1][idx2], permutation[idx2][idx1] = 1 permutation[idx1][idx1], permutation[idx1][idx1] = 0 return np.matrix(permutation) def bunch_kaufmann(mtx, alpha): n_th = [] sum = 0 while sum < mtx.shape[0]: # i-th step of algorithm: we work with n[i]-th mtx_shiffer mtx_shiffer = mtx[sum: mtx.shape[0], sum: mtx.shape[1]] # find index with largest diagonal abs idx = max([(abs(mtx_shiffer[j][j]), j) for j in xrange(mtx_shiffer.shape[0])], key=itemgetter(0))[1] # make permutation matrix for (0, idx) transposition permutation = trans_matrix(mtx_shiffer.shape[0], 0, idx) # conjugate M' with permutation matrix mtx_shiffer = np.array(permutation*np.matrix(mtx_shiffer)*permutation) # find index for larger column abs and this abs [lambda_val, idx] = max([(abs(mtx_shiffer[j][0]), j) for j in xrange(mtx_shiffer.shape[0])], key=itemgetter(0)) n_k = 0 if abs(mtx_shiffer[0][0]) >= alpha*lambda_val: n_k = 1 # make identity permutation matrix permutation = trans_matrix(mtx_shiffer.shape[0], 0, 0) else: # same as above, for idx-th column of mtx_shiffer (without diagonal element) [sigma_val, j_idx] = max([(abs(mtx_shiffer[j][idx]), j) for j in xrange(mtx_shiffer.shape[0]) if j != idx], key=itemgetter(0)) if sigma_val*abs(mtx_shiffer[0][0]) >= alpha*lambda_val**2: n_k = 1 # make identity permutation matrix permutation = trans_matrix(mtx_shiffer.shape[0], 0, 0) else: if abs(mtx_shiffer[idx][idx]) >= alpha*sigma_val: n_k = 1 # make non-identity permutation matrix (0, idx) for mtx_shiffer permutation = trans_matrix(mtx_shiffer.shape[0], 0, idx) # conjugate M' with permutation matrix mtx_shiffer = permutation*np.matrix(mtx_shiffer)*permutation else: n_k = 2 permutation = trans_matrix(mtx_shiffer.shape[0], 2, idx)*trans_matrix(1, j_idx) # conjugate M' with permutation matrix mtx_shiffer = permutation*np.matrix(mtx_shiffer)*permutation sum += n_k n_th.append(n_k) return mtx #mtx = np.loadtxt('matrix.txt') mtx = hilb(10) f = open("result.txt", "w") f.write(str(mtx)) print bunch_kaufmann(mtx, 0.5) ```