function [x,rnorms,errornorms] = gradient_descent(A,b,tol,maxit,x0,plotphifun)
%
% Gradient Descent algorithm to minimize 0.5 x'Ax - b'x, assuming that
%  A is symmetric and positive definite.  Solution is x = A^{-1} b.
% Could be useful if cannot store A, but only have access to matrix-vector
% multiplies, but this is just for illustrative purposes: the 
% Conjugate Gradient method is almost always much better in this situation.
% However, Gradient Descent generalizes to non-quadratic functions, while
% it is not clear how to best generalize CG to non-quadratic functions.
%
% tol, maxit, x0, plotphifun are all optional. Setting plotphifun to 1 is
% very expensive!  
%
% rnorms is the vector of residual norms
% errornorms is the vector of error norms which requires solving xsol=A\b
%  to compute, so this is expensive!
%
% It is interesting to run this with A set to the discrete Laplacian, and
% to compare the rate of convergence of rnorms and errornorms to the
% formula (kappa-1)/(kappa+1) where kappa is the condition number of A.
% Not so interesting for A = B'*B, with B random, which is very ill conditioned.
%
if nargin < 3
    tol = 1e-6;
    
end
if nargin < 4
    maxit = 10*length(b);
end
if nargin < 5
    x0 = zeros(size(b));
end
if nargin < 6
    plotphifun = 0;
end
%
x = x0;
r = b - A*x;
rtr = r'*r;
iter = 0;
if nargout >= 2
    rnorms = []; % so defined if loop does not execute
end
if nargout >= 3
    errornorms = [];
    xsol = A\b; % just so we can compute the error
end
while rtr > tol^2 && iter < maxit
    iter = iter + 1;
    s = A*r;
    rtAr = r'*s; %  WHY SHOULD THIS BE POSITIVE? SHOULD WE CHECK THIS?
    alpha = rtr/rtAr; % minimizes phi along x + alpha r: see p.183
    if plotphifun  % expensive! but may give insight
        getphifunplot(A,b,x,r,alpha)
    end
    x = x + alpha*r; % new x iterate
    %%%  r = b - A*x;  % definition of new residual
                  % there is another way of writing the new residual
                  % we have rnew = b - A*xnew 
                  %              = b - A*(x + alpha*r)
                  %              = b - A*x - alpha*A*r
                  %              = r - alpha*A*r 
                  %              = r - alpha*s
                  % which SAVES ONE MATRIX-VECTOR MULTIPLY.... try it.
                  % This trick is used in the CG algorithm too.
    r = r - alpha*s;
    rtr = r'*r;
    if nargout >= 2
        rnorms(iter) = sqrt(rtr);
    end
    if nargout >= 3
        errornorms(iter) = norm(x - xsol);
    end
end

end % end of function gradient_descent

function getphifunplot(A,b,x,r,alpha)
%
m = 101; % plot phi at m points from 0 to 2*alpha
t = linspace(0,2*alpha,m);
for i=1:length(t)
    y = x + t(i)*r;
    phifun(i) = 0.5*y'*A*y - b'*y;
end
clf
plot(t,phifun)
hold on
plot(0,phifun(1),'k*')
plot(alpha,phifun((m-1)/2),'r*'); % plot the function at t=alpha
legend('phi','current point','new iterate')
shg
fprintf('press any key to continue\n')
pause

end % end of function getphifunplot