% function [xk,fk,gk] = steepest_descent(x0,f,fpars,opts)
% 
%  Uses steepest descent (BUG: with fixed step size)
%  to find a minimizer of the function f, with starting point x0
%
%  Inputs
%   x0    starting point for the iterations
%   f     handle to function to optimize, see quad.m for sample interface
%   fpars (optional) parameters for the function
%   opts  structure specfiying options for optimization procedure
%    .tol    relative tolerance for convergence (DEFAULT: 1e-8)
%    .maxit  max number of iterations (DEFAULT: 100)
%
%  Outputs
%   xk      minimizer
%   fk, gk  f(xk) and \grad f(xk)
%
% sample usage:
%  opts.tol = 1e-5
%  steepest_descent([1,1],@quad,[],opts)
function [xk,fk,gk] = steepest_descent(x0,f,fpars,opts)

% process arguments
if (nargin<4)
 opts=[];   % default options for the optimization procedure
 if (nargin<3)
  fpars=[]; % default parameters for the function
 end;
end;
if ~isfield(opts,'tol') opts.tol=1e-8; end;   % default tolerance
if ~isfield(opts,'maxit') opts.maxit=100; end; % default max iterations

% calculate the f, gradf at initial point 
xk=x0; k=0; [fk,gk] = f(xk,fpars);

% we look at the norm of the gradient to detect convergence
% the scaling by norm(x0) 
while (norm(gk) >= opts.tol * norm(x0) & k < opts.maxit)
 tk = 1e-1; % compute step size
 xk = xk - tk*gk; k=k+1; % take steepest descent step
 [fk,gk] = f(xk);        % function evaluation at x_{k+1}
 fprintf('k=%4d, f(x) = %15g, ||gradf(x)|| = %15g\n',k,fk,norm(gk));
end;