function [localX, flag, relres, iter, resvec, resveccg] = hMinresImpl(codistr, A, b, tol, maxit, M1, M2, x)
% HMINRESIMPL Implementation of minres function with minimum communication using
% algorithm described in Codistributed_sparse_iterative_solvers_RFA on the wiki.
% b is codistributed 1d by rows; A.fun is the local part of A distributed 1d
% by rows, using the same partition scheme as b; 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
% 5  Preconditioner M was not symmetric positive definite

%   Copyright 2017-2022 The MathWorks, Inc.

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

normB = distributedutil.vector2NormFromLocalParts(localB);
flag = 1;
iter = 0;
maxit = max(maxit, 1);
n = codistr.Cached.GlobalSize(1);

resvec = zeros(maxit+1,1);         % Preallocate vector for norms of residuals
if nargout >= 6
    resveccg = zeros(maxit+1, 1);
end
xmin = localX;                          % Iterate which has minimal residual so far
imin = 0;                          % Iteration at which xmin was computed
tolb = tol*normB;                    % Relative tolerance

% Compute with initial guess
solution = iterApplyA(A, localX, localX, 'notransp', partition);
% Step 2
r = localB - solution;
% Check if we have a solution
normRes = distributedutil.vector2NormFromLocalParts(r);
resvec(1) = normRes;                  % resvec(1) = norm(b-A*x0)
normResMin = normRes;
resveccg(1) = resvec(1);
relresmin = resvec(1)/normB;
if relresmin <= tol                 % Initial guess is a good enough solution
    flag = 0;
    relres = relresmin;
    resvec = resvec(1);
    resveccg = resveccg(1);
    return
end

% allocate local arrays
nlocal = size(r, 1);
v = repmat(r, 1, 3);
m = zeros(nlocal, 3);
Am = zeros(nlocal, 3);
index = 1:3;
[new, older] = deal(index(1), index(3));
% Step 3
[v(:, new), beta] = iPrecondition('mldivideWithRcondCheck', r, M1, M2, codistr);
if beta <= 0
    flag = 5;
    relres = relresmin;
    resvec = resvec(1);
    resveccg = resveccg(1);
    return
end
% Step 4
beta = sqrt(beta);
% Step 5
vv = v(:, new)/beta;
% Steps 6-8
[Am(:, older), alpha] = iterApplyA(A, vv, vv, 'notransp', partition);
% Step 9
v(:, new) = Am(:, older) - (alpha/beta)*r;
% Step 10
% Local reorthogonalization
numer = distributedutil.vectorProdFromLocalParts(vv, v(:, new));
denom = distributedutil.vectorProdFromLocalParts(vv, vv);
v(:, new) = v(:, new) + (numer/denom)*vv;
% Step 11
snprod = beta;
cs = -1;
sn = 0;
deltabar = beta;
epsilon = 0;
stag = 0;
moreSteps = 0;
maxMoreSteps = min([floor(n/50),5,n-maxit]);
% Step 12
for ii = 1:maxit
    % Step 12.1, 12.5, 12.22
    index = mod(index+1, 3)+1;
    [new, old, older] = deal(index(1), index(2), index(3));
    % Step 12.3
    betaold = beta;
    % Step 12.2:3
    [v(:, new), beta] = iPrecondition('mldivide', v(:, old), M1, M2, codistr);
    if beta < 0
        flag = 5;
        break
    end
    beta = sqrt(beta);
    % Step 12.4
    delta = cs*deltabar + sn*alpha;
    % Step 12.6
    m(:, new) = vv - delta*m(:, old) - epsilon*m(:, older);
    % Step 12.7
    Am(:, new) = Am(:, new) - delta*Am(:, old) - epsilon*Am(:, older);
    % Step 12.8
    gammabar = sn*deltabar - cs*alpha;
    % Step 12.9
    epsilon = sn*beta;
    % Step 12.10
    deltabar = -cs*beta;
    % Step 12.11
    gamma = sqrt(gammabar^2 + beta^2);
    % Step 12.12
    m(:, new) = m(:, new)/gamma;
    Am(:, new) = Am(:, new)/gamma;
    % Step 12.13
    cs = gammabar/gamma;
    sn = beta/gamma;
    sncs = snprod*cs;
    norms = distributedutil.vector2NormFromLocalParts([{m(:, new)} {localX}]);
    stag = iterCheckStag(stag, norms, sncs);
    % Step 12.14
    localX = localX + sncs*m(:, new);
    % Step 12.15
    if ~iterIsEmptyArg(M1)
        % Step 12.15.1
        r = r - sncs*Am(:, new);
        % Step 12.15.2
        snprod = snprod*sn;
        % Step 12.15.3
        normRes = distributedutil.vector2NormFromLocalParts(r);
    else
        % Step 12.16.1
        snprod = snprod*sn;
        % Step 12.16.2
        normRes = abs(snprod);
    end
    resvec(ii+1,1) = normRes;
    % Step 12.17
    if nargout >= 6
        if (cs == 0)
            % It's possible that this cs value is zero (CG iterate does not
            % exist)
            normRescg = Inf;
        else
            rcg = r - snprod*(sn/cs)*Am(:, new);
            normRescg = distributedutil.vector2NormFromLocalParts(rcg);
        end
        resveccg(ii+1, 1) = normRescg;
    end
    % Step 12.18
    [resvec, flag, stag, moreSteps, imin, xmin, normResMin, normRes] =...
        iterCheckConv(A, localB, localX, resvec, normRes, ii, tolb, stag, moreSteps, ...
        maxMoreSteps, imin, xmin, normResMin, flag, partition);
    if flag ~= 1 || ii == maxit
        iter = ii;
        break
    end
    % Step 12.19
    vv = v(:, new)/beta;
    % Step 12.20
    Am(:, older) = iterApplyA(A, vv, vv, 'notransp', partition);
    % Step 12.22
    v(:, new) = Am(:, older) - (beta/betaold)*v(:, older);
    % Step 12.23
    alpha = distributedutil.vectorProdFromLocalParts(vv, v(:, new));
    % Step 12.24
    v(:, new) = v(:, new) - (alpha/beta)*v(:, old);
end
% Returned solution is first with minimal residual
if(flag == 0)
    relres = normRes/normB;
else
    resComp = localB - iterApplyA(A, xmin, xmin, 'notransp', partition);
    normResComp = distributedutil.vector2NormFromLocalParts(resComp);
    if normResComp <= normRes
        localX = xmin;
        iter = imin;
        relres = normResComp/normB;
        resvec(imin+1) = normResComp;
    else
        iter = ii;
        relres = normRes/normB;
    end
end

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

function [localZ, rho] = iPrecondition(op, localR, M1, M2, codistr)
% IPRECONDITION applies the preconditioners (if supplied) to r.  It
% uses local parts for input and output.  It also calculates the residual.
hasAPreconditioner =  ~iterIsEmptyArg(M1) || ~iterIsEmptyArg(M2);
if hasAPreconditioner
    r = codistributed.build(localR, codistr, 'noCommunication');
    if ~iterIsEmptyArg(M1)
        z = distributedutil.iterapp(op, M1.fun, M1.type, ...
                           r, 'notransp');
    else
        z = r;
    end
    if ~iterIsEmptyArg(M2)
        z = distributedutil.iterapp(op, M2.fun, M2.type, ...
                           z, 'notransp');
    end
    localR = getLocalPart(r);
    localZ = getLocalPart(z);
else
    localZ = localR;
end
rho = distributedutil.vectorProdFromLocalParts(localR, localZ);
end
