function [P,L,R] = LR_mit_pivot(A)
%
%  Zerlegung einer Matrix A in die Form PA=LR
%
%  Input
%     A   invertierbare quadratische Matrix
%
%  Output
%     P   Permutationsmatrix
%     L   linke untere Dreiecksmatrix
%     R   rechte obere Dreiecksmatrix
%           so dass PA=LR
%
%

  n =     size(A,1);
  L =     eye(n);
  P =     eye(n);
  pivot = 1:n;  % enthaelt die Nummern der Zeilen der urspruengl. Matrix A
                % in der Reihenfolge der aktuellen Permutation
 
  for k = 1 : n   % fuer alle Spalten:
    [m,mind] = max(abs(A(k:n,k)));  % mind liefert den Zeilen-Index 
                                    % des Pivotelements
    if mind ~= 1  % Spalten werden vertauscht, belege Vektor pivot neu
       pivot([k mind+k-1])=pivot([mind+k-1 k]);
    end
     
    for i = k+1 : n    % Gauss-Elimination
      l             = A(pivot(i),k) / A(pivot(k),k);
      A(pivot(i),k) = l;     % Speichern der Eintraege von L unterhalb der Hauptdiagonalen
      A(pivot(i),k+1:n) = A(pivot(i),k+1:n) - l * A(pivot(k),k+1:n); 
    end
  end
  
  R = triu(A(pivot,1:n));
  L = A(pivot,1:n) - R + eye(n);
  P = P(pivot,1:n);
