function [L,U,P] = lu(A)
% Embedded MATLAB Library function.
%
% Limitations:
% 1) No NaN & Inf inputs.
% 2) Can only return three outputs.

% $INCLUDE(DOC) toolbox/eml/lib/matlab/lu.m $
% Copyright 2002-2004 The MathWorks, Inc.
% $Revision: 1.1.6.4 $  $Date: 2004/12/16 21:55:21 $

eml_assert(nargin > 0, 'error', 'Not enough input arguments.');
eml_assert(nargout == 3,'error',...
    'Embedded MATLAB only supports [L,U,P] = LU(A)');
eml_assert(isfloat(A), 'error', ['Function ''lu'' is not defined for values of class ''' class(A) '''.']);

[X,pivot] = eml_lu(A);
[L,U] = breakuplu(X);
P = getPermMatrix(pivot,class(X));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Breakup LU into separate L and U matrices
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [L,U] = breakuplu(X)
    m = size(X,1);
    n = size(X,2);

    if (isreal(X))
        if(m <= n)
            L = zeros(m, m, class(X));
            U = zeros(m, n, class(X));
        else
            L = zeros(m, n, class(X));
            U = zeros(n, n, class(X));
        end
    else
        if(m <= n)
            L = complex(zeros(m, m, class(X)));
            U = complex(zeros(m, n, class(X)));
        else
            L = complex(zeros(m, n, class(X)));
            U = complex(zeros(n, n, class(X)));
        end
    end

    for i = 1:m
        for j = 1:n
            if(i > j)
                L(i,j) = X(i,j);
            end
            if(i == j)
                L(i,j) = 1;
                U(i,j) = X(i,j);
            end
            if(i < j)
                U(i,j) = X(i,j);
            end
        end
    end

   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Convert pivot into permutation matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function p = getPermMatrix(pivot, cls)
    N = length(pivot);
    p = zeros(N,cls);
    for i = 1:N
        p(i,pivot(i)) = 1;
    end

    
function [X,pivot] = eml_lu(A)
% Embedded MATLAB Library function.
%
% Limitations:
% 1) No NaN & Inf inputs.

% $INCLUDE(DOC) toolbox/eml/lib/matlab/lu.m $
% Copyright 2002-2004 The MathWorks, Inc.
% $Revision: 1.1.6.4 $  $Date: 2004/12/16 21:55:21 $

eml_assert(nargin > 0, 'error', 'Not enough input arguments.');
eml_assert(isfloat(A), 'error', ['Function ''lu'' is not defined for values of class ''' class(A) '''.']);

m = size(A,1);
n = size(A,2);

if isempty(A) || m == 1
    X = A;
    pivot = 1:m;
    return;
end

% Init
X = A;

% Initialize row-pivot indices
pivot = 1:m;

for k = 1:min(m-1,n)
    p = k;    
    
    maxval = abs(X(k,k));
    % Scan the lower triangular part of this column only
    % MATLAB idiom:
    % [maxval,p] = max(abs(X(k:m,k)));
    % p = p+k-1;
    for i = k+1:m
        if(abs(X(i,k)) > maxval)
            maxval = abs(X(i,k));
            p = i;
        end
    end
    
    % Swap rows if max value if not in row k
    if(p ~= k)
        % MATLAB idiom:
        % X([p,k],:) = X([k,p],:); 
        for a = 1:n
            temp = X(p,a);
            X(p,a) = X(k,a);
            X(k,a) = temp;
        end
        % Swap pivot row indices
        % MATLAB idiom:
        % pivot([p,k]) = pivot([k,p]); 
        ptemp = pivot(k);
        pivot(k) = pivot(p);
        pivot(p) = ptemp;
    end

    % Divide lower triangular part of column by max
    Adiag = X(k,k);
    if Adiag ~= 0
        % MATLAB idiom:
        % X(k+1:m,k) = X(k+1:m,k) / X(k,k);
        for i = k+1:m
            X(i,k) = X(i,k) / Adiag;
        end
        % Subtract multiple of column from remaining columns
        % MATLAB idiom:
        % X(k+1:m,k+1:n) = X(k+1:m,k+1:n) - X(k+1:m,k)*X(k,k+1:n);
        for j = k+1:n
            for i = k+1:m
                X(i,j) = X(i,j) - X(i,k)*X(k,j);
            end
        end
    end
end
% [EOF]