function blatt5_hrassnigg

clc; clear;

%% Achtung Matrix eingefügt:

N=5;
M=N;
N1 = N+1;
N1N1=N1^2;
h=1/N;
tau=1/M;

d = 1;
r = 1;


aa=linspace(4,4,(N+1)^2);
for i=1:N+1
    aa(i*(N+1))=3;
end

D=diag(aa)-diag(ones((N+1)^2-1,1),1)-diag(ones((N+1)^2-1,1),-1)...
     -diag(ones((N+1)^2-N-1,1),N+1) -diag(ones((N+1)^2-N-1,1),-N-1);
 
 for i=linspace(N1,N1N1-N1,N1-1)
     D(i+1,i)=0;
     D(i,i+1)=0;
     
 end
 

I2=eye(N1N1);

B=d/h^2*D+r*I2;
A=(I2+tau*B);
n=size(A,1);

norm(A-A');

% Ende Matrix eingefügt

%% Auf Bidiagonalform bringen:
B=A;
for k=1:n
    
    x=B(k:n,k);
    e    =  zeros(size(x));
    e(1) = 1;
    u=x+sign(x(1))*norm(x,2)*e;
    u=u/norm(u,2);
    B(k:n,k:n)=B(k:n,k:n)-2*u*(u'*B(k:n,k:n));
    
    if k<=n-2
        x=B(k,k+1:n);
        x=x';
        e=zeros(size(x));
        e(1)=1;
        v=x+sign(x(1))*norm(x,2)*e;
        v=v/norm(v,2);
        B(k:n,k+1:n)=B(k:n,k+1:n)-2*(B(k:n,k+1:n)*v)*v';
    end
end



%% Singulärwertzerlegung

Z=B'*B;

EigZ=QRzerleg(Z);

EigZ=sqrt(EigZ);

%% Zusammenfassung / Ausgabe:

disp([num2str(EigZ')  '   =   Singulärwerte von A mit QR-Zerlegung berechnet'])
disp([num2str(svd(A)') '  =   Singulärwerte von A mit der in Matlab implementierten Funktion berechnet'])

end




%% Funktionen

function [Q,H]=tridiag(A)


n=size(A,1);
I=eye(n);
H=A;
Q=I;


 
for k=1:(n-2)
    u = zeros(n,1);
    
    for i=1:(n-k)
        u(i)=H(k+i,k);
    end

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

function [a]=QRzerleg(H)
n=size(H,1);
A  = H;
a  = zeros(n,1);
b  = zeros(n-1,1);
ep = 1.0e-10;
for i=1:n
    a(i)=A(i,i);
end

for i=1:n-1
    b(i)=A(i,i+1);
end
    kmax=1000;
    x  =  zeros(n,1);
    g  =  x;
    c  =  x;
    d  =  x;
    s  =  x;
    e  =  x;
    y  =  x;
    
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)<ep
        break
    end
end
end
