%% QR-Zerlegung
% In diesem Matlab Programm wird die schrittweise Durchführung der Givens- und Householder-Verfahren zur Bestimmung einer QR-Zerlegung
%         QA=R   (also A=Q^TR)
% einer Matrix A \in \mathbb{R}^{m\times n} illustriert.
%% Die Beispielmatrizen:
% Beispiel 3.24 (Buch)   (m=n)
% Matrix, die bei Gauss-Elimination Pivotisierung erfordert  (m=n)
% Singuläre Matrix   (m=n)
% Matrix  mit m>n
% Matrix mit m<n
% Sparse-Matrix mit voller erster Spalte
% Sparse-Matrix

bsp = 4
format short;
A=example_A(bsp)

% Orthogonale Transformationen
% Sie können als orthogonale Transformation eine der beiden Methoden (Givens- oder Householderverfahren) wählen:
% Givens-Rotationen
% Householder-Spiegelungen

method=1;

% Schrittweiser Durchlauf
% Wenn das Verfahren nicht schrittweise durchlaufen werden soll, ist der folgende Parameter stepwise auf 0 zu setzen. Ansonsten ist der Wert 1 zuzuweisen.

stepwise=1;

%% Die QR-Zerlegung
% Im folgenden wird nun die QR-Zerlegung der Matrix A berechnet  (QA=R).

[rowsA,columnsA] = size(A);
Q = eye(rowsA,rowsA);

if (rowsA > columnsA) 
    jstop=columnsA; 
else 
    jstop=rowsA-1; 
end

switch method
    case 1 % Givens
        for j=1:jstop
           disp('-------------------'); 
            for i=j+1:rowsA
                [Q,A]=givens_schritt(A,Q,j,i,stepwise);
            end
        end
    case 2 % Householder
    for j=1:jstop
        [Q,A]=householder_schritt(A,Q,j,stepwise);
    end
end    

% Es folgt die Ausgabe des Ergebnisses der Zerlegung. Zur Kontrolle wird noch die Fehlermatrix E=A-Q^TR bestimmt.

 disp('-------------------'); 
fprintf('\n\nEs folgen die berechneten Matrizen der Zerlegung:  \n \n')
R=A
Q
 disp('-------------------'); 
fprintf('\n Die Fehlermatrix A-(Q^T)*R wird berechnet:   \n');
Fehler = example_A(bsp) - Q'*R  

%
%% Aktivität:
% Wenden Sie beide Methoden auf die verschiedenen Beispielmatrizen an.
%

function A = example_A(k)  
    switch k
        case 1
            A =[2 -1 -3 3; 4 0 -3 1; 6 1 -1 6; -2 -5 4  1];
        case 2
            A = [1 2 3 4; -2 1 0 -1;1 -3 -3 5; 0 10 21 24];
        case 3
            A = [1 3 8 4;2 6 22 15;-1 -3 1 -2;3 9 29 22];
        case 4
            A =[2 -1 -3; 4 0 -3; 6 1 -1; -2 -5 4; 2 8 -1];
        case 5
            A =[2 -1 -3 3 -2; 4 0 -3 1 5; 6 1 -1 6 1];
        case 6
            A = eye(8,4);
            A(2,1)=1;                        
            A(3,1)=-2;
            A(4,1)=0.5;                        
            A(5,1)=0.2;            
            A(6,1)=-0.5;            
            A(7,1)=1;
            A(8,1)=-2;
            A(2,4)= 3;
            A(6,3)=-6;                        
            A(3,4)=-6;
            A(7,1)=2;
            A(7,3)=10;               
        case 7
            A = eye(8,3);
            A(4,1)= 3;
            A(6,2)=-6;
            A(1,3)=2;
            A(7,3)=10;
    end
end


function [Q,A]=givens_schritt(A,Q,j,i,stepwise)
    if abs(A(i,j)) >= 10^(-15)  % nur Eintraege ``ungleich Null'' werden eliminiert
        [rowsA,columnsA]=size(A);
        Qs=eye(2);
        r=sqrt(A(j,j)*A(j,j)+A(i,j)*A(i,j));
        c=A(j,j)/r;
        s=A(i,j)/r;
        Qs=[c s; -s c]; % Givens-Matrix
        A([j i],j:columnsA) = Qs*A([j i],j:columnsA);
        Q([j i],1:rowsA) =Qs*Q([j i],1:rowsA);
        if stepwise == 1
            fprintf('\nAnnulierung des Elements A(%d,%d)',i,j);
            A
        end
    end
end

function [Q,A]=householder_schritt(A,Q,j,stepwise)
        [rowsA,columnsA]=size(A);
        alpha=norm(A(j:rowsA,j),2) * sign(A(j,j));
        v=A(j:rowsA,j);
       if norm(v) >= 10^(-15) %teste of v ``ungleich Null'' ist
        v(1)=v(1)+alpha;
        vs=2/(v'*v);        
        Qv=eye(length(v))-vs*(v*v'); % Householder-Matrix
        % Beachte: Householder-Matrix wird in der Praxis nicht
        % berechnet, sondern auf Vektoren angewandt (Rang-1-Update)
        A(j:rowsA,j:columnsA) = Qv*A(j:rowsA,j:columnsA);
        Q(j:rowsA,1:rowsA)=Qv*Q(j:rowsA,1:rowsA);
        if stepwise == 1
            fprintf('\nAnnulierung der Elemente A(%d:%d,%d)',j+1,rowsA,j);
            A
        end
       end 
    end
