function [u,uk,res]=SteepestDescent(A,f,u0,tol,m)
% STEEPESTDESCENT solves Au=f using the steepest descent method
%   [u,uk,res]=SteepestDescent(A,f,u0,tol,m); solves Au=f using steepest 
%   descent starting at the initial guess u0 up to a tolerance tol using 
%   at most m iterations. A has to be symmetric positive definite.
%   SteepestDescent returns in the matrix uk the iterates, in u the 
%   solution computed, and in res the history of the norm of the residuals.

if nargin<5, m=100; end          % default values
if nargin<4, tol=1e-6; end
r=f-A*u0;
res(1)=norm(r);
uk(:,1)=u0;
k=0;
while res(k+1)/norm(f)>tol && k<=m
  k=k+1;
  al=r'*r/(r'*A*r);
  uk(:,k+1)=uk(:,k)+al*r;
  r=f-A*uk(:,k+1);
  res(k+1)=norm(r);
end
u=uk(:,k+1);
