function [X, rCond] = hMldivideAndMrdivide(transposeFlagA, A, B, throwMrdivideErrors)
%hMldivideAndMrdivide is a common implementation of mldivide and mrdivide
%   for (co)distributed matrices. Optional second output is reciprocal
%   condition number (rCond).
%   


%   Copyright 2012-2022 The MathWorks, Inc.

import distributedutil.TransposeFlags

nargoutchk(1, 2);

rCond = [];

% Treat logicals and chars as doubles
if ischar(A) || islogical(A)
    A = double(A);
end
if ischar(B) || islogical(B)
    B = double(B);
end

sizeA = size(A);
if isa(A, 'codistributed')
    % First, check incompatible sizes, and incompatible types and sparsity.
    iErrorForIncompatibleDimension(A, B, transposeFlagA);
    if (issparse(A) && isUnderlyingType(B, 'single')) || ...
            (issparse(B) && isUnderlyingType(A, 'single'))
        if throwMrdivideErrors
            ME = MException(message('MATLAB:mrdivide:sparseSingleNotSupported'));
        else
            ME = MException(message('MATLAB:mldivide:sparseSingleNotSupported'));
        end
        throwAsCaller(ME);
    end
    
    % Try to make B full for use within SCALAPACK, and make sure to throw a
    % nice error message if it is not possible.
    if issparse(B) && ~issparse(A)
        try % Error must be thrown on all workers
            B = distributedutil.syncOnError(@full, B);
        catch
            ME = MException(message('parallel:distributed:SparseOutOfMemory'));
            throwAsCaller(ME);
        end
    end
    
    % Handle mixed double-single inputs.
    if isUnderlyingType(A, 'double') && isUnderlyingType(B, 'single')
        A = single(A);
    elseif isUnderlyingType(A, 'single') && isUnderlyingType(B, 'double')
        B = single(B);
    end
    
    % For rectangular matrices, we need to transpose or conjugate
    % explicitly, if needed.
    if size(A,1) ~= size(A,2)
        if transposeFlagA == TransposeFlags.TRANSPOSE
            A = transpose(A);
        elseif transposeFlagA == TransposeFlags.CONJUGATE_TRANSPOSE
            A = ctranspose(A);
        end
    end
    
    % DECOMPOSITION does not support Banded and Tridiagonal with transpose
    % or ctranspose because ScaLAPACK is broken. Work around this by
    % transposing explicitly in these cases.
    if transposeFlagA ~= TransposeFlags.NO_TRANSPOSE && nnz(A) > 0
        % In this case, we would have called bandwidth anyway.
        [kl, ku] = bandwidth(A);
        nnzA = nnz(A);
        n = size(A,1);
        nnzFullBand = n*(kl+ku+1)-kl*(kl+1)/2-ku*(ku+1)/2;
        banddensity = 0.5; % default value
        if nnzA/nnzFullBand > banddensity && (kl+ku+1) < n/spmdSize
            % The matrix is a banded one, we transpose explicitly.
            if transposeFlagA == TransposeFlags.TRANSPOSE
                A = transpose(A);
            elseif transposeFlagA == TransposeFlags.CONJUGATE_TRANSPOSE
                A = ctranspose(A);
            end
            % And we change the transpose flag.
            transposeFlagA = TransposeFlags.NO_TRANSPOSE;
        end
    end
    
    % Now call DECOMPOSITION and SOLVE.
    decA = decomposition(A);
    if size(A,1) ~= size(A,2) % rectangular
        X = decA\B;
    else % square
        X = underlyingSolveWithTransposeFlag(decA, B, transposeFlagA);
    end
    
    % If B has NaN and rcond is empty we should return NaN so that a
    % warning is thrown. Otherwise, it means DECOMPOSITION has already thrown a
    % warning if needed. Note DECOMPOSITION supports RCOND for square only.
    if size(A,1) == size(A,2)
        rCond = gather(rcond(decA));
        if hasNaN(B) && isempty(rCond)
            rCond = NaN;
        end
    end
else
    % A is not codistributed and B is codistributed.
    bDist = getCodistributor(B);
    if ~( isa(bDist, 'codistributor1d') && (bDist.Dimension == 2) )
        B = redistribute(B, codistributor1d(2));
        bDist = getCodistributor(B);
    end
    % Since A is not codistributed and B is distributed by columns, the
    % workers all have different linear systems to solve.  A new codistributor
    % will be needed if A is not square.  In any case, the results of the
    % divide operation will have the same column distribution as B, and
    % the number of rows will be the number of columns in the input A.
    % Warnings are issued by the core MATLAB function and rCond is not
    % modified.
    
    % The first dimension of the result depends on the size of A and the
    % used transpose flag of A.
    if transposeFlagA == TransposeFlags.NO_TRANSPOSE
        mRhs = sizeA(2);
    else
        mRhs = sizeA(1);
    end
    aDist = codistributor1d(2, bDist.Partition, [mRhs size(B, 2)]);
    localX = iSolveLocalSystem(transposeFlagA, A, getLocalPart(B), throwMrdivideErrors);
    X = codistributed.pDoBuildFromLocalPart(localX, aDist);
end
end

function iErrorForIncompatibleDimension(A, B, transposeFlagA)
if transposeFlagA == distributedutil.TransposeFlags.NO_TRANSPOSE
    dimToCompare = 1;
else
    dimToCompare = 2;
end

if size(A, dimToCompare) ~= size(B, 1)
    ME = MException(message('MATLAB:dimagree'));
    throwAsCaller(ME);
end
end

function X = iSolveLocalSystem(transposeFlagA, A, B, throwMrdivideErrors)
% Solve A*X = B, but throw any MLDIVIDE errors as MRDIVIDE errors if needed
import distributedutil.TransposeFlags
try
    switch transposeFlagA
        case TransposeFlags.NO_TRANSPOSE
            X = A\B;
        case TransposeFlags.TRANSPOSE
            X = A.'\B;
        case TransposeFlags.CONJUGATE_TRANSPOSE
            X = A'\B;
    end
catch ME
    if throwMrdivideErrors
        % This string replacement technique works as long as all the error
        % IDs continue to use the same theme (identical in all respects
        % except for rdivide/ldivide). Otherwise this code should change to
        % a switch statement.
        newMessageID = strrep( ME.identifier, 'ldivide', 'rdivide' );
        error( message( newMessageID ) );
    else
        rethrow(ME);
    end
end
end
