function varargout = hIluColImplLegacy(codistrA, localA, setup)
%HILUCOLIMPL Calculate the incomplete LU decomposition of a sparse matrix
% distributed by columns and return in single-matrix form.  This algorithm
% uses the KIJ variant of gaussian elimination.  

% Blocks of rows are processed according to the partition scheme.
% One worker will have a square block with the main diagonal running down
% the middle of it.  Other workers will be entirely in L or U.
%   Inputs: A     - distributed sparse square matrix to be decomposed
%           setup - structure containing control parameters
%   Single Output - decomposed matrix with L and U packed into one array
%   Two outputs   - local parts of L and U

%   Copyright 2019-2022 The MathWorks, Inc.

distDim = codistrA.Dimension;
partitionA = codistrA.Partition;
blockLimits = cumsum(partitionA);
blockLimits = [0 blockLimits(1:end-1); blockLimits]';
firstDiag = blockLimits(spmdIndex, 1);
mwTag = 2^15 - mod(sum(mfilename), 2^15) + (1:spmdSize);
[numRows, numCols] = size(localA);
fillin = [];
dim2sum = 0;
sparseMap = localA ~= 0;

% Update sparseMap to allow fillin on diagonal of A
diagLimits = blockLimits(spmdIndex, 1) + 1:blockLimits(spmdIndex, 2);
sparseMap(diagLimits, :) = sparseMap(diagLimits, :) | speye(numCols);

% Set dimension along which to sum fillin values if we are doing modified
% ilu.  Leave as 0 if we're not doing modified ilu.
if isfield(setup, 'milu')
    switch(setup.milu)
        case 'col'
            dim2sum = 1;
            fillin = zeros(1, numCols);
        case 'row'
            dim2sum = 2;
            fillin = zeros(numRows, 1);
    end
end
zeroPivot = false;
            
% Three phases: entirely in U, on diagonal, entirely in L
for kBlock = 1:spmdSize
    kLimits = blockLimits(kBlock, :);
    ks = kLimits(1)+1:kLimits(2);
    if isempty(ks)
        localL = 0;
    end
    for k = ks
        if spmdIndex > kBlock % in U
            if dim2sum == distDim
                spmdPlus(fillin(k), kBlock);  % result needed only on diagonal
            end
            if k < numRows
                [localA, fillin] = iUpdateU(k, kBlock, mwTag, localA, sparseMap, ...
                    dim2sum, fillin);
            end
        elseif spmdIndex == kBlock % block on diagonal
            localCol = k - firstDiag;
            if dim2sum == distDim
                fillin(k) = spmdPlus(fillin(k), kBlock); %#ok<AGROW>
                localA(k, localCol) = localA(k, localCol) - fillin(k);
            elseif dim2sum == 3 - distDim
                localA(k, localCol) = localA(k, localCol) - fillin(localCol);
            end
            if k == numRows
                continue
            end
            [localA, fillin, pivot] = iUpdateDiagBlock(k, localCol, mwTag, localA, ...
                 sparseMap, dim2sum, fillin);
            zeroPivot = zeroPivot || pivot == 0;
        elseif dim2sum == distDim  % block is in L and dim2sum == distDim otherwise, nothing to do
            spmdPlus(fillin(k), kBlock); % result needed only on diagonal
        end
    end
    % All workers check for zero pivot so that we can throw a coordinated
    % error.
    if any(spmdCat(zeroPivot))
        % Throw error.
        ex = MException(message("parallel:distributed:IluZeroPivot"));
        throwAsCaller(ex);
    end
    if spmdIndex == kBlock 
        % Store the part of localA that is in L in localL and save index
        % for assembling U for the output.
        localL = tril(localA(ks, :), -1);
        localU = triu(localA(ks, :), 0);
    end
end

% Prepare output
if nargout == 1
    varargout{1} = localA;
else
    emptyPartOfL = sparse(blockLimits(spmdIndex, 1), numCols);
    L = [emptyPartOfL; localL + speye(numCols);...
        localA((blockLimits(spmdIndex, 2)+1):end, :)];
    varargout{1} = L;
    emptyPartOfU = sparse(numRows - blockLimits(spmdIndex, 2), numCols);
    U = [localA(1:blockLimits(spmdIndex, 1), :); localU;...
        emptyPartOfU];
    varargout{2} = U;
end
end

function fillin = iUpdateFillin(fillin, filled, dim2sum, numRows, numCols, update)
% This function updates the fillin vector for cases where it is being
% maintained

[ifilled, jfilled] = find(filled);
fillinUpdate = sum(sparse(ifilled, jfilled, update(filled), numRows,...
    numCols), dim2sum);
fillin = fillin + fillinUpdate;
end

function [localA, fillin] = iUpdateU(k, kBlock, mwTag, localA, sparseMap, ...
    dim2sum, fillin)
% Function to update a column in the upper triangular part of A

[numRows, numCols] = size(localA);
bufferA = spmdReceive(kBlock, mwTag(spmdIndex) + k); % spmdSend is in diagonal block.
update = [sparse(k, 1); bufferA]*localA(k, :);
sparseMapUpdate = update~=0;
filled = sparseMapUpdate > sparseMap;
if any(filled, 'all')
    if dim2sum
        fillin = iUpdateFillin(fillin, filled, dim2sum, numRows,...
            numCols, update); 
    end
    update(filled) = 0;
end
localA = localA - update;
end

function [localA, fillin, pivot] = iUpdateDiagBlock(k, localCol, mwTag, localA, ...
    sparseMap, dim2sum, fillin)
% Function to update a block that spans the main diagonal of A

[numRows, numCols] = size(localA);
pivot = full(localA(k, localCol));
localA(k+1:end, localCol) = localA(k+1:end, localCol)/...
    localA(k, localCol);
if isnan(pivot)  % May be fill in localCol
    filled = (localA(:, localCol) ~= 0) > sparseMap(:, localCol);
    localA(filled, localCol) = 0;
end
for recvr = spmdIndex+1:spmdSize
    spmdSend(localA(k+1:end, localCol), recvr, mwTag(recvr) + k);
    % Corresponding spmdReceives are in the workers processng U blocks
    % to the right.
end
update = [sparse(k, 1); localA(k+1:end, localCol)]*...
    [sparse(1, localCol) localA(k, localCol+1:end)];
sparseMapUpdate = update~=0;
filled = sparseMapUpdate > sparseMap;
if any(filled, 'all')
    if dim2sum
        fillin = iUpdateFillin(fillin, filled, dim2sum, numRows,...
            numCols, update); 
    end
    update(filled) = 0;
end
localA = localA - update;
end
