function [L,U,P] = lu(A,opt)
%MATLAB Code Generation Library Function

%   Copyright 2002-2018 The MathWorks, Inc.
%#codegen

coder.columnMajor;
if nargout == 3
    coder.inline('always'); % For shrinktype in the vector p case.
end
coder.internal.assert(nargin > 0,'MATLAB:minrhs');
coder.internal.assert(isa(A,'float'), ...
    'Coder:toolbox:unsupportedClass','lu',class(A));
coder.internal.assert(ismatrix(A),'Coder:toolbox:lu_inputMustBe2D');
if nargin == 2
    coder.internal.prefer_const(opt);
    coder.internal.assert(coder.internal.isTextRow(opt), ...
        'MATLAB:factorFlags:invalidPermFlag');
    len = max(1,strlength(opt));
    VECTOR = strncmpi(opt,'vector',len);
    coder.internal.assert(VECTOR || strncmpi(opt,'matrix',len), ...
        'MATLAB:factorFlags:invalidPermFlag');        
else
    VECTOR = false;
end
coder.internal.assert(nargin == 1 || VECTOR || strcmpi(opt,'matrix'), ...
    'MATLAB:factorFlags:invalidPermFlag', ...
    'IfNotConst','Fail');
if nargout <= 1
    L = LU(A);
elseif nargout == 2
    [L,U] = LU(A);
elseif VECTOR
    [L,U,iP] = LUp(A);
    P = double(iP);
else
    [L,U,P] = LUP(A);
end

%--------------------------------------------------------------------------

function [L,U] = LU(A)
ONE = coder.internal.indexInt(1);
m = coder.internal.indexInt(size(A,1));
n = coder.internal.indexInt(size(A,2));
if nargout == 1
    L = coder.internal.lapack.xgetrf(m,n,A,ONE,m);
else
    [A,ipiv] = coder.internal.lapack.xgetrf(m,n,A,ONE,m);
    [L,U] = expandlu(A,eml_ipiv2perm(ipiv,m));
end

%--------------------------------------------------------------------------

function [L,U,P] = LUp(A)
ONE = coder.internal.indexInt(1);
m = coder.internal.indexInt(size(A,1));
n = coder.internal.indexInt(size(A,2));
[A,ipiv] = coder.internal.lapack.xgetrf(m,n,A,ONE,m);
P = double(eml_ipiv2perm(ipiv,m));
[L,U] = expandlu(A);

%--------------------------------------------------------------------------

function [L,U,P] = LUP(A)
ONE = coder.internal.indexInt(1);
m = coder.internal.indexInt(size(A,1));
n = coder.internal.indexInt(size(A,2));
[A,ipiv] = coder.internal.lapack.xgetrf(m,n,A,ONE,m);
[L,U,P] = expandlu(A,eml_ipiv2perm(ipiv,m));

%--------------------------------------------------------------------------

function [L,U,P] = expandlu(X,pivot)
% Form L, U, and possibly P matrices.  If pivot is supplied and nargout<3,
% the pivot vector is applied to L. Otherwise, P is the corresponding
% permutation matrix.
m = size(X,1);
n = size(X,2);
mn = eml_min(m,n);
L = eml_expand(coder.internal.scalarEg(X),[m,mn]);
U = eml_expand(coder.internal.scalarEg(X),[mn,n]);
% Extract the upper triangular matrix U.
for j = 1:mn
    for i = 1:j
        U(i,j) = X(i,j);
    end
end
for j = mn+1:n
    for i = 1:mn
        U(i,j) = X(i,j);
    end
end
if nargin == 2 && nargout < 3
    % Form the lower triangular matrix L and apply the pivot vector.
    for j = 1:mn
        L(pivot(j),j) = 1;
        for i = j+1:m
            L(pivot(i),j) = X(i,j);
        end
    end
else
    % Extract the lower triangular matrix L.
    for j = 1:mn
        L(j,j) = 1;
        for i = j+1:m
            L(i,j) = X(i,j);
        end
    end
    % Form the permutation matrix P if required.
    if nargout == 3
        P = zeros(m,'like',real(X));
        if nargin == 2
            for j = 1:m
                P(j,pivot(j)) = 1;
            end
        else
            for j = 1:m
                P(j,j) = 1;
            end
        end
    end
end

%--------------------------------------------------------------------------
