  function ShowSteepestDescent()
% Illustrate how steepest descent with exact line search can be used
% to solve the symmetric positive definite system Ax = b by
% minimizing (x'Ax)/2 - x'b

% Generate an example...
A = delsq(numgrid('C',35));
[n,n] = size(A);
b = randn(n,1);

% Take this many steps at most
m = 100;

% No preconditioner...
[xStar,resVec] = SD(A,b,m);
% Incomplete Cholesky predocndioner with different drop
% tolerances...
[xStar0,resVec0] = PSD(A,b,m,cholinc(A,'0')');
[xStar1,resVec1] = PSD(A,b,m,cholinc(A,.01)');
[xStar2,resVec2] = PSD(A,b,m,cholinc(A,.005)');
[xStar3,resVec3] = PSD(A,b,m,cholinc(A,.001)');
close all
semilogy((1:m)',resVec,(1:m)',resVec0,(1:m)',resVec1,(1:m)',resVec2,(1:m)',resVec3)
legend('G = I', 'cholinc(A,''0'')', 'cholinc(A,.01)','cholinc(A,.005)','cholinc(A,.001)','location','best')
title(sprintf('Steepest Descent with Various Preconditioners    n = %1d',n),'fontsize',14)
xlabel('Iterations','fontsize',14)
ylabel('|| A*xk - b ||','fontsize',14)
shg



  function [xStar,resVec] = SD(A,b,itMax)
% A is nxn symmetric positive definite.
% b is nx1.
% Carries out itMax steps steepest descent applied to phi(x) = (x'Ax)/2 - x'*b.
% The starting vector is x = zeros(n,1)
% resVec is a column m-vector, resVec(k) = norm(A*xx - b) where xk is the
%    kth iterate
% xStar is the final iterate.
n = length(b);
x = zeros(n,1);
resVec = zeros(itMax,1);
for k=1:itMax
    % The gradient...
    g = A*x-b;
    % The step length...
    mu = (g'*g)/(g'*A*g);
    % The update...
    x = x - mu*g;
    resVec(k) = norm(A*x - b);
end
xStar = x;


  function [xStar,resVec] = PSD(A,b,itMax,G)
% A is nxn symmetric positive definite.
% b is nx1.
% G is an approximate Cholesky factor of  A.
% Carries out itMax steps of steepest descent applied to 
%      phi(z) = (z'Atilde z)/2 - z'*btilde
% where Atilde = inv(G)*A*inv(G') and btilde = inv(G)*b. Notice that if
% Atilde*z = btilde, then Ax = b provided G'x = z.
% The starting vector is z = zeros(n,1)
% resVec is a column m-vector, resVec(k) = norm(A*xk - b) where xk = G'\zk
% where zk is the kth iterate
% xStar = G'\zStar where zStar is the final z iterate.
n = length(b);
z = zeros(n,1);
btilde = G\b;
errVec = zeros(itMax,1);
for k=1:itMax
    % The gradient...
    g =  G\(A*(G'\z)) - btilde;
    % The step length...
    gtilde = G'\g;
    mu = (g'*g)/(gtilde'*A*gtilde);
    % The update...
    z = z- mu*g;
    resVec(k) = norm(A*(G'\z) - b);
end
xStar = G'\z;
