function[p,posdef,k] = pcgr(DM,DG,g,kmax,tol,mtxmpy,H,R,pR,callerflag,pcoptions,varargin)
%

%PCGR	Preconditioned conjugate gradients
%
% [p,posdef,k] = PCGR(DM,DG,g,kmax,tol,mtxmpy,H,R,pR) apply
% a preconditioned conjugate gradient procedure to the quadratic
% 
%         q(p) = .5p'Mp + g'p, where
%
% M = DM*H*DM + DG. kmax is a bound on the number of permitted
% CG-iterations, tol is a stopping tolerance on the residual (default
% is tol = .1), mtxmpy is the function that computes products
% with the Hessian matrix H,  
% and R is the cholesky factor of the preconditioner (transpose) of
% M. So, R'R approximates M(pR,pR), where pR is a permutation vector.
% On output p is the computed direction, posdef = 1 implies
% only positive curvature (in M) has been detected; posdef = 0
% implies p is a direction of negative curvature (for M).
% Output parameter k is the number of CG-iterations used (which
% corresponds to the number of multiplications with H).
%

%   Copyright 1990-2022 The MathWorks, Inc.

% Initializations.
n = length(DG);
r = -g; 
p = zeros(n,1); 

% Precondition .
z(pR,1) = R\(R'\r(pR));

stoptol = tol*norm(z);
inner1 = r'*z;
posdef = 1;

kmax = max(kmax,1);  % kmax must be at least 1

switch callerflag
    case 'hessprecon'
        % Nothing
    case 'jacobprecon'
        % Add flag - No transpose
        varargin = [{0},varargin];
    otherwise
        error(message('optimlib:pcgr:InvalidCallingFcn'))
end

% First direction is z
d = z;
% PRIMARY LOOP.
for k = 1:kmax
   ww = DM*d;
   
   % matrix product
   w = feval(mtxmpy,H,ww,varargin{:});
   
   ww = DM*w +DG*d;
   denom = d'*ww;
   if denom <= 0
      if norm(d) == 0
         p = d;
      else      
         p = d/norm(d);
      end
      posdef = 0;
      break
   else
      alpha = inner1/denom;
      p = p + alpha*d;
      r = r - alpha*ww;
   end

   z(pR,1) = R\(R'\r(pR));
   
   % Exit?
   if norm(z) <= stoptol 
      break; 
   end
   inner2 = inner1;
   inner1 = r'*z;
   beta = inner1/inner2;
   d = z + beta*d;
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(pcoptions)
    k = 0;
end
