function [Q,R,E] = eml_qr(A,varargin)
%MATLAB Code Generation Library Function

%   QR decomposition.  Optional arguments may be 0, for economy mode,
%   'matrix' to generate a permutation matrix for the third output, or
%   'vector' to generate a permutation vector for the third output.

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

coder.internal.assert(nargin > 0,'MATLAB:minrhs');
coder.internal.assert(isa(A,'float'), ...
    'Coder:toolbox:unsupportedClass','qr',class(A));
coder.internal.assert(ismatrix(A),'Coder:toolbox:qr_inputMustBe2D');
coder.internal.prefer_const(varargin);
[ECONOMY_MODE,PERMUTATION_MATRIX] = parseopts(varargin{:});
if nargout == 3
    if ECONOMY_MODE
        [Q,R,jpvt] = qr_econ(A);
    else
        [Q,R,jpvt] = qr_full(A);
    end
    if PERMUTATION_MATRIX
        % Create a permutation matrix from jpvt.
        E = zeros(size(A,2),class(A));
        for k = 1:size(A,2)
            E(jpvt(k),k) = 1;
        end
    else
        % Return integer jpvt vector.
        E = jpvt;
    end
elseif nargout == 2
    if ECONOMY_MODE
        [Q,R] = qr_econ(A);
    else
        [Q,R] = qr_full(A);
    end
else
    Q = coder.internal.lapack.xgeqrf(A);
end

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

function [ECONOMY_MODE,PERMUTATION_MATRIX] = parseopts(varargin)
% Check the optional inputs to determine whether this is economy mode and
% whether to generate a permutation matrix for the third output.
coder.internal.prefer_const(varargin);
coder.internal.assert(nargin <= 2,'Coder:toolbox:eml_qr_3');
ECONOMY_MODE = false;
% The default for PERMUTATION_MATRIX is false in economy mode, otherwise
% true.  If any inputs are supplied, either this is economy mode or one
% the optional inputs will override the following initialization.
PERMUTATION_MATRIX = (nargin == 0);
for k = coder.unroll(1:nargin)
    coder.internal.assert(eml_is_const(varargin{k}),'Coder:toolbox:eml_qr_4');
    if coder.internal.isTextRow(varargin{k})
        len = max(1,strlength(varargin{k}));
        VECTOR = strncmpi(varargin{k},'vector',len);
        PERMUTATION_MATRIX = strncmpi(varargin{k},'matrix',len);
        coder.internal.assert(VECTOR || PERMUTATION_MATRIX, ...
            'MATLAB:qr:unknownOptionForSecondInput');
    elseif coder.internal.isBuiltInNumeric(varargin{k}) && ...
            isscalar(varargin{k}) && varargin{k} == 0
        ECONOMY_MODE = true;
    else
        coder.internal.assert(false,'MATLAB:qr:unknownOptionForSecondInput');
    end
end

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

function [Q,R,jpvt] = qr_full(A)
% Full QR with or without pivoting. To return a permutation vector for E,
% pass in PERMUTATION_VECTOR as the second input.  The default OPT is
% PERMUTATION_MATRIX.
ONE = ones(coder.internal.indexIntClass);
coder.inline('always');
m = cast(size(A,1),coder.internal.indexIntClass);
n = cast(size(A,2),coder.internal.indexIntClass);
if ~eml_is_const(n)
    coder.varsize('Q','tau');
end
% Q is m x m, R is m x n
Q = coder.nullcopy(eml_expand(coder.internal.scalarEg(A),[m,m]));
R = coder.nullcopy(A);
if m > n
    % Q = [A,zeros(m,m-n)];
    for j = 1:n
        for i = 1:m
            Q(i,j) = A(i,j);
        end
    end
    for j = n+1:m
        for i = 1:m
            Q(i,j) = 0;
        end
    end
    if nargout == 3
        [Q,tau,jpvt1] = coder.internal.lapack.xgeqp3(Q);
        jpvt = coder.nullcopy(zeros(1,n,class(jpvt1)));
        for j = 1:n
            jpvt(j) = jpvt1(j);
        end
    else
        [Q,tau] = coder.internal.lapack.xgeqrf(Q);
    end
    for j = 1:n
        for i = 1:j
            R(i,j) = Q(i,j);
        end
        for i = coder.internal.indexPlus(j,1):m
            R(i,j) = 0;
        end
    end
    Q = unpackQR(m,m,n,Q,ONE,m,tau,ONE);
else
    if nargout == 3
        [A,tau,jpvt] = coder.internal.lapack.xgeqp3(A);
    else
        [A,tau] = coder.internal.lapack.xgeqrf(A);
    end
    for j = 1:m
        for i = 1:j
            R(i,j) = A(i,j);
        end
        for i = coder.internal.indexPlus(j,1):m
            R(i,j) = 0;
        end
    end
    for j = coder.internal.indexPlus(m,1):n
        for i = 1:m
            R(i,j) = A(i,j);
        end
    end
    A = unpackQR(m,m,m,A,ONE,m,tau,ONE);
    for j = 1:m
        for i = 1:m
            Q(i,j) = A(i,j);
        end
    end
end

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

function [Q,R,jpvt] = qr_econ(A)
% Economy size QR with pivoting.
coder.inline('always');
ONE = ones(coder.internal.indexIntClass);
if nargout == 3
    [A,tau,jpvt] = coder.internal.lapack.xgeqp3(A);
else
    [A,tau] = coder.internal.lapack.xgeqrf(A);
end
m = cast(size(A,1),coder.internal.indexIntClass);
n = cast(size(A,2),coder.internal.indexIntClass);
minsz = eml_min(m,n);
R = coder.nullcopy(eml_expand(coder.internal.scalarEg(A),[minsz,n]));
% if m >= n, Q is m x n, R is n x n
% if m <= n, Q is m x m, R is m x n
for j = 1:minsz
    for i = 1:j
        R(i,j) = A(i,j);
    end
    for i = (j + 1):minsz
        R(i,j) = 0;
    end
end
for j = (m + 1):n
    for i = 1:minsz
        R(i,j) = A(i,j);
    end
end
A = unpackQR(m,minsz,minsz,A,ONE,m,tau,ONE);
Q = coder.nullcopy(eml_expand(coder.internal.scalarEg(A),[m,minsz]));
for j = 1:minsz
    for i = 1:m
        Q(i,j) = A(i,j);
    end
end

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

function A = unpackQR(m,n,k,A,ia0,lda,tau,itau0)
% Dispatch to the proper function to unpack the QR decomposition based on the
% complexity of A
coder.inline('always');
if isreal(A)
    A = coder.internal.lapack.xorgqr(m,n,k,A,ia0,lda,tau,itau0);
else
    A = coder.internal.lapack.xungqr(m,n,k,A,ia0,lda,tau,itau0);
end

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