% Musterlösung
%
% Numerische Mathematik, SS 2013
%
% LR_Zerlegung vollständige Pivotsuche (complete pivot)
%
% Ingolf Busch, Juni 2013

% Input: A     quadratische Matrix


% Output: S  struct mit allen Informationen 
%         S = lr(A)
%
%         L linke untere normierte Dreiecksmatrix
%         R rechte obere Dreiecksmatrix
%         P Permutationsmatrix für Saltenpivotsuche
%         Q zusätzliche Permutationsmatrix für vollständige Pivotsuche


function [L,R,P,Q] = my_lr_cp(A)

n=size(A,1);

pp = 1:n; qq = 1:n;

for k = 1:(n-1)
    
    %größtes Element bestimmen
    %erst alle Maxima in jeder spalte
    [colmaxima, rowindices] = max( abs(A(k:n, k:n)) );
    %dann das größte Spaltenmaximum
    [~, colindex] = max(colmaxima);
    %Bestimmen der Indizes des Pivotelementes
    row = rowindices(colindex)+k-1; col = colindex+k-1;
    
    % Permutiere das Pivotelement(samt zeilen und spalten) auf Hauptdiagonale
    A( [k, row], : ) = A( [row, k], : );
    A( :, [k, col] ) = A( :, [col, k] );
    %Merken der gemachten Permutation
    pp( [k, row] ) = pp( [row, k] );
    qq( [k, col] ) = qq( [col, k] );
    
    %Führe Elemination nach Schulbuchversion aus
    A((k+1):n,k) = A((k+1):n,k)/A(k,k);
    A((k+1):n,(k+1):n) = A((k+1):n,(k+1):n) - A((k+1):n,k) * A(k,(k+1):n);
end

if(nargout==1)
    L.mode = 'cp';
    L.L    = eye(n)+tril(A,-1);
    L.R    = triu(A);
    L.P    = eye(n);
    L.P    = L.P(pp,:);
    L.Q    = eye(n);
    L.Q    =L.Q(:,qq);
else
    L    = eye(n)+tril(A,-1);
    R    = triu(A);
    P    = eye(n);
    P    = P(pp,:);
    Q    = eye(n);
    Q    = Q(:,qq);
end