function y = do_vecdim (hfcn, x, vecdim)

  ## Input validation
  ## Remove dimensions larger than actual array
  vecdim(vecdim > ndims (x)) = [];
  if (isempty (vecdim))
    y = x;
    return;
  endif

  ## Calculate permutation vector and permute
  remdims = 1:ndims (x);    # all dimensions
  remdims(vecdim) = [];     # delete ones specified by vecdim to leave remaining
  perm = [remdims, vecdim];

  y = permute (x, perm);

  ## Reshape to put all vecdims in final dimension
  szy = size (y);
  nremd = numel (remdims);
  sznew = [szy(1:nremd), prod(szy(nremd+1 : end))];
  y = reshape (y, sznew);

  ## Call function handle on single, squashed dimension
  y = hfcn (y, nremd+1);

  ## Inverse permute back to correct dimensions
  y = ipermute (y, perm);

endfunction
