function [U, L] = gaussian_elimination_with_complete_pivoting(A)

[n, m] = size(A);
if n ~= m
    error('Matrix is not squared!')
end
    
if det(A) == 0
    error('Matrix is not nonsingular!')
end

A

% for k = 1 : n-1
for k = 1 : n-1
    i = k:n;
    j = k:n;
    A(i, j)
    [max_val, rows_of_max_in_col] = max(abs(A(i, j)));
    [max_val, max_col] = max(max_val);
    max_row = rows_of_max_in_col(max_col);
    % Assing value of mi and lambda in respect to the main A matrix
    [mi, lm] = deal(max_row+k-1, max_col+k-1)
    A([k mi], 1:n) = deal(A([mi k], 1:n))
    A(1:n, [k lm]) = deal(A(1:n, [lm k]))
    p(k) = mi
    q(k) = lm
    % Perform Gaussian elimination with the greatest pivot
    
    if A(k, k) ~= 0
        rows = k+1 : n;
        A(rows, k) = A(rows, k)/A(k, k);
        A(rows, rows) = A(rows, rows) - A(rows, k) * A(k, rows);
    end
end

U = triu(A);
L = tril(A, -1) + eye(n);
p
I = eye(n);
P = I(p, :)
q
Q = I(:, q)