function dmp = mpder(mp,xi)
%MPDER derivative of a multivariate polynomial in mp-form
%
%        DMP = MPDER(MP,XI)
%
% applies the directional derivatives D_XI(:,j) to the polynomial in MP.
%
%   See also: MPMAK

% cb: 9aug97; regress to v4 13sep97
%   Copyright (c) 1997 by C. de Boor

% MP contains (see MPMAK), among other things, the COEFS for the representation
% 
%   p(x) = sum_a x^a (|a| choose a) COEFS(nn(a)) 
%
% Differentiation in the direction  xi  gives  
%   D_xi ()^a/a! = sum_i xi(i) ()^{a-e_i}/(a-e_i)!
% hence
%
%  D_xi p(x) = sum_a x^a (|a| choose a) (|a|+1) \sum_i xi(i)*COEFS(nn(a+e_i))

[coefs,d,k,center,dimpjd,mm] = mpbrk(mp);

[dxi,nxi] = size(xi);
if dxi~=d
   error(['The direction(s) must be in R^',num2str(d),'.']), end

[df,nc] = size(coefs);
for ixi=1:nxi
%v5   x = reshape(xi(:,ixi),[1,d,1]);
   x = xi(:,ixi);
   if k==0
      coefs(:,1) = zeros(df,1);
   else
%v5      for j=0:k-1
%v5         ranj = dimpjd(1+j)+1:dimpjd(2+j); lj = length(ranj);
%v5         coefs(:,ranj) = (j+1)*reshape(sum( ...
%v5                         reshape(coefs(:,dimpjd(j+2)+mm(:,ranj)),[df,d,lj]) ...
%v5                                 .*x(ones(df,1),:,ones(1,lj)),2),[df,lj]);
%v5      end
      for ff=1:df
         for j=0:k-1
            ranj = dimpjd(1+j)+1:dimpjd(2+j); lj = length(ranj);
            coefs(ff,ranj) = (j+1)*sum( ...
                            reshape(coefs(ff,dimpjd(j+2)+mm(:,ranj)),d,lj) ...
                                    .*x(:,ones(1,lj)));
         end
      end
      mm(:,dimpjd(1+j)+1:dimpjd(k+2)-lj) = []; % adjust MM to the lowered degree
      k = k-1;
   end
end

dmp = mpmak(coefs(:,1:dimpjd(k+2)),d,k,center,dimpjd(1:k+2),mm);