%% LR-Zerlegung
% Von der  n\times n-MatrixA wird mit Gauß-Elimination mit Spaltenpivotisierung  eine Zerlegung PA=LR bestimmt.  Die Zwischenresultate können angezeigt werden, und die berechneten Matrizen P, L und R werden ausgegeben.
%% Die Beispielmatrizen
% Beispiel 3.24 (Buch)
% Reguläre Matrix, Pivotisierung notwendig (sonst Pivotelement gleich Null)
% Singuläre Matrix

bsp =3;
A=example_A(bsp)
n= length(A);

%% 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=0;

%% Erzeugung der Dreiecksgestalt
% Hier wird nun die Dreiecksgestalt von A erzeugt. Unterhalb der Diagonalen werden die L-Einträge abgespeichert.

piv= 1:n; % speichert die Reihenfolge der Zeilen
for i=1:n-1              
    [A,piv]=pivotisierung(A,piv,i,stepwise);
    if   A(i,i)~=0  % testet ob Diagonaleintrag ungleich Null ist
        A=annullierung(A,i,stepwise);
    end   
end

% Es folgt die Ausgabe des Ergebnisses der Zerlegung.

 disp('-------------------');
fprintf('\n\nEs folgen die berechneten Matrizen der Zerlegung:  \n \n')
L = eye(n) + tril(A,-1)
R = triu(A)
P = eye(n);
P = P(piv,:)
 disp('-------------------');
fprintf('Berechnete Zerlegung, (P^T)*L*R');
Resultat=P'*L*R  % Test: sollte die urspruengliche Matrix ergeben
fprintf('stimmt mit A überein?');

%
%% Aktivität:
% Bestimmen Sie die LR-Zerlegung für die drei Beispiele.  Interprerieren Sie die Ergebnisse.  Beachten Sie, dass das Verfahren auch für die singuläre Matrix (Beispiel 3) funktioniert.
%

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];
    end
end

function [A,piv] = pivotisierung(A,piv,i,stepwise)
     disp('-------------------');
% Die Pivotisierung, wie sie hier durchgeführt wird, sucht 
% in der Spalte vom behandelten Diagonalelement bis zur letzten 
% Zeile den betragsmäßig größten Eintrag    
        n=length(A);   
% Matlab liefert hier das Maximum und die Zeile des Maximums.        
        [amax,pos] = max(abs(A(i:n,i)));
        if amax==0
            fprintf('Pivotelement ist 0!\n');
            pos=1;
        end
% Der lokale Index pos liegt im Bereich 1...n-i+1 und muss
% in den globalen Index p im Bereich i...n umgerechnet werden.
        p = pos+i-1;
        if p>i
% Im Falle, dass das Diagonalelement nicht das größte Element ist,
% wird die Zeile mit dem behandelten Diagonalelement mit der Zeile 
% in der das betragsmäßige größte Element steht, vertauscht.            
           A([p i],:) = A([i p],:);
           piv([p i])= piv([i p]);
           if stepwise==1            
               fprintf('Zeile %d vertauscht mit Zeile %d',i,p)            
               A
           end
        end
end


function A = annullierung(A,i,stepwise)
    n=length(A);
    for k=i+1:n
        faktor = A(k,i)/A(i,i); % = L(k,i) aus LR-Zerlegung
        A(k,i)= faktor; % L-Eintrag unterhalb der Diagonalen speichern
        A(k,i+1:n) = A(k,i+1:n) - faktor*A(i,i+1:n);
        
        if stepwise==1
            fprintf('\nAnnulierung des Elements (%d,%d) durch L-Eintrag %f\n',k,i,faktor)
            A
        end
    end
end