clear
clc
%----------- Einf. i. d. Numerik Übungsblatt 5 -----------
% Marlene Lepuschitz

%--------------a)---------------
n=5;
A=rand(n);
%Erzeugen von SPD Matrix
A=A*A'+eye(n);

Asav = A;
Psav = eye(n);

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


%Tridiagonalform
alpha = -mysign(A(2,1))*sqrt(sum(A(2:n,1).^2));
r     = sqrt((1/2)*(alpha^2-A(2,1)*alpha));
v1    = zeros(n);
v1(2) = (A(2,1)-alpha)/(2*r);
for k=3:n
    v1(k)=(A(k,1))/(2*r); 
end

P1=eye(n)-2*v1*v1';
A=P1*A*P1;

Psav = P1*Psav;

for k=2:n-2
    alpha = -mysign(A(k+1,k))*sqrt(sum(A(k+1:n,k).^2));
    r     = sqrt((1/2)*(alpha^2-A(k+1,k)*alpha));
    v     = zeros(n);
    K     = k+1;
    v(K)  =(A(k+1,k)-alpha)/(2*r);
    for j=k+2:n
        v(j)=A(j,k)/(2*r);
    end
    P=eye(n)-2*v*v';
Psav = P*Psav;
    A=P*A*P;
end
A_trid=A

norm(A_trid-Psav*Asav*Psav')
input('?');

%----------b)------------------
%Eigenwerte
n     = size(A,1);
a     = diag(A);
b     = diag(A,1);
tol   = 1e-5;
k_max = 350;

%Speicherplatz vorallokieren
x = zeros(n,1);
y = zeros(n-1,1);
r = zeros(n,1);
c = zeros(n-1,1);
s = zeros(n-1,1);
d = zeros(n-1,1);
e = zeros(n-2,1);
k = 0;
 
while (k < k_max && norm(b,1) > tol) 
    x(1) = a(1);
    y(1) = b(1);
     
    for j = 1:(n-1)
        r(j) = sqrt(x(j)^2 + b(j)^2);
        c(j) = x(j)/r(j);
        s(j) = b(j)/r(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
    
    r(n) = x(n);
    a(1) = s(1) * d(1) + c(1)*r(1);
    b(1) = s(1) * r(2);
    
    for j = 2: (n-1)
        a(j) = s(j)*d(j)+c(j)*c(j-1)*r(j);
        b(j) = s(j)*r(j+1);
    end
  
    a(n) = c(n-1)*r(n);  
    k = k+1;      
end
Eigenwerte = a

Test=eig(Asav)
