ffffffffffffffffffffffffffffffffffff

function [L, U, P] = pivotedOuterProductLU(A)

    n = size(A, 1);
    U = A;              % We'll modify U in place
    L = eye(n);
    P = eye(n);

    for k = 1 to n-1
        
        % -- 1) Find the pivot in row k, from columns k..n:
        [~, pivotCol] = max( abs(U(k, k:n)) );
        pivotCol = pivotCol + (k - 1);

        % -- 2) Swap columns in U, P, and the part of L
        if pivotCol ~= k
            % Swap columns k and pivotCol in U
            tempCol = U(:, k);
            U(:, k) = U(:, pivotCol);
            U(:, pivotCol) = tempCol;
            
            % Swap columns k and pivotCol in P
            tempCol = P(:, k);
            P(:, k) = P(:, pivotCol);
            P(:, pivotCol) = tempCol;
            
            % For L, only swap columns k and pivotCol below row k
            % (strictly speaking, for i>k, L(i,k) ↔ L(i,pivotCol))
            if k > 1
                tempCol = L(:, k);
                L(:, k) = L(:, pivotCol);
                L(:, pivotCol) = tempCol;
            end
        end

        % -- 3) The pivot is now U(k,k). Use it to define multipliers in L
        pivotVal = U(k, k);
        for i = k+1 to n
            L(i, k) = U(i, k) / pivotVal;
        end

        % -- 4) Outer-product (rank-1) update of the submatrix
        for i = k+1 to n
            for j = k+1 to n
                U(i, j) = U(i, j) - L(i, k)*U(k, j);
            end
        end

    end

    % (At the end, U is upper-triangular, L is lower-triangular with 1s on diagonal.)
    % If desired, you can set L(i, i) = 1.0 explicitly for all i.
    
end