function[H] = sfd(x,grad,H,group,alpha,funfcn,varargin)
%SFD	Sparse Hessian via finite gradient differences
%
% H = sfd(x,grad,H,group,fdata,fun) returns the
% sparse finite difference approximation H of a Hessian matrix
% of function 'fun'  at current point x.
% Vector group indicates how to use sparse finite differencing:
% group(i) = j means that column i belongs to group (or color) j.
% Each group (or color) corresponds to a finite gradient difference.
% fdata is a data array (possibly) needed by function 'fun'.
%
% H = sfd(x,grad,H,group,fdata,fun,alpha) overrides the default
% finite differencing stepsize.
% 

%   Copyright 1990-2000 The MathWorks, Inc.
%   $Revision: 1.3 $  $Date: 2000/06/16 22:25:42 $

xcurr = x(:);  % Preserve x so we know what funfcn expects
[m,n] = size(H); 
v = zeros(n,1);
ncol = max(group); epsi = sqrt(eps);
if isempty(alpha), alpha = ones(ncol,1)*sqrt(eps); end
H = spones(H); d = zeros(n,1);
for k = 1:ncol
   d = (group == k);
   if nargin < 5
      xnrm = norm(xcurr(d));
      xnrm = max(xnrm,1);
      alpha(k) = alpha(k)*xnrm;
   end
   y = xcurr + alpha(k)*d;
   
   % Make x conform to user-x
   x(:) = y;
   %[dummy,v] = feval(fun,y);
   switch funfcn{1}
   case 'fun'
      error('should not reach this')
   case 'fungrad'
      [dummy,v(:)] = feval(funfcn{3},x,varargin{:});
      %OPTIONS(11)=OPTIONS(11)+1;
   case 'fun_then_grad'
      %newval = feval(funfcn{3},x,varargin{:}); 
      v(:) = feval(funfcn{4},x,varargin{:});
      % OPTIONS(11)=OPTIONS(11)+1;      
   otherwise
      error('Undefined calltype in FMINUNC');
   end
   
   
   w = (v-grad)/alpha(k);
   cols = find(d); 
   lpoint = length(cols);
   A = sparse(m,n);
   A(:,cols) = H(:,cols);
   H(:,cols) = H(:,cols) - A(:,cols);
   [i,j,val] = find(A);
   [p,ind] = sort(i);
   val(ind) = w(p);
   A = sparse(i,j,full(val),m,n);
   H = H + A;
end
H = (H+H')/2;   % symmetricize



