function L = ichol(A, varargin)
%ICHOL  Sparse Incomplete Cholesky factorization of codistributed array
%   
%   L = ichol(A)
%   L = ichol(A,SETUP)
%   
%   Limitations:
%   The "type" field of the structure "SETUP", if provided, must be "nofill".
%   The input matrix must be real symmetric or complex hermitian.
%   
%   Example:
%   spmd
%       A = codistributed(gallery('poisson', 1600));
%       setup.type = 'nofill';
%       L = ichol(A, setup);
%   
%   See also ICHOL, CODISTRIBUTED, CODISTRIBUTED/ILU, GALLERY.
%   


%   Copyright 2020-2022 The MathWorks, Inc.

narginchk(1, 2);
nargoutchk(0, 1);

% Make sure all arguments other than A are gathered.
if nargin > 1
    distributedutil.CodistParser.verifyNonCodistributedInputs([{A}, varargin]);
    varargin = distributedutil.CodistParser.gatherElements(varargin);
end

% If A isn't codistributed, call MATLAB function.
if ~iscodistributed(A)
    % Only the nondata argument was codistributed.
    L = ichol(A, varargin{:});
    return
end

% Check syntax
% setup = iCheckArgs(A, varargin);
try
    setup = iCheckArgs(A, varargin);
catch E
   throw(E);
end

codistrA = getCodistributor(A);
codistributed.pVerifyUsing1d("ilu", A);
distDim = codistrA.Dimension;
localA = getLocalPart(A);

% If we already know that A is symmetric (or hermitian), we let the
% computational routines know to avoid one extra check.
isMatrixKnownSymmetric = ...
    A.Metadata.Symmetry.isKnown() && ...
    (issymmetric(A) || ishermitian(A));

if distDim == 2
    % A is distributed by columns
    localL = codistrA.hIcholColWise(localA, setup, isMatrixKnownSymmetric);
else
    % A is distributed by rows
    localL = codistrA.hIcholRowWise(localA, setup, isMatrixKnownSymmetric);
end
L = codistributed.pDoBuildFromLocalPart(localL, codistrA);
end

function setup = iCheckArgs(A, argList)
% Check that A is a sparse square matrix and that setup, if given, has
% field values we support
if ~issparse(A)
    error(message("MATLAB:ichol:FirstArgNotSparse"));
end
if ~isUnderlyingType(A, "double")
    error(message("MATLAB:ichol:FirstArgNotDouble"));
end
[m, n] = size(A);
if m ~= n
    error(message("MATLAB:square"));
end
% Ensure consistency for the 1 worker edge case, triangular inputs lead to
% mismatch when creating the communication pattern.
if spmdSize == 1
    if (istriu(A) || istril(A)) && ~ishermitian(A)
        error(message("parallel:distributed:IcholNotSymmetric"));
    end
end

% Use ICHOL itself on some in-memory data to get free argument parsing!
% Note that we are already try/catch-ing in the upper layer here.
ichol(speye(1), argList{:});

if ~isempty(argList)
    setup = argList{1};
    if isfield(setup, "type") && ~strcmpi(setup.type, "nofill")
        error(message("parallel:distributed:IluIcholUnsupportedSetupTypeOption", "ICHOL"));
    end
else
    setup = struct('type', 'nofill', ...
        'michol', 'off', ...
        'shape', 'lower');
end
% Partially filled setup.
if ~isfield(setup, "type")
    setup.type = "nofill";
end
if ~isfield(setup, "michol")
    setup.michol = "off";
end
if ~isfield(setup, "shape")
    setup.shape = "lower";
end

end
