order - fuse - order kernel logic#74
Conversation
Codecov Report✅ All modified and coverable lines are covered by tests.
🚀 New features to boost your workflow:
|
| g = 8 * sizeof(Int) - leading_zeros(M + 1) | ||
| importance = 2 .* (1 .<< (g .* (N .- indexorder(strides[1])))) | ||
| for k in 2:M | ||
| importance = importance .+ (1 .<< (g .* (N .- indexorder(strides[k])))) | ||
| end |
There was a problem hiding this comment.
I have to admit that I forgot exactly how this cost model works 😄
There was a problem hiding this comment.
Something like, the importance of a dimension is
sum over arrays of 2^(log2(M+2) * (N - indexorder of dimension in array))
So if an index is the first index in each of the arrays, its importance is
(M+1) * (M+2)^((N-1))
There was a problem hiding this comment.
I also didn't really look into this in detail. I think realistically, any stride that goes beyond 2 cacheline widths probably doesn't actually matter all that much, so the main part is that we want as many contiguous writes and then contiguous reads?
To be totally honest, I was looking at the permutation part of this last week and since I have a machine with 512b SIMD registers, it's actually possible to speed up this kernel by quite a lot by doing an in-register transpose, for example reading in 8 lanes of 8 contiguous Float64s, then SIMD shuffling them so they're transposed, and then writing the 8 lanes of 8 contiguous floats. The Julia compiler does not seem to be autovectorizing this.
The caveat is of course that you need 2 stride-1 dimensions that are large enough, the exact kernel sizes have to be tuned for the machines, and that this works mainly for larger SIMD widths, as the benefit decreases for larger eltypes or smaller lanewidths.
Anyways, might at some point try and include such a kernel, but probably don't want to do this in Julia, that seems a compile-time nightmare :)
There was a problem hiding this comment.
Oh yes, the strided kernels are far from optimal, especially for very specific tasks like transpose. At some point I wanted to invest time to use LoopVectorization for writing a dedicated add_transpose implementation, but never got to it. Given the current states of LoopVectorization, it might not have been the best idea anyway.
Closes #72.
The previous approach was
fuse-order, where thefuseonly works if at least one of the input arrays is ordered in increasing strides. The pathological case beingpermutedims(A, (2,1)) .= permutedims(B, (2, 1)), where fusion is trivial but never happened.Here, I simply first order, then fuse, and then finally order again to make sure the leftover dim 1's from the fusion are migrated to the end.