function Adanc

n=5;
A=randn(n,n);
I=eye(n,n);
A = A*A'+ eye(n);
ASym = true;

% Pseudo -Code zur Transformation auf Hessenberg oder
% Tridiagonale Form
H=A;
Q=I;

s0 = @(t) double(t>=0) - double(t<0);

for k=1:n-2
    u = H(k+(1:(n-k)),k);
    theta = norm(u,inf);
    u = u/theta;
    theta = s0(u(1))* norm(u,2);
    u(1) = u(1) + theta;
    theta = theta*u(1); 
    for l=k:n
        sm=0;
        for i=1:n-k
            sm=sm+H(k+i,l)*u(i);
        end
        for i=1:n-k
            H(k+i,l)=H(k+i,l)-sm*u(i)/theta;
        end
    end
    for l=1:n
        sm=0;
        for i=1:n-k
            sm=sm+Q(l,k+i)*u(i);
        end
        for i=1:n-k
            Q(l,k+i)=Q(l,k+i)-sm*u(i)/theta;
        end
    end
    if (ASym)
        for l=k:n
            sm=0;
            for i=1:n-k
                sm=sm+H(l,k+i)*u(i);
            end
            for i=1:n-k
                H(l,k+i)=H(l,k+i)-sm*u(i)/theta;
            end
        end
    end
end
norm(A-Q*H*Q')

tol=1e-8;
kmax=100;

a = zeros(n,1);
for i=1:n
    a(i)=H(i,i);
end
b = zeros(n-1,1);
for i=1:(n-1)
    b(i)=H(i,i+1);
end
x=zeros(n,1); g=x; y=zeros(n-1,1); d=y; e=zeros(n-2,1);
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)=c(j)*a(j+1)-s(j)*y(j);
        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) < tol*norm(a))
        break;
    end
end
[sort(a),sort(eig(A))]
norm(sort(a)-sort(eig(A)))
