4 ms·

Figured I may as well start writing a notes page on Transpose instead of just replying here: https://mlochbaum.github.io/BQN/implementation/primitive/transpose

Figured I may as well start writing a notes page on Transpose instead of just replying here: https://mlochbaum.github.io/BQN/implementation/primitive/transpose.html https://mlochbaum.github.io/BQN/implementation/primitive/tra... (EDIT: missed a 0 in the first version and got the wrong timing for Dyalog) For the particular case you give, the axis permutation swaps two axes, so it's invertible and the same left argument applies in both J and APL (may require adjustment for ⎕IO of course). I measure Dyalog the same as J, 0.3s. This array has a large fixed axis at the end, meaning that chunks of 77 floats (over half a kB) stay together, and the best way is just to move each of these with memcpy. Seems both languages miss this optimization? 0.3s is only 2.2GB/s. The array's definitely too large to fit in even L3 cache but this seems too slow for an in-memory copy.

mlochbaum 4y ago
No, I think they're both using memcpy. perf says they're spending all the time in libc at least. The total time is a little slower than operations like rotate that I know are using memcpy (0.18s), so the chunk size does seem to introduce some overhead.
- moonchild 4y ago
  Glibc strings functions are ok, but not great (but one-size-fits-all is hard, and for all I know they could be fine here). While 512 bytes seems like a lot, the core loop is likely 4x unrolled avx; that is 128 bytes per iteration, so only 4 iterations. And it probably tries to align the dst, so annoying dispatch overhead (which you could avoid by working in batches of k at a time, eg k=4 for double floats and avx in the general case).
  - mlochbaum 4y ago
    The big difference is the access pattern: see the benchmarks below. Index does the small memcpys, but it speeds up if the indices are in order (my earlier benchmark used ⌽, not ⊖, because I don't remember APL any more). So prefetching might help. I guess it's possible that a larger-scale blocking would too? But 5GB/s for ⊖ isn't great either. In an application that uses these huge arrays and needs the best performance (most don't!), it needs to be split up, ideally into sections that fit in L1, so that multiple array operations can be applied to those chunks and stay CPU-bound (at least for transpose, maybe not for arithmetic). That's why I wouldn't be too interested in optimizing this case. ⎕IO←0 ⋄ ar←,[0 1 2]a←?27 1000 40 77⍴0 ci←(⊢⍴∘⍳×/)3↑⍴a ⍝ cell indices i←⍉ci ⋄ cmpx '2 1 0 3⍉a' 'ar[i;]' 2 1 0 3⍉a → 2.8E¯1 | 0% ⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕ ar[i;] → 2.8E¯1 | 0% ⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕ i←5⊖ci ⋄ cmpx '5⊖a' 'ar[i;]' 5⊖a → 1.2E¯1 | 0% ⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕ ar[i;] → 1.7E¯1 | +46% ⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕
moonchild 4y ago
Thanks. > cache associativity This one is annoying. Nt accesses might help, though probably not for the in-place case. You mention padding; the paper I linked uses padding as well, to shorten cycle length for inplace+parallel. Although, I would guess that dimensions which are nice to the cache will be not-nice to the blocking, and vice versa > Of course larger loops are possible as well. But most of the time it's really the base case that's important. The base case also handles two axes most of the time, but can incorporate all the 2D optimizations like blocking. Unless I am misunderstanding--reversing the axes of an array of shape 1e6 1e3 2 2 will have an inner loop dealing with the two innermost axes, and some extra dispatch outside of that. Which is not ideal, obviously--the inner loop only deals with 4 elements.

mlochbaum 4y ago

No, I think they're both using memcpy. perf says they're spending all the time in libc at least. The total time is a little slower than operations like rotate that I know are using memcpy (0.18s), so the chunk size does seem to introduce some overhead.

moonchild 4y ago
Glibc strings functions are ok, but not great (but one-size-fits-all is hard, and for all I know they could be fine here). While 512 bytes seems like a lot, the core loop is likely 4x unrolled avx; that is 128 bytes per iteration, so only 4 iterations. And it probably tries to align the dst, so annoying dispatch overhead (which you could avoid by working in batches of k at a time, eg k=4 for double floats and avx in the general case).
- mlochbaum 4y ago
  The big difference is the access pattern: see the benchmarks below. Index does the small memcpys, but it speeds up if the indices are in order (my earlier benchmark used ⌽, not ⊖, because I don't remember APL any more). So prefetching might help. I guess it's possible that a larger-scale blocking would too? But 5GB/s for ⊖ isn't great either. In an application that uses these huge arrays and needs the best performance (most don't!), it needs to be split up, ideally into sections that fit in L1, so that multiple array operations can be applied to those chunks and stay CPU-bound (at least for transpose, maybe not for arithmetic). That's why I wouldn't be too interested in optimizing this case. ⎕IO←0 ⋄ ar←,[0 1 2]a←?27 1000 40 77⍴0 ci←(⊢⍴∘⍳×/)3↑⍴a ⍝ cell indices i←⍉ci ⋄ cmpx '2 1 0 3⍉a' 'ar[i;]' 2 1 0 3⍉a → 2.8E¯1 | 0% ⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕ ar[i;] → 2.8E¯1 | 0% ⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕ i←5⊖ci ⋄ cmpx '5⊖a' 'ar[i;]' 5⊖a → 1.2E¯1 | 0% ⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕ ar[i;] → 1.7E¯1 | +46% ⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕

moonchild 4y ago

Thanks. > cache associativity This one is annoying. Nt accesses might help, though probably not for the in-place case. You mention padding; the paper I linked uses padding as well, to shorten cycle length for inplace+parallel. Although, I would guess that dimensions which are nice to the cache will be not-nice to the blocking, and vice versa > Of course larger loops are possible as well. But most of the time it's really the base case that's important. The base case also handles two axes most of the time, but can incorporate all the 2D optimizations like blocking. Unless I am misunderstanding--reversing the axes of an array of shape 1e6 1e3 2 2 will have an inner loop dealing with the two innermost axes, and some extra dispatch outside of that. Which is not ideal, obviously--the inner loop only deals with 4 elements.