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