function testTD
clear; clc;
%
% Transformiert Matrix auf Tridiagonalform
% mit Householder Reduktion
%
n = 5;
%disp('Ausgangsmatrix:');
A = rand(n,n); % Zufallsmatrix
A = A*A'+eye(n); % symmetrisch + SPD

%
H = A;
Q = eye(n); % Einheitsmatrix
s = @(t) double(t>=0) - double(t<0);

for k=1:n-2,
 %   [the, si, u] = householder(h)
 u = H(k+(1:(n-k)),k);
 t = norm(u,inf);
 u = u/t;
 t = s(u(1))* norm(u,2);
 u(1) = u(1) + t;
 t = t*u(1);    
 
    for l = k:n
       sm = dot(H(k+(1:n-k),l),u(1:n-k));
      H(k+(1:n-k),l)=H(k+(1:n-k),l) - sm*u((1:n-k))/t;
    end  
    for l= 1:n
        sm = dot(Q(l,k+1:n),u(1:n-k));
        Q(l,k+(1:n-k)) = Q(l,k+(1:n-k)) - sm*u(1:n-k)'/t;             
    end
    for l = k:n
        sm = dot(H(l,k+1:n), u(1:n-k));
        H(l,k+1:n) = H(l,k+1:n) - sm/t*u(1:n-k)';        
    end
    
end
H
norm(A-Q*H*Q')
[sort(eig(A)),sort(eig(H))]


% Eigenwerte von A

disp('Eigenwerte:');
eigenw1_A = eig(A);
eigenw2_A = myqr(H);

[sort(eigenw1_A),sort(eigenw2_A)]
norm(sort(eigenw1_A)-sort(eigenw2_A))



end

function ew = myqr(A)
n = size(A, 1);
a = diag(A);
b = diag(A,1);
tol = 1e-10;
k_max = 2000;
k = 0;
x = zeros(n,1);
y = zeros(n-1,1);
g = zeros(n,1);
c = zeros(n-1,1);
s = zeros(n-1,1);
d = zeros(n-1,1);
e = zeros(n-2,1);

while (k < k_max && norm(b,1) > tol) 
    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);        
k = k+1;      
       
end
ew = a;
end

