function [u,g,f_vals] = grad_method(fname,u0,options)
% Usage: [u,g,f_vals] = grad_method(fname,u0,options)
% Purpose: example optimization routine for solving the problem
%   min_u fname(u).
%
% Input:
%   fname: name of file for computing function value (and the gradient).
%   u0: initial value for the solution
%   options: structure of user given values for
%     maxits: maximum number of iterations
%     tol: stopping tolerance
%     stepsize: fixed stepsize
%     grad: 'on' = fname returns the gradient also, 'off' = use FD gradient
% Output:
%   u: solution of the optimization problem
%   g: gradient of fname at u
%   f_vals: values of fname during the optimization iterations
% Author: TK 3.11.2001 (late in the dark evening...)
%
maxiters = options.maxits;
tolerance = options.tol;
ustepsize = options.stepsize;
usergradient = options.GradObj;
%
n = length(u0);
u = u0;
f = feval(fname,u);
f_vals = [f];
d = ones(size(u0));
%
k = 0;
while (k < maxiters) & (ustepsize*norm(d) > tolerance)
  k = k + 1;
  switch (usergradient)
    case ('off') % use the finite difference gradient
      g_old = fd_grad(fname,u);
      % Testing the accuracy of FD gradient (just because we can do it here):
      [f,g] = feval(fname,u);
      norm(g - g_old)
    otherwise
      [f,g_old] = feval(fname,u);
  end%switch
  d = - g_old;
  u = u + ustepsize*d;
  f = feval(fname,u)
  if (f_vals(k) < f)
    disp('Function value increasing...')
  end%if
  f_vals = [f_vals; f];
  pause
end
%
[f,g] = feval(fname,u);
