function grdd = grad(f,x0,para);
%
% Purpose: To calculate the gradient of the function f at the point x0
%
% Format: 
%
% If f takes parameters, then the format is
%
% gradient = grad('f',x0,para)
%
% If not, the format is
%
% gradient = grad('f',x0)
%
% where f is a function and f(x) is an mxn matrix
%
% x0 is a kx1 vector
%
% Output: the gradient, an mnxk matrix defined as in Magnus & Neudecker (1988).
%
%

if ~isreal(x0);
   if abs(imag(x0))>eps;
      error('Not implemented for complex matrices');
   else
      x0 = real(x0);
   end
end
if nargin>2;
   f0 = vec(feval(f,x0,para));
else
   f0 = vec(feval(f,x0));
end
n = size(f0,1);
k = size(x0,1);
grdd = zeros(n,k);

ax0 = abs(x0);
if ~(x0==0);
   dax0 = x0./ax0;
else
   dax0 = 1;
end
dh = 10^(-8)*max([ax0,(1e-2)*ones(k,1)]')'.*dax0;
xdh = x0+dh;
dh = xdh-x0;
res = x0*ones(1,k);
arg = diag(xdh) + res - diag(diag(res));
for i = 1:k;
   if nargin>2
      grdd(:,i)=vec(feval(f,arg(:,i),para));
   else
      grdd(:,i)=vec(feval(f,arg(:,i)));
   end
end
grdd = (grdd-f0*ones(1,k))./(ones(n,1)*dh');