function varargout = hIluRowWise(codistrA, localA, setup)
%hIluRowWise incomplete LU of a sparse matrix distributed row wise.
%
%   LU = hIluRowWise(CODISTRA, LOCALA, SETUP) returns "L+U-speye(size(A))"
%
%   [L,U] = hIluRowWise(CODISTRA, LOCALA, SETUP) returns unit lower
%   triangular L and upper triangular U

%   Copyright 2020-2022 The MathWorks, Inc.

partitionA = codistrA.Partition;
blockLimits = cumsum(partitionA);
blockLimits = [0 blockLimits(1:end-1); blockLimits]';
firstDiag = blockLimits(spmdIndex, 1);
if isfield(setup, 'milu')
    assert(strcmpi(setup.milu, "off") || strcmpi(setup.milu, "row"), ...
        "setup.milu must be 'off' or 'row' for hIluRowImpl");
    if strcmp(setup.milu, 'row')
        % This is equivalent to 'col' for hIluColWise
        setup.milu = 'col';
    end
end

% A is distributed by rows, we call hIluColImplCpp on A.' and transpose
% back the result. As localA is stored by columns (CSC format), this is
% faster than implementing the row-wise counterpart of hIluColWise.
codistrA = codistrA.pTransposeCodistributor();
[localU, localL] = codistrA.hIluColWise(localA.', setup); % Swap L and U
localL = localL.';
localU = localU.';
% L is supposed to be unit-diagonal and U owns the diagonal entries.
% As we swapped L and U we have to rescale and unscale accordingly.
if isrow(localL)
    d = localL(1, firstDiag + 1);
else
    d = diag(localL, firstDiag);
end
localU = localU .* d;
d = spmdCat(d.', 2);
localL = localL ./ d;

if nargout == 1
    varargout = {tril(localL, -firstDiag-1) + localU};
else
    varargout = {localL, localU};
end
end
