function Pieber
clc;

n = 20;
M = rand(n,n);
M = M.'*M;

% a) A -> Tridiagonale Matrix T
Q = eye(n);
A = M;
for k = 1:(n-2)
	u    = A((k+1):end, k);
	in   = norm(u, inf);
	u    = u/in;
	tau  = vorzeichen(u(1))*norm(u, 2);
	u(1) = u(1) + tau;
	P    = [eye(k), zeros(k,n-k); ...
           zeros(n-k,k), eye(n-k) - u*u.'/(tau*u(1))];
    
    A = P*A*P;
	Q = P*Q;
end

norm(M-Q'*A*Q)
input('?');

% Probe
T = Q*M*Q.';
%disp('Probe: ');
%disp(num2str(T));

% b) Eigenwerte
[eigenwerte, it] = QR_Algorithmus(T);

norm(sort(eigenwerte(:))-sort(eig(T) ))
input('?');

disp('Probe Eigenwerte:');
disp(['norm(eigenwerte-eig(A), 2): ',  ...
    num2str(norm(eigenwerte-eig(A), 2))]);
disp(['Iterationen: ', num2str(it)]);

return

% ----------------------------------------
% -----    vorzeichen    -----------------
% ----------------------------------------
function s = vorzeichen(x)
	if (x == 0)
		s = 1;
	else
		s = x/abs(x);
	end
return

% ----------------------------------------
% -----    QR_Algorithmus    -------------
% ----------------------------------------
function [a, it] = QR_Algorithmus(A)
it_max = 10000;
n = size(A, 1);
a = diag(A);
b = zeros(n-1,1);
for i=1:(n-1)
	b(i) = A(i+1,i);
end

x = zeros(n,1); 
y = zeros(n-1,1);
c = zeros(n-1,1);
s = zeros(n-1,1);
g = zeros(n-1,1);
d = zeros(n-1,1);
e = zeros(n-2,1);
it = 1;
while ((it < it_max) && (norm(b,1) > 1e-14))
	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);
    it = it + 1;
end
return
