C = tmult(A, B) is equivalent to:
sz = [size(B) 1];
sz(1) = size(A, 1);
C = zeros(sz);
for a = 1:prod(sz(3:end))
C(:,:,a) = A(:,:,a) * B(:,:,a);
end

but is completely vectorized, so much faster. Tmult also supports bsxfun-style expansion of singular dimensions where appropriate, such that tmult(rand(4, 3, 10), rand(3, 2)) yields a 4x2x10 output.