% cgpdq  -  Conjugate Gradient minimizer for Positive-Definite Quadratics
%  Minimizes function  f(x) = (1/2)x'Qx - a'x;
%  equivalently, solves  Qx = a.
%
%  [x,evalct] = cgpdq(Q,a,x0,eta)
%       Q - symmetric, positive-definite matrix
%       a - column vector, sized to match Q
%      x0 - initial guess at solution
%     eta - desired tolerance on residual Qx-a
%
%       x - returned column vector answer
%  evalct - number of 1D steps taken
%
%  VERBOSITY - global workspace variable, value controls printing levels.
%
%  Theory says that x should be determined exactly once evalct reaches
%  the number of rows in Q, but in practice more steps can make
%  a dramatic decrease in a residual that's already quite small.
%  If user asks for an eta that's too small, routine keeps taking
%  more steps, to a maximum of twice the theoretical limit.  Look
%  at the return variable evalct when it comes back:  if it's smaller
%  than 2n, the tolerance has been achieved; if it equals 2n, the
%  tolerance hasn't been realized, because it was unreasonably small.

% Philip D. Loewen
% 26 Feb 98 - original
%  7 Mar 01 - automatic restart, cosmetics
%  1 May 01 - general beautification

function [x,evalct] = cgpdq(QQ,a,x0,eta)

Q = (QQ+QQ')/2; % Enforce symmetry of $Q$.
x = x0;
xo = x;

v  = Q*x-a;         % Transposed gradient of  f(x)=(1/2)x'Qx-a'x  at point x
rc = v'*v;          % Norm-squared of gradient at "current" point
k  = 0;             % Counts steps away from given point x0.
fo = 0.5*x'*Q*x-a'*x;
fx = 0.5*x'*Q*x-a'*x;

while (rc > eta^2) & (k<2*max(size(Q))),   % Second test is a safeguard.

  vp(1,' ');
  vp(2,['k=',sprintf('%2.0f',k),':  x''=',sprintf(' %11.3e',x'),'.']);

  % Define direction in which to seek new point:
  fulllaps = k/max(size(Q));
  if k==0 | fulllaps==round(fulllaps),
    % Fresh start, assuming no previous direction info.
    vp(1,['k=',sprintf('%2.0f',k),':  Restart--compute residual directly, put d = -grad f(x).']);
    v  = Q*x-a;      % Transposed gradient of  f(x) = x' Q x - a' x  at point x
    rc = v'*v;       % Norm-squared of gradient at "current" point
    d  = -v;
    if rc<=eta^2, break, end;  % Sometimes exact residual is extra small
  else
    b = rc/ro;     % b = coefficient in previous search direction.
    d = -v + b*d;  % New direction: steepest descent plus correction.
  end;
  
  fx = 0.5*x'*Q*x - a'*x;    % Used only in diagnostics.
  vp(1,['k=',sprintf('%2.0f',k),':  Residual norm =',sprintf('%9.1e',sqrt(rc)),'.']);
  if k>0,
    vp(1,['k=',sprintf('%2.0f',k),':  Length of step to this x =', ...
          sprintf('%9.1e',norm(x-xo)),'.']);
    vp(2,['k=',sprintf('%2.0f',k),':  Relative f-reduction (wrt step ',int2str(k-1),') =', ...
          sprintf('%9.1e',(fo-fx)/(eps+abs(fx))),'.']);
    vp(2,['k=',sprintf('%2.0f',k),':  Relative change in x (wrt step ',int2str(k-1),') =', ...
          sprintf('%9.1e',norm(x-xo)/(eps+norm(x))),'.']);
  end

  vp(3,['k=',sprintf('%2.0f',k),':  d'' =',sprintf(' %11.3e',d'),'.']);
  vp(3,['k=',sprintf('%2.0f',k),':  v''d=',sprintf(' %11.3e',v'*d),'.']);

  w = Q*d;         % Here is the one matrix-vector product per iteration.
  t = rc / (d'*w); % In pure quadratic case, step length t known exactly.
  v = v + t*w;     % Analytic update formula for gradient is available too.

  ro = rc;         % Save "old" value of gradient's norm-squared.
  rc = v'*v;
  xo = x;          % Remember previous x-value.
  x = x + t*d;     % Update to new point.

  fo = fx;
  fx =0.5*x'*Q*x-a'*x;
  vp(3,['k=',sprintf('%2.0f',k),':  t  =',sprintf('  %10.4e',t),'.']);

  k = k+1;
end;            % End of main loop.

evalct = k;

vp(1,' ');
fx = 0.5*x'*Q*x-a'*x;
vp(1,['k=',sprintf('%2.0f',k),':  Stopping with x shown below.']);
for j=1:length(x),
  vp(1,sprintf([blanks(ceil(7+log10(k+0.5))),'%19.11e'],x(j)));
end;
vp(1,['k=',sprintf('%2.0f',k),':  Residual norm predicted by CG  =', ...
      sprintf('%9.1e',sqrt(rc)),'.']);
vp(1,['k=',sprintf('%2.0f',k),':  Residual norm from direct calc =', ...
      sprintf('%9.1e',norm(Q*x-a)),'.']);
vp(1,['k=',sprintf('%2.0f',k),':  Length of step to this x =', ...
      sprintf('%9.1e',norm(x-xo)),'.']);
vp(2,['k=',sprintf('%2.0f',k),':  Relative f-reduction (wrt step ',int2str(k-1),') =', ...
      sprintf('%9.1e',(fo-fx)/(eps+abs(fx))),'.']);
vp(2,['k=',sprintf('%2.0f',k),':  Relative change in x (wrt step ',int2str(k-1),') =', ...
      sprintf('%9.1e',norm(x-xo)/(eps+norm(x))),'.']);

vp(1,' ');

return
