function [maxJ, err, index] = checkjacobi(fun, fpar, x, h)
%CHECKJACOBI : Check the user's implementation of a nonlinear 
% (vector) function and its Jacobian (gradient).  
% The implementation must be a MATLAB function with declaration
%            function  [f, J] = fun(x, fpar)
% fpar  can e.g. be an array with coordinates of data points, or 
% it may be dummy.  
% f  and  J  should be the vector function and its Jacobian 
% evaluated at  x .
%  
% Call
%      [maxJ, err, index] = checkjacobi(fun, fpar, x, h)
%
% Input parameters 
% fun :  String with the name of the function.
% fpar:  Parameters of the function.  May be empty.
% x   :  The point where we wish to check the Jacobian.
% h   :  Steplength used in difference approximation to J .
%
% Output parameters
% maxJ  :  The maximal absolute value of the elements in  J .
% err   :  Vector with three elements.  The maximal absolute deviation
%          between an element in  J  and the corresponding difference
%          approximation, computed by
%          error(1):  Forward difference with steplength  h
%          error(2):  Backward difference with steplength  h/2 
%          error(3):  "Extrapolated difference", see below. 
% index :  3*2 array, with  index(k,:)  giving the position in  J
%          where  error(k)  occurs. 
%
% If the function is "smooth" and  h  is sufficiently small (but not 
% so small that rounding errors dominate), then we can expect
%       error(1) "=" A*h,   error(2) "=" -A*(h/2)
% where "=" means "approximately equal to" and  A  is a constant, that
% depends on the function and  x , but not on  h .  The "extrapolated
% approximation" is  
%       (forward approximation + 2(backward approximation))/3
% and its error is of the order of magnitude  error(2)^2 .
%
% The algorithm is further discussed in "Checking Gradients", which
% can be downloaded from   http://www.imm.dtu.dk/~hbn/Software/#AUX

% Hans Bruun Nielsen,  IMM, DTU.  99.09.30 / 00.09.13
% Modified to work with scalar functions.  00.11.30

  %  Initialize
  sx = size(x);   n = max(sx);
  if  min(sx) ~= 1,  error('x  must be a vector'), end
  [f J] = feval(fun, x, fpar);                        % offset point
  sf = size(f);   m = max(sf);
  if  min(sf) ~= 1,  error('f  must be a vector'), end
  sJ = size(J);  
  if  m == 1  % Scalar case.  Allow both row and column
    ok = (min(sJ) == 1) & (max(sJ) == n);
    J = J(:).';
  else        % Vector case.  Insist on  m*n
    ok = norm([m n] - sJ) == 0;
  end
  if  ~ok
    error('The sizes of  x , f  and  J  do not match'), end
  
  err = zeros(3,1);   index = zeros(3,2);
  maxJ = max(max(abs(J)));   df = zeros(m,3);

  % loop through the coordinate directions.
  for j = 1 : n 
    %    Make differences for all function components
    xx = x;   xx(j) = x(j) + h;                      % forward point 
    dx = xx(j) - x(j);           % true difference in j'th direction        
    ff = feval(fun, xx, fpar);
    df(:,1) = (ff(:) - f(:))/dx;            % forward approximation
    xx(j) = x(j) - h/2;   dx = x(j) - xx(j);        % backward point       
    fb = feval(fun, xx, fpar);
    df(:,2) = (f(:) - fb(:))/dx;           % backward approximation
    df(:,3) = (2*df(:,2) + df(:,1))/3;  % extrapolated approximation
    %  Store errors
    for  k = 1 : 3
      dif = df(:,k) - J(:,j); 
      [md, i] = max(abs(dif));    % largest difference and its index
      if  md > abs(err(k))                           % New maximum
        err(k) = dif(i);   index(k,:) = [i j];
      end
    end;
  end  % j-loop
% ==========  end of  checkjacobi