classdef DenseCOD < matlab.mixin.internal.Scalar
    % DENSECOD   COD decomposition of a dense matrix
    %
    %   FOR INTERNAL USE ONLY -- This feature is intentionally undocumented.
    %   Its behavior may change, or it may be removed in a future release.
    %
    
    %   Copyright 2017-2020 The MathWorks, Inc.
    
    properties (Access = private)
        m_
        n_

        QR_
        tau1_
        perm_
        tau2_
    end
    
    properties (GetAccess = public, SetAccess = private)
        rank_ = [];
        ranktol_ = [];
    end
    
    methods
        function f = DenseCOD(A, tol)
            [f.m_,f.n_] = size(A);
            
            if f.m_ < f.n_
                A = A';
            end
            
            if isempty(tol)
                % Setting tol < 0 causes qrFactor to compute default tolerance from A.
                % Default tolerance is defined as:
                % If A is real:
                % tol = min(max(size(A))*eps(class(A)), sqrt(eps(class(A)))) * abs(R(1, 1))
                %
                % If A is complex:
                % tol = min(10*max(size(A))*eps(class(A)), sqrt(eps(class(A)))) * abs(R(1, 1))

                tol = -2;
            end
            
            % Construct object.
            [f.QR_,f.tau1_,f.perm_,f.rank_,f.ranktol_] = ...
                matlab.internal.decomposition.builtin.qrFactor(A, tol);
            
            if f.rank_ < min(f.m_, f.n_)
                [f.QR_, f.tau2_] = ...
                    matlab.internal.decomposition.builtin ...
                    .codFactor(f.QR_, f.rank_);
            end
            
        end
        
        function x = solve(f,b,transposed)
            
            if f.m_ < f.n_
                % The internally computed decomposition always has m >= n
                transposed = ~transposed;
            end
            
            if ~transposed
                % Solve Qlong * R * Qshort' * x = b, using formula
                % x = (Qshort * (R \ (Qlong' *b) ) )
                x = matlab.internal.decomposition.builtin.codSolve(f.QR_, f.tau1_, f.tau2_, f.rank_, b, transposed);
                
                % Apply permutation
                x(f.perm_, :) = x;
            else
                % Apply permutation
                b = b(f.perm_, :);
                
                % Solve Qshort * R' * Qlong' * x = b, using formula
                % x = (Qlong * (R' \ (Qshort' *b) ) )
                x = matlab.internal.decomposition.builtin.codSolve(f.QR_, f.tau1_, f.tau2_, f.rank_, b, transposed);
            end

        end
    end
end
