function [x,rnorm,enorm]=algSD(A,b,rtol,maxit,x0,xex)
%Steepest Descent Algorithm.  Algorithm 5.2 page 138
% input:
%   A     = sparse nxn (symmetric, positive definite)
%   b     = rhs nx1
%   rtol  = stopping criteria for the relative residual, ||r||/||r_0||
%   maxit = max # of iterations allowed
%   x0    = init guess nx1
%   xex   = exact solution, A\b, nx1
%
% nonvargin: A,b
%
% default:
%   rtol  = 1e-6
%   maxit = 100
%   x0    = rand(size(b));
%   xex   = []
%
% usage:
%   n=100;
%   A=gallery('poisson',n,n); b=zeros(n,1); x0=rand(n,1);
%
%   x=algSD(A,b)
%   x=algSD(A,b,rtol)
%   x=algSD(A,b,rtol,maxit)
%   x=algSD(A,b,rtol,maxit,x0)
%   x=algSD(A,b,rtol,maxit,x0,xex)
%   [x,r]=algSD(A,b,...)
%   [x,r,e]=algSD(A,b,...)
%

% Author: L. Olson
%       Initial: Fri 26 Jan 2007 04:22:32 PM CST
% Last Modified: 01/29/2007 07:02:55

n=length(b);

if nargin<3 || isempty(rtol)
  rtol=1e-6;
end
if nargin<4 || isempty(maxit)
  maxit=100;
end
if nargin<5 || isempty(x0)
  x0=rand(n,1);
end
if nargin<6 || isempty(xex)
  xex=zeros(size(b));
end

% initialize
x=x0;

r=b-A*x;
rnorm=zeros(n,1);
enorm=zeros(n,1);
rnorm(1)=norm(r);
enorm(1)=1.0;

%-----------------------------------------
% main loop
for iter=1:maxit
  p=A*r;
  alph=dot(r,r)/dot(p,r);
  x=x+alph*r;
  r=r-alph*p;
  %
  rnorm(iter)=norm(r);
  if ~isempty(xex)
    enorm(iter)=norm(x-xex);
  end
  rnorm(iter)=norm(r);
  if (rnorm(iter)/rnorm(1) < rtol)
    break;
  end
end
%-----------------------------------------

rnorm=rnorm(1:iter);
enorm=enorm(1:iter);