function [LR P] = lrpivot(A)
% function [LR P] = lr(A)
% berechnet die LR-Zerlegung der Matrix A mit Spaltenpivotsuche
%
% Eingabe:
% A...quadratische Matrix
%
% Ausgabe:
% LR...quadratische Matrix.
%      - Unterhalb der Diagonalen stehen die Eintraege von L
%        von der Diagonalen aufwaerts stehen die Eintraege von R
% P....Permutationsvektor. 
%      - Der i-te Eintrag in P gibt an, mit welcher Zeile
%        die urspruenglich i-te Zeile in A vertauscht worden ist.
%
% Es gilt:
% A(P,:) = L*R mit
%               L=eye(size(A))+tril(LR,-1)
%        und    R=triu(LR);
%
% Markus Grasmair,
% Wien, April 2011

% Teste, ob die Eingabematrix quadratisch ist
[m n] = size(A);
if(n~=m)
    error('Die Eingabematrix ist nicht quadratisch.');
end    
    
% Initialisierung
LR = A;
% Der Permutationsvektor wird als Integervektor abgespeichert
P = uint16(1:n);

% Iteration
for i=1:(n-1)
    % Pivotsuche:
    % bringe die Zeile an die erste Stelle, fuer die der i-te Eintrag im
    % Vergleich zur ganzen Zeile relativ am groessten ist.

    temp = abs(LR(i:end,i))./sum(abs(LR(i:end,i:end)),2);
    % Wenn eine Zeile in R Null wird, dann ist die Matrix nicht
    % invertierbar. Dies ist genau dann der Fall, wenn in temp ein Eintrag
    % der Form 0/0 vorkommt - diese werden in Matlab als NaN (not a number)
    % abgespeichert
    if(any(isnan(temp)))
        % Teste ob irgendein Eintrag von temp gleich NaN ist.
        % Wenn ja, bricht der Algorithmus ab.
        error('Die Eingabematrix ist nicht invertierbar.');
    else
        % Die erste Ausgabe der Funktion max gibt das Maximum an,
        % die zweite, an welcher Position das Maximum angenommen wurde
        % Wir brauchen eigentlich nur die Position; um auf diese zugreifen
        % zu koennen, muessen wir aber auch den Wert des Maximums auslesen,
        % den wir in die (nie benoetigte) Variable temp2 schreiben
        [temp2 pivot] = max(temp);
    end
    % Da wir fuer die Bestimmung des Pivotelements nur die Untermatrix der
    % letzten (n+1-i) Zeilen von R verwendet haben, muessen wir zur
    % Bestimmung der tatsaechlich Position noch i addieren.
    pivot=uint16(pivot+i-1);
    
    % Wenn das Pivotelement nicht in der aktuellen Zeile steht,
    % vertausche die Zeilen in R und berechne den neuen Permutationsvektor
    % durch Vertauschen der Eintraege in P
    if(pivot~=i)
        LR([i pivot],:) = LR([pivot i],:);
        P([i pivot]) = P([pivot i]);
    end
    
    % LR Update:
    % Dividiere die i-te Spalte unterhalb der Diagonalen durch das
    % Diagonalelement
    LR((i+1):n,i) = LR((i+1):n,i)/LR(i,i);
    % Fuer j,k=i+1,...,n setze
    % LR(j,k) = LR(j,k) - LR(j,i)*LR(i,k)
    LR((i+1):n,(i+1):n) = LR((i+1):n,(i+1):n)-LR((i+1):n,i)*LR(i,(i+1):n);
end