function [localX, flag, relres, iter, resvec] = hBicgImpl(codistr, A, b, tol, maxit, M1, M2, x)
% Implementation of bicg function with minimum communication using
% algorithm described in Sparse_distributed_iterative_solvers on the wiki.
% b is codistributed 1d by rows; A.fun and A.tfun are the local parts of A and
% A' distributed 1d by rows, resp. columns, using the same partition scheme
% as b; and M1 and M2, if they are matrices, are distributed by rows. The
% "step" comments refer to lines in the algorithm on the wiki.

% The flag variable indicates the progress of the solution:
% 0  Converged
% 1  Not converged yet
% 2  There was an infinite value in the output from the preconditioner
% 3  Iteration has stagnated
% 4  A scalar parameter has become too small or too large

%   Copyright 2015-2022 The MathWorks, Inc.

% Message-passing info
[partitionStart, partitionEnd] = globalIndices(codistr, 1);
partition = matlab.internal.ColonDescriptor(partitionStart, partitionEnd);
n = codistr.Cached.GlobalSize(1);
localB = getLocalPart(b);
localX = getLocalPart(x);

normB = distributedutil.vector2NormFromLocalParts(localB);
tolb = tol*normB;
flag = 1;

resvec = zeros(maxit+1,1);         % Preallocate vector for norms of residuals
xmin = localX;                          % Iterate which has minimal residual so far
imin = 0;                          % Iteration at which xmin was computed

% Compute with initial guess
% Step 2
r = localB - iterApplyA(A, localX, localX, 'notransp', partition);
% Step 3
rt = r;                            % Shadow residual
% Check if we have a solution
normRes = distributedutil.vector2NormFromLocalParts(r);
resvec(1) = normRes;               % resvec(1) = norm(b-A*x0)
if normRes <= tolb
    resvec = resvec(1);
    flag = 0;
    relres = normRes/normB;
    iter = 0;
    return
end

% Steps 4 & 5
[p, pt, rho] = iPrecondition('mldivideWithRcondCheck', r, rt, M1, M2, codistr);
flag = iterCheckParam(rho, flag);
normResMin = normRes;
normResAct = normRes;
stag = 0;
moreSteps = 0;
maxMoreSteps = min([floor(n/50), 5, n-maxit]);

% Step 6
for ii = 1:maxit
    % Loop Step 1
    qt = iterApplyA(A, pt, pt, 'transp', partition, n);
    [q, ptq] = iterApplyA(A, p, pt, 'notransp', partition);
    flag = iterCheckParam(ptq, flag, 1);
    if flag>1
        break
    end
    alpha = rho/ptq;
    flag = iterCheckParam(alpha, flag, 2);
    if flag > 1
        break
    end
    norms = distributedutil.vector2NormFromLocalParts([{p} {localX}]);
    stag = iterCheckStag(stag, norms, alpha);
    % Loop Step 2
    localX = localX + alpha*p;
    % Loop Step 3
    r = r - alpha*q;
    % Loop Step 4
    rt = rt - conj(alpha)*qt;
    % Loop Step 5
    rho1 = rho;
    [z, zt, rho, normRes] = iPrecondition('mldivide', r, rt, M1, M2, codistr);
    resvec(ii+1) = normRes;
    [resvec, flag, stag, moreSteps, imin, xmin, normResMin, normResAct] =...
        iterCheckConv(A, localB, localX, resvec, normRes, ii, tolb, stag, moreSteps, ...
        maxMoreSteps, imin, xmin, normResMin, flag, partition);
    if flag ~= 1
        break
    end
    flag = iterCheckParam(rho, flag);
    if flag>1
        break
    end
    beta = rho/rho1;
    flag = iterCheckParam(beta, flag);
    if flag>1
        break
    end
    % Loop Step 6
    p = z + beta*p;
    % Loop Step 7
    pt = zt + conj(beta)*pt;
end
if isempty(ii)
    ii = 0;
end

% Returned solution is first with minimal residual
if flag == 0
    iter = imin;
    relres = normResAct/normB;
else
    r = localB - iterApplyA(A, xmin, xmin, 'notransp', partition);
    if flag == 1
        resComp = localB - iterApplyA(A, localX, localX, 'notransp', partition);
        norms = distributedutil.vector2NormFromLocalParts([{r} {resComp}]);
        normRes = norms(1);
        normResAct = norms(2);
    else
        normRes = distributedutil.vector2NormFromLocalParts(r);
    end
    if normRes <= normResAct
        localX = xmin;
        iter = imin;
        relres = normRes/normB;
    else
        iter = ii;
        relres = normResAct/normB;
    end
end 

% Truncate the zeros from resvec
if flag <= 1 || flag == 3
    resvec = resvec(1:ii+1);
else
    resvec = resvec(1:ii);
end
end

function [localZ, localZt, ptq, res] = iPrecondition(op, localR, localRt, M1, M2, codistr)
% IPRECONDITION applies the preconditioners (if supplied) to r and rt.  It
% uses local parts for input and output.  It also calculates the norm of
% the input residual.
if nargout > 3
    res = distributedutil.vector2NormFromLocalParts(localR);
end
hasAPreconditioner =  ~iterIsEmptyArg(M1) || ~iterIsEmptyArg(M2);
if hasAPreconditioner
    r = codistributed.build(localR, codistr, 'noCommunication');
    rt = codistributed.build(localRt, codistr, 'noCommunication');
    if ~iterIsEmptyArg(M1)
        z = distributedutil.iterapp(op, M1.fun, M1.type, ...
                           r, 'notransp');
    else
        z = r;
    end
    if ~iterIsEmptyArg(M2)
        zt = distributedutil.iterapp(op, M2.tfun, M2.type, ...
                           rt, 'transp');
        z = distributedutil.iterapp(op, M2.fun, M2.type, ...
                           z, 'notransp');
    else
        zt = rt;
    end
    if ~iterIsEmptyArg(M1)
        zt = distributedutil.iterapp(op, M1.tfun, M1.type, ...
                           zt, 'transp');
    end
    localZ = getLocalPart(z);
    localZt = getLocalPart(zt);
    localRt = getLocalPart(rt);
else
    localZ = localR;
    localZt = localRt;
end
ptq = distributedutil.vectorProdFromLocalParts(localRt, localZ);
end
