function [A, p] = scalaChol(A, enumUpLow)
% SCALACHOL  ScaLAPACK Cholesky factorization

%   Copyright 2006-2023 The MathWorks, Inc.

narginchk(2, 2);
nargoutchk(2, 2);

% Special shortcut for scalar case
if isscalar(A)
    if isreal(A) && A > 0
        A = sqrt(A);
        p = 0;
    else % complex or nonpositive
        A = zeros(0, 0, "like", real(A));
        p = 1;
    end
    return
end

% Redistribute input according to the default 2dbc distribution scheme.
dimToRedistribute = [];
if ~isa(getCodistributor(A), "codistributor2dbc")
    codA = getCodistributor(A);
    dimToRedistribute = codA.Dimension;
    A = redistribute(A, codistributor2dbc());
end

distA = getCodistributor(A);
lbgrid = distA.LabGrid;
descA = arraydescriptor(distA);
localA = hGetLocalPartWithCorrectComplexity(A);

% Note the internal builtin requires the flag to be passed as a char vector.
[localA, p] = parallel.internal.codistextern.denseCholFactor(localA, descA, ...
                                                  lbgrid(1), lbgrid(2), ...
                                                  distA.Orientation, ...
                                                  char(enumUpLow));
codistr = codistributor2dbc(lbgrid, distA.BlockSize,  ...
                            distA.Orientation, size(A));

A = codistributed.pDoBuildFromLocalPart(localA, codistr);

if ~isempty(dimToRedistribute)
    % We redistributed A before the computation and we need to redistribute
    % the output back to the original codistributor1d.
    A = redistribute(A, codistributor1d(dimToRedistribute));
end

% The leading N-by-N triangular part of A contains the factored matrix,
% (upper by default, but lower if there are two input arguments)
% and its strictly (lower, resp. upper) triangular part is not referenced.  
% Extract the factor of interest from A.
if enumUpLow == "upper"
    A = triu(A);
else % lower
    A = tril(A);
end
end
