function A = gaussLR(A)
    % Gaussian Elimination for reforming the matrix into left-lower and right-upper triangular matrix
    % A: A R^(n \times n) 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)

    for i = 1:len
        prefactor = A(i,i)
        if prefactor == 0
            error("LR-Zerlegung failed. Please use pivotization, or check that the matrix is regular!")
        end

        % elimination (find R)
        factors = zeros(len, 1);
        for j = (i+1):len
            factors(j) = A(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(i,i:len))
        A(:,i:len) = A(:,i:len) + multmatrix;

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