Assig5part2

import numpy as np
import scipy
from scipy import linalg

def eig(A):
    k = 0
    #Use the householder reduction to hessenberg form
    Ak = tridiagonalization(A)
    m, n = np.shape(Ak)
    #If it is a 1x1 matrix return
    if(m < 2):
        return Ak
    #If the matrix is 2x2 run the algorithm for 5 times (gives good results)
    #If it is bigger run untill we get 0 element in sub/super diag
    max = 10000
    if m == 2:
        max = 5
    while k < max:
        #Shift
        mu = Ak[m-1][m-1]
        #QR factorisation 
        Q, R = scipy.linalg.qr(Ak-mu*np.eye(m))
        #Reconstruction of the matrix for the next step.
        Ak = R@Q+mu*np.eye(m)
        #Test if the matrix has zero element on the sub/super diag
        #If so call the function recursively 
        if m > 2:
            for i in range(1,m-1):
                if np.abs(Ak[i-1][i]) < 10**(-9):
                    A1 = eig(Ak[0:i,0:i])
                    A2 = eig(Ak[i:m,i:m])
                    return np.block([[A1, np.zeros((i, m-i))],[np.zeros((m-i,i)), A2]])
        k += 1
    return Ak

def tridiagonalization(A):
    m,n = np.shape(A)
    for k in range(0,m-2):
        x = A[k + 1 : m, k]
        size = np.size(x)
        e1 = np.eye(1,size,0)
        vk = np.sign(x[0])*scipy.linalg.norm(x)*e1 + x
        vk = vk/scipy.linalg.norm(vk)
        A[k+1:m,k:m] = A[k+1:m,k:m] - 2*vk.T@(vk@A[k+1:m,k:m])
        A[0:m,k+1:m] = A[0:m,k+1:m] - 2*(A[0:m,k+1:m]@vk.T)@vk
    return A

#Dim of A
n = 100
#Create a random symetric matrix
A = np.random.rand(n,n)
A = 1/2*(A+A.T)
#Call our eig to calculator
eigs = np.diag(eig(A))
eigs = np.sort(eigs)
print(eigs)
#Use scipy eig to calculate 
a,_ = scipy.linalg.eig(A)
a = np.real(a)
a = np.sort(a)
print(a)
print(scipy.linalg.norm(a-eigs))