I think that matIdx = t(0:(n-1)) @ j(1,k,k) + j; might be a little faster. I wasn't really expecting this, but perhaps it is because you avoid subtracting one from all the elements of row(m)? I tested on a 200*40 matrix and the permutation operation is about 30% quicker.