function ds = dimsum(A,d,varargin);
% Matrix multiply a vector using a multi-dimensional array.
% This allows weighted summation in a 3-D/4-D arrays over
% a single index.
%
% >> size(T)
% ans =
%   128 64 10 23
% >> size(dz)
% ans =
%   10 1
% >> vintegral=dimsum(T,dz,3);
% >> size(vintegral)
% ans =
%   128 64 23

error(nargchk(2,3,nargin))

if size(d,2) ~= 1, error('Second argument must be a vector'); end

if nargin==3
  dim=varargin{1};
else
  if sum( size(A)==size(d,1) )==1
   dim=find( size(A)==size(d,1));
  else
   error('One and only one dimension must match unless you specify the dimension to sum')
  end
end

if size(A,dim) ~= size(d,1)
 error('Dimensions of matrix and vector do not match')
end

dims=[size(A) 1];
nd=size(d,1);
ip=[dim:ndims(A)+1 1:dim-1];
A=permute(A,ip);
A=A(1:nd,:);
ds=A'*d;

dims(dim)=1;
dims=dims(ip);
ds=reshape(ds,dims);
ds=squeeze(ipermute(ds,ip));