%% Beispiel 6 Blatt 5
clc; clear; 

N = 1;
N1 = N+1;
n = N1*N1;
D0 = -ones(n,1); 
A = spdiags(D0,-N1,n,n) + spdiags(D0,N1,n,n) + spdiags(D0,1,n,n) ...
    + spdiags(D0,-1,n,n) + speye(n)*4; D(n,n) = 3;

%% a) Householder Transformation
H = A;
Q = eye(size(A));
for k = 1:n-1
    s = 0;
    for i = 1:n-k
        s = s + H(k+i,k)^2;
    end
alpha = -sign(H(k+1,k))*sqrt(s);
r = sqrt(0.5*(alpha^2-H(k+1,k)*alpha));
v = zeros(n,1);
v(k+1) = (H(k+1,k)-alpha)/(2*r);
v(k+2:n)= H(k+2:n,k)/(2*r);
P = eye(n) - 2*v*(v');
H = P*H*P;
Q = Q*P;
end

disp(strcat('||A-Q*H*Q''||=',num2str(norm(A-Q*H*Q'))))

%% b) Bestimmung der Eigenwerte (QR-Algorithmus)

a = diag(H);
b = zeros(1,n);
x = b; y = b; g = b; c = b; d = b; s = b; e = b;
for i = 1:n-1
    b(i) = H(i,i+1);
end
kmax = 100;
eps = 1*e-3;
for k = 1:kmax;
    x(1) = a(1); y(1) = b(1);
    for j = 1:n-1
        g(j) = sqrt(x(j)^2 + b(j)^2);
        c(j) = x(j)/g(j); s(j) = b(j)/g(j);
        d(j) = c(j)*y(j) + s(j)*a(j+1);
        x(j+1) = -s(j) * y(j) + c(j)*a(j+1);
        if j ~= n-1
            e(j) = s(j)*b(j+1); y(j+1) = c(j)*b(j+1);
        end
    end
    g(n) = x(n);
    a(1) = s(1)*d(1)+c(1)*g(1); b(1)= s(1)*g(2);
    for j = 2:n-1
        a(j) = s(j)*d(j) + c(j)*c(j-1)*g(j);
        b(j) = s(j) * g(j+1);
    end
    a(n) = c(n-1)*g(n);
    if(norm(b)< eps) 
        break
    end
end

disp(strcat('||a-eig(A)||=',num2str(norm(sort(a(:))-sort(eig(A))))))
