function x = lsqminnorm(A, b, varargin)
%LSQMINNORM Minimum-norm solution of least-square system
%    X = LSQMINNORM(A, B) returns a vector X that minimizes norm(A*X - B).
%    If there are many solutions X to this problem, then the solution with
%    minimal norm(X) is returned. If B has multiple columns, the previous
%    statements are true for each column of X and B, respectively.
%
%    X = LSQMINNORM(A, B, TOL) additionally specifies the tolerance, which
%    is used to determine the rank of A. By default, TOL is computed based
%    on the QR decomposition of A.
%
%    LSQMINNORM(A, B, TOL) is a possible alternative to PINV(A, TOL) * B.
%    It is supported for sparse matrices, and is typically more efficient.
%    The two functions are not exactly equivalent, since LSQMINNORM uses
%    the Complete Orthogonal Decomposition (COD) to find a low-rank
%    approximation of A, while PINV uses the Singular Value Decomposition
%    (SVD).
%
%    X = LSQMINNORM(..., RANKWARN) specifies whether LSQMINNORM should
%    produce a warning when the matrix A is detected to be of low rank.
%    RANKWARN can be:
%        'nowarn' - (default) No warning is given if A has low rank.
%          'warn' - A warning is given if A has low rank.
%
%    See also: PINV, DECOMPOSITION, PAGELSQMINNORM

%   Copyright 2017-2023 The MathWorks, Inc.

narginchk(2, 4);

tol = -2; % this signals to internal classes to use default tolerance
checkCondition = false;

if ~isfloat(A)
    error(message('MATLAB:lsqminnorm:InvalidA'));
end

if ~isfloat(b)
    error(message('MATLAB:decomposition:InvalidB'));
end

if ~ismatrix(A) || ~ismatrix(b)
    error(message('MATLAB:lsqminnorm:inputMustBe2D'))
end

if nargin > 2
    offset = 1;
    if isnumeric(varargin{offset})
        tol = varargin{1};
        if ~isscalar(tol) || ~isreal(tol) || ~isfloat(tol) || ~(tol >= 0)
            error(message('MATLAB:lsqminnorm:InvalidTol'));
        end
        offset = offset+1;
    end

    if offset <= length(varargin)
        name = varargin{offset};
        if matlab.internal.math.partialMatch(name, 'warn')
            checkCondition = true;
        elseif matlab.internal.math.partialMatch(name, 'nowarn')
            checkCondition = false;
        else
            error(message('MATLAB:lsqminnorm:InvalidWarn'));
        end
    end

    if offset + 1 <= length(varargin)
        error(message('MATLAB:lsqminnorm:InvalidOption'));
    end
end

castBack = @(x) x;
if isa(A, 'single') ~= isa(b, 'single')
    if issparse(A) || issparse(b)
        error(message('MATLAB:mldivide:sparseSingleNotSupported'));
    end
    if isa(b, 'single')
        % Note: A\b would cast matrix A to single - but here the decomposition
        % is already done in double, so return x = single( dA\double(b) ) instead.
        b = double(b);
        castBack = @single;
    else
        b = single(b);
    end
end

if ~issparse(A) && issparse(b)
    b = full(b);
end

if size(b, 1) ~= size(A, 1)
    error(message('MATLAB:decomposition:mldivide'));
end

[m, n] = size(A);
if ~issparse(A)

    if m < n
        A = A';
    end

    % QR decomposition of A: A = Q*R*P
    [QR,Q1extra,perm,k,ranktol] = matlab.internal.decomposition.builtin.qrFactor(A, tol);

    if k < min(m, n)
        % Second QR factorization of R(1:k, :)'
        % A = Qlong * R(1:k, 1:k) * Qshort' * P
        [QR, tau2] = matlab.internal.decomposition.builtin.codFactor(QR, k);
    else
        tau2 = [];
    end
    % If we transposed A, then we
    %   Solve Qlong * R * Qshort' * x = b, using formula
    %   x = (Qshort * (R \ (Qlong' *b) ) )
    %   Next, we apply the permutation to x
    %   x(perm,:) = x
    % Otherwise, we do the following
    %   Apply our permutation to b
    %   b = b(perm,:)
    %   Then, solve Qshort * R' * Qlong' * x = b, using formula
    %   x = (Qlong * (R' \ (Qshort' *b) )
    x = matlab.internal.decomposition.builtin.codSolve(QR, Q1extra, tau2, k, b, m < n, perm);

    if checkCondition && k < min(size(A))
        warning(message('MATLAB:rankDeficientMatrix',sprintf('%d',k),sprintf('%13.6e',ranktol)));
    end

else % issparse(A)

    if m < n
        dA = matlab.internal.decomposition.SparseCOD(A, double(tol), true);

        x = solve(dA, b, false);

        if checkCondition && dA.rank_ < min(size(A))
            warning(message('MATLAB:rankDeficientMatrix',sprintf('%d',dA.rank_),sprintf('%13.6e',dA.ranktol_)));
        end
    else

        % More efficient variant without using internal decomposition objects
        % This avoids having to store the Householder vectors even implicitly
        [R, colperm1, QTb, rank_, ranktol_] = ...
            matlab.internal.math.sparseQRnoQ(A, true, double(tol), b);

        if rank_ == min(m, n)
            x = QTb(1:rank_, :);

            % x2 = R(1:rank_, 1:rank_) \ x;
            x2 = matlab.internal.decomposition.builtin.sparseTriangSolve(R(1:rank_,1:rank_), x, 'upper', false);

            x = zeros(n, size(b, 2), 'like', b);

            x(colperm1, :) = x2;

        else

            if checkCondition
                warning(message('MATLAB:rankDeficientMatrix',sprintf('%d',rank_),sprintf('%13.6e',ranktol_)));
            end

            M = R;
            M(rank_+1:end, :) = [];

            [H2, tau2, rowperm2, R, colperm2] = ...
                matlab.internal.math.implicitSparseQR(M', true, 0);

            rowperm2(colperm1) = rowperm2;

            x = QTb(colperm2, :);

            %x2 = R' \ x;
            x2 = matlab.internal.decomposition.builtin.sparseTriangSolve(R, x, 'upper', true);

            x = zeros(n, size(b, 2), 'like', b);

            x(1:rank_, :) = x2;

            if rank_ < min(m, n)
                if isreal(H2) && isreal(tau2) && ~isreal(x)
                    % Real Q applied to complex x not supported in built-in
                    x = matlab.internal.math.applyHouseholder(H2, tau2, rowperm2, real(x), false) + ...
                        1i*matlab.internal.math.applyHouseholder(H2, tau2, rowperm2, imag(x), false);
                else
                    x = matlab.internal.math.applyHouseholder(H2, tau2, rowperm2, x, false);
                end
            end
        end
    end
end

x = castBack(x);
