function [s,posdef,cgiters,znrm,exitflag,lambda,msg,PT,LZ,UZ,pcolZ,PZ] = ...
    ppcgr(grad,H,mtxmpy,Aeq,b,options,defaultopt,verb,computeLambda,tolA,PT, ...
    LZ,UZ,pcolZ,PZ,varargin)
%PPCGR	Preconditioned conjugate gradients with linear equalities.
%
% [s,posdef,cgiters] = PPCGR(grad,H,A,b,options)  apply
% a preconditioned conjugate gradient procedure to the linearly
% constrained quadratic
%
%         q(s) = .5s'Hs + grad's: As = b.
%
% On output s is the computed direction, posdef = 1 implies
% only positive curvature (in H constrained to null(A))
% has been detected; posdef = 0
% implies s is a direction of negative curvature (for M in null(A)).
% Output parameter cgiters is the number of CG-iterations used (which
% corresponds to the number of multiplications with H).
% A preconditioner will be based on the form
%
%            HH  AA'
%    M =
%            AA  0
%
% where in general HH is a SPD banded approximation
% to H, the nonzeros of AA are a subset of the nonzeros of A (based
% on a stopping tolerance). Parameters in the named parameter
% list options can be used to override default values to define M,
% including defining M (H) implicitly.
%
% [s,posdef,cgiters] = PPCGR(grad,H,A,b,options,tolA)
% Override the dropping tolerance to define AA from A.
%
% [s,posdef,cgiters] = PPCGR(grad,H,A,b,options,tolA,PT,LZ,UZ,pcolZ,PZ);
% PT is a permutation matrix such that A*PT has a leading nonsingular
% m-by-m matrix, A_1, where m is the number of rows of A
% PZ*A_1(:,pcolZ) = LZ*UZ, the sparse LU-factorization on A_1.
%
% [s,posdef,cgiters,znrm,exitflag,PT,LZ,UZ,pcolZ,PZ]=PPCGR(grad,H,A,b,options)
% outputs the sparse LU-factorization of A_1. znnrm is the scaled norm
% of the residual, exitflag is the termination code.

%   Copyright 1990-2020 The MathWorks, Inc.

% Import the warning free mldivide to avoid disabling / enabling the nearly
% singular warnings
import matlab.internal.math.nowarn.mldivide

if nargin < 3
    error(message('optimlib:ppcgr:NotEnoughInputs'));
end
[m,n] = size(Aeq);

if m > n
    exitflag = -10; s = []; lambda = []; posdef = []; cgiters = []; znrm = [];
    msg = [];
    % m > n not handled by ppcgr; abort
    return
end

exitflag = [];

if nargin < 4 || isempty(b)
    b = zeros(m,1);
end

if size(b,1) ~= m
    error(message('optimlib:ppcgr:SizeMismatch'));
end

if nargin < 7
    computeLambda = false;
    if nargin < 6
        verb = 0;
        if nargin < 5
            options = [];
        end
    end
end

% In case the defaults were gathered from calling: optimset('quadprog'):
numberOfVariables = n;
% Pre-allocate the lambda structure for the constraint types not handled
% here (or set to empty if lambda is not required).
if computeLambda
    lambda.ineqlin = zeros(0,1); % An empty column-vector
    lambda.lower = zeros(numberOfVariables,1);
    lambda.upper = zeros(numberOfVariables,1);
else
    lambda = [];
end

% Add Preconditioner to defaultopt. It's not a documented option, and is
% not added to defaultopt at the user-facing function level.
defaultopt.Preconditioner = @preaug;

pcmtx = optimget(options,'Preconditioner',defaultopt,'fast');
pcf = optimget(options,'PrecondBandWidth',defaultopt,'fast');
if nargin >= 9 && ~isempty(tolA)
    pcf(2,1) = tolA; % tolA not currently used
end
tol = optimget(options,'TolPCG', defaultopt,'fast') ;
defaultopt.MaxPCGIter = n; % not the default so override
kmax = optimget(options,'MaxPCGIter', defaultopt,'fast') ;
if ischar(kmax) && strcmpi(kmax,'max(1,floor(numberofvariables/2))')
    kmax = max(1,floor(numberOfVariables/2));
end
kmax = max(kmax,1);   % kmax must be at least 1

% save initial grad for computing Lagrange multipliers
gradsave = grad;

% Get feasible point
[y, flag] = feasibl(Aeq,b);
% We don't yet handle this case of dependent rows; abort
if flag == -1
    exitflag = -10; s = []; lambda = []; posdef = []; cgiters = []; znrm = [];
    msg = [];
    return
end
w = feval(mtxmpy,H,y,varargin{:});
grad = grad + w;
% If Aeq is a (square) non-singular matrix
if (m == n)
    posdef = []; cgiters = 0; znrm = 0;
    exitflag = 1;
    msg = qpTrustRegionMsg('optimlib:ppcgr:ExitSolvedSquare','one_feasible_point',verb > 0,'');
    s = y;
    if computeLambda
        rhs = -gradsave-w;
        lambda.eqlin = Aeq'\rhs;
    end
    return
end


% Transform to get fundamental null basis
if nargin < 10 || isempty(PT)
    PT = findp(Aeq);
end

A = Aeq*PT';
grad = PT*grad;

% For clarity, setup logical flags for all calls to fzmult.
doTranspose = true;
noTranspose = false;

% Project grad
if nargin < 14 || isempty(LZ) % last 4 arguments are a group (excluding varargin)
    [g,LZ,UZ,pcolZ,PZ] = fzmult(A,grad,doTranspose);
else
    g = fzmult(A,grad,doTranspose,LZ,UZ,pcolZ,PZ);
end
r = -g;

nz = length(g);
s = zeros(nz,1);

% Compute preconditioner
HH = [];
if isequal(size(H),[n,n])
    HH = PT*H*PT';
end
[L,U,P,pcol] = feval(pcmtx,HH,pcf,A,varargin{:});

% Precondition
rr = [zeros(m,1);r];
rhs = [rr;zeros(m,1)];
zz = L\(P*rhs);
zt(pcol,1) = U\zz;
z = zt(m+1:n,1);
znrm = norm(z);
stoptol = tol*znrm;
inner1 = r'*z;
posdef = 1;

% Flag for printing exit message
printMsg = verb > 0;

d = z;
% PRIMARY LOOP.
for k = 1:kmax
    dd = fzmult(A,d,noTranspose,LZ,UZ,pcolZ,PZ);
    dd = PT'*dd;
    w = feval(mtxmpy,H,dd,varargin{:});
    ww = PT*w;
    w = fzmult(A,ww,doTranspose,LZ,UZ,pcolZ,PZ);
    denom = d'*w;
    if denom <= 0
        if norm(d) > 0
            s = d/norm(d);
        else
            s = d;
        end
        s = y + PT'*fzmult(A,s,noTranspose,LZ,UZ,pcolZ,PZ);
        posdef = 0;
        if computeLambda
            w = feval(mtxmpy,H,s,varargin{:});
            rhs = -gradsave-w;
            lambda.eqlin = Aeq'\rhs;
        end
        if denom == 0 && sum(d)== 0 % No curvature left: we have a (singular) solution
            exitflag = 4;
            msg = qpTrustRegionMsg('optimlib:ppcgr:ExitLocalOpt','qp_zero_curvature',printMsg,'');
        else % either denom < 0, or denom = 0 and sum(d) ~= 0
            exitflag = -3;
            if denom < 0
                msg = qpTrustRegionMsg('optimlib:ppcgr:ExitUnboundedNegCurv','qp_decrease_indefinitely',printMsg,'');
            else
                msg = qpTrustRegionMsg('optimlib:ppcgr:ExitUnboundedZeroCurv','qp_decrease_indefinitely',printMsg,'');
            end
        end
        % If the Newton step is computed via the direct factorization, then it is
        % "exact" and no actual PCG iterations took place. This is done for display
        % purposes only.
        if isinf(pcf(1))
            cgiters = 0;
        else
            cgiters = k;
        end
        return
    else % denom > 0, continue
        alpha = inner1/denom;
        s = s + alpha*d;
        r = r - alpha*w;
    end
    
    % Precondition
    rr =[zeros(m,1);r];
    rhs = [rr;zeros(m,1)];
    zz = L\(P*rhs);
    zt(pcol,1) = U\zz;
    z = zt(m+1:n,1);
    
    % Exit?
    znrm = norm(z);
    if znrm <= stoptol
        s = y +  PT'*fzmult(A,s,noTranspose,LZ,UZ,pcolZ,PZ);
        if computeLambda
            w = feval(mtxmpy,H,s,varargin{:});
            rhs = -gradsave-w;
            lambda.eqlin = Aeq'\rhs;
        end
        exitflag = 1;
        msg = qpTrustRegionMsg('optimlib:ppcgr:ExitSolved','qp_equality_good',printMsg,'');
        
        % If the Newton step is computed via the direct factorization, then it is
        % "exact" and no actual PCG iterations took place. This is done for display
        % purposes only.
        if isinf(pcf(1))
            cgiters = 0;
        else
            cgiters = k;
        end
        return;
    end
    inner2 = inner1;
    inner1 = r'*z;
    % For next iteration
    beta = inner1/inner2;
    d = z + beta*d;
end
cgiters = k;

if cgiters >= kmax
    exitflag = 0;
    msg = qpTrustRegionMsg('optimlib:ppcgr:ExitMaxIter','solver_stopped_premature',printMsg,'');
end
s = y +  PT'*fzmult(A,s,noTranspose,LZ,UZ,pcolZ,PZ);
if computeLambda
    w = feval(mtxmpy,H,s,varargin{:});
    rhs = -gradsave-w;
    lambda.eqlin = Aeq'\rhs;
end
% If the Newton step is computed via the direct factorization, then it is
% "exact" and no actual PCG iterations took place. This is done for display
% purposes only. This condition is not expected to be met in the case that
% the Newton step is computed via the direct factorization. This is done to
% maintain consistency.
if isinf(pcf(1))
    cgiters = 0;
end