lu solver using scaled partial pivoting

lu_spp.m

function [L, U, P] = lu_spp(A)
% [L, U, P] = lu_spp(A) computes a unit lower triangular matrix L,
% an upper triangular matrix U, and a permutation matrix P such that
% P * A = L * U.
% The LU factorization is computed using scaled partial pivoting.

% for the n x m matrix A, n == m
size_A = size(A);
if(size_A(1) ~= size_A(2))
    return
end
n = size_A(1);

s = max(abs(A), [], 2); % the scalar coefficient vector

L = eye(n);
P = eye(n);
% iterate row numbers k = 1, ..., n - 1
for k = 1 : n - 1
    % find the max abs value divided by the corresponding scalar
    % coefficient in column k
    [~, j] = max(abs(A(k : n, k)) ./ s(k : n));
    % normalize the index
    j = j + k - 1;
    % swap rows if the index of the max abs value divided by the
    % corresponding scalar coefficient is not k
    if (j > k)
        % swap rows in the itermediate matrix A
        A([j, k], :) = A([k, j], :);
        % swap rows in the permutation matrix P
        P([j, k], :) = P([k, j], :);
        % swap elements in the lower triangular matrix L
        if (k > 1) % important
            L([j, k], 1 : k - 1) = L([k, j], 1 : k - 1);
        end
    end
    % assign the scalar coefficient to the lower triangular matrix L
    L(k + 1 : n, k) = A(k + 1 : n, k) / A(k, k);
    % subtract the scalar coefficient, multiplied by the pivoting
    % row k, from the current row i
    A(k + 1 : n, k : n) = A(k + 1 : n, k : n) - L(k + 1 : n, k) * A(k, k : n);
end
% assign the itermediate matrix A to the upper triangular matrix U
U = A;