%% Einführung in die numerische Mathematik

%% Blatt 5

%% Bsp_6

clear; clc;

%%nicht relevant
%%beliebige SPD Testmatrix A  wird erstellt (ähnlich wie bei Blatt 4, Bsp 6.)
M=2;
N=2;
N1=N+1;
N1N1=N1*N1;
h=1/N;
tau=1/M;
r=2;
d=3;
D=diag([4*ones(N1N1-1,1);3])-diag(ones(N1N1-1,1),1)-diag(ones(N1N1-1,1),-1)-diag(ones(N1N1-N1,1),N1)-diag(ones(N1N1-N1,1),-N1);
for k=N1:N1:N1N1-N1  
    D(k:k+1,k)=[3,0];
    D(k,k+1)=0;
    F(k+1,k)=0;    
end
B=(d/(h^2))*D+r*eye(N1N1);
A=eye(N1N1)+tau*B;

n=N1N1;


%% Bsp 6

%%-----------------------------------------------------------------------

%%Transformation der Matrix auf obere bidiagonale Form

%% A SPD Matrix, H Matrix zur Speicherung der bidiagonalen Form von A
%% n  Anzahl der Zeilen bzw. Spalten der Matrix A
%% Pm, Qm zur Speicherung der orthogonalen Transformationsmatrizen
%% ad, my, sigma, e Hilfsgrößen zur Housholdertransformation P=I-u*u'/sigma
%% c, tau2, sigma2, e Hilfsgrößen zur Householdertransformation
%% Q=I-u2*u2'/tau

Pm=eye(n);
Qm=eye(n);

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

H=A;
for j=1:n-1
    
ad=H(j:n,j);

my=norm(ad);

e=[1;zeros(n-j,1)];
 
u=ad+mysign(ad(1))*my*e;

sigma=(1/2)*(norm(u)^2);
 
u=[zeros(j-1,1);u];

P=eye(n)-(u*u')/sigma;

C=H(j:n,j+1:n);

c=P(j:n,j:n)*C;
c=c(1,:)';
tau2=norm(c);
e=[1;zeros(n-(j+1),1)];

u2=c+tau2*e*mysign(c(1));
sigma2=(1/2)*(norm(u2)^2);
u2=[zeros(j,1);u2];

Q=eye(n)-(u2*u2')/sigma2;

H=P*H*Q;

Pm=P*Pm;
Qm=Qm*Q;

end

norm(Pm*A*Qm-H)
input('?');

%Bidiagonal_A=H

B=H'*H;

%% Bestimmung der Singulärwerte von A SPD (Eigenwerte von A'*A)

%% A SPD Testmatrix, B Bidiagonale Form von A
%% n  Anzahl der Zeilen bzw. Spalten der Matrix A
%% weitere Bezeichungen analog zum Pseudo-Code der VO

for i=1:n
    
    a(i)=B(i,i);
    
end

for i=1:n-1
   
    b(i)=B(i,i+1);
end

kmax=100000;          %% maximale Anzahl der Iterationen(zur Absicherung)

    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)<10e-23
           break
       end
    end

    
    %% Ausgabe der Singulärwerte von A
    
%    Singulaerwerte_A=a
    %% Test u=eig(A'*A)
norm(sort(sqrt(a(:)))-sort(svd(A)))
