function [A, P] = lrPivot(A)
    % Gaussian Elimination for reforming the matrix into left-lower and right-upper triangular matrix
    % A: A R^(n \times n) matrix
    % P: Permutations-Matrix

    % actually not used here (but it was planned to be used): bsxfun, which makes this possible: bsxfun(@times, [1,2;3,4], [1;2]) = [1,2;6,8]. It multiplies the first row by 1, and the second one by 2.
    % newer (and actually (still not) used here): [1,2;3,4].*[1;2] = [1,2;6,8]

    % for convenience
    len = length(A)
    pvec = 1:len


    % Generate LR-Zerlegung
    for i = 1:len

        % pivot
        [maxVal, maxValPos] = max(abs(A(pvec(i:end), i)))

        % swap
        temp = pvec(i)
        pvec(i) = pvec(maxValPos + i - 1)
        pvec(maxValPos + i - 1) = temp

        prefactor = A(pvec(i),i)
        if prefactor == 0
            error("LR-Zerlegung failed. Please check that the matrix is regular!")
        end

        % elimination (find R)
        factors = zeros(len, 1)
        for j = (i+1):len
            factors(pvec(j)) = A(pvec(j), i) / prefactor
        end

        % A(:,i:len) describes the matrix A with only the last i cols (so, for example, a 3x3 matrix with (:,2,len) with len=3) is the last 2 cols of the matrix, and with that a 2x3 matrix 
        
        % extra step for traceability, could also be compressed to a single line
        multmatrix = (-factors * A(pvec(i),i:len))
        A(:,i:len) = A(:,i:len) + multmatrix

        % save factors for L
        A(:,i) = A(:,i) + factors
    end

    P = eye(len)(:,pvec)
end