4 ms·

It would be useful to add a recursively defined algorithm like trapezoidal integration, which is something that basically every autovectorization tool will fall

by qppo 6y ago

It would be useful to add a recursively defined algorithm like trapezoidal integration, which is something that basically every autovectorization tool will fall over y[n] = (x[n] + x[n-1]) / 2 + y[n-1] Which will quickly show you that you can't just replace the numeric operations with their intrinsic counterparts and get speed up (and if you did, you may wind up with incorrect code!). It also shows you that sometimes more is less. To vectorize that code you wind up doing more multiplication than you would normally require for the scalar, and depending on the conditions where you execute it, it might not be able to improve over the scalar implementation. There are more formal methods for converting these kinds of operations (linear, time-invariant recursive systems/filters) to vector equivalents using state space matrix forms, but then you start having to benchmark whether or not you get any benefits.

dragontamer 6y ago
I got nerd-sniped. Lets see... y[n] = (x[n] + x[n-1]) / 2 y = prefixSum(y); I'm pretty sure prefixsum over Y would solve the problem (see: https://en.wikipedia.org/wiki/Prefix_sum#Algorithm_1:_Shorter_span,_more_parallel https://en.wikipedia.org/wiki/Prefix_sum#Algorithm_1:_Shorte... for the Vectorized solution to prefix sum). Without the addition of y[n-1], the first line is now trivially parallelizable by most auto-vectorizers. The prefix-sum probably needs to be written with intrinsics by hand, but once written it generalizes to many operations. (Prefix-max, prefix-min, prefix-XOR). "Prefix-sum" is sometimes called "scan" in the literature (https://developer.nvidia.com/gpugems/gpugems3/part-vi-gpu-computing/chapter-39-parallel-prefix-sum-scan-cuda https://developer.nvidia.com/gpugems/gpugems3/part-vi-gpu-co...) ----------- EDIT: A surprising number of "recursive" definitions can become a prefix-operation (or scan). Associativity is the main requirement. (Note that floating-points are not associative: error builds up non-associatively). So you'd need to redo your error analysis if these were floats.
- Const-me 6y ago
  > the first line is now trivially parallelizable by most auto-vectorizers It’s not that trivial to vectorize, would be something like this for SSE: const __m128 x = _mm_loadu_ps( px ); // Rotate the lanes __m128 xPrev = _mm_shuffle_ps( x, x, _MM_SHUFFLE( 2, 1, 0, 3 ) ); // Insert the highest lane of previousStepX into lowest lane of xPrev // The IMM number means extract lane 3, insert into lane 0, don't zero out any lanes xPrev = _mm_insert_ps( xPrev, previousStepX, 0b11000000 ); const __m128 halfSum = _mm_mul_ps( _mm_add_ps( x, xPrev ), _mm_set1_ps( 0.5f ) ); previousStepX = x; Which compiler do you think would do that automatically?
  - dragontamer 6y ago
    Why would you need to rotate the lanes? const __m128 x = _mm_loadu_ps( px ); const __m128 xPrev = _mm_loadu_ps( px + 4 ); _mm_loadu_ps: load-unaligned. There's no requirement for addresses to load/unload. There are probably penalties for loading across a cache line, but its perfectly valid to do this. ----------- > Which compiler do you think would do that automatically? GCC simply does the unaligned load: https://godbolt.org/z/2HwYMi https://godbolt.org/z/2HwYMi Clang does the shuffle successfully: https://godbolt.org/z/Y3hbQn https://godbolt.org/z/Y3hbQn The shuffle is probably optimal... I'm surprised CLang is good enough to recognize the opportunity. But simply doing 2x loads from memory is probably fine, especially because the 2nd load is guaranteed to be from L1 cache. Besides, assuming 128-bit vectors, 2x 128-bit vectorized loads every 4-elements (averaging 1/2 a load per element) is still superior to scalar code (1-load per element). Honestly, Clang surprised me by generating the superior shuffle code. I was expecting both GCC and Clang to generate the suboptimal (but simpler) unaligned loads.
    - Const-me 6y ago
      > Why would you need to rotate the lanes? To half the count of RAM access instructions. Even L1D cache hit is still slower than these two instructions for shuffle and insert, both are quite fast, they deliver the result on the next cycle. > GCC simply does the unaligned load Right, that’s twice as many loads as necessary. > Clang does the shuffle successfully I wouldn’t call that code “successfully”. Clang uses shufps instructions to shuffle integer vectors. On many CPUs, especially older than a couple of years, there’s a non-trivial latency, a few cycles, to pass data between integer and floating-point execution units. The approach is correct with these shuffles, but not the implementation. For integers it should have used at least _mm_shuffle_epi32 = pshufd or _mm_slli_si128 = pslldq, but ideally _mm_alignr_epi8 = palignr from SSSE3.
      - dragontamer 6y ago
        I guess we have different measures of "success". I don't expect any compiler to beat me in intrinsics or hand-crafted assembly. After all, I can use the compiler. But the compiler can't use my brain. So in all circumstances, I will always at least tie the compiler. Its an "auto-vectorization success", because the pseudo-code I've laid out auto-vectorizes to SSE 128-bit, AVX 256-bit, and AVX512 512-bit code. Furthermore, CLang supports many architectures, and probably can convert my simple for-loop into SIMD on other architectures without any intervention (ie: ARM Neon) The code doesn't need to be optimal, it just needs to be faster than the scalar form and portable.
- qppo 6y ago
  That's an elegant way of doing it. The cool thing about these kinds of algorithms is that there are multiple forms with different numerical and performance behaviors, and it can be tough to pick the best one for the given context. The more classical formal methods are the state-space canonical forms [0] which can be derived for higher order sequences and may be faster than a prefix sum implementation, or might not be, or may be better behaved numerically. [0] https://www.engr.mun.ca/~millan/Eng6825/canonicals.pdf https://www.engr.mun.ca/~millan/Eng6825/canonicals.pdf
cjhanks 6y ago
It is possible to SIMD parallelize multiple linear problems when you can guarantee convergence of each subproblem and compute iterative bounds. But sure, SIMD is no panacea. Especially when the latency between the vector instruction shares transistors with the previous instruction.

by qppo 6y ago

dragontamer 6y ago

I got nerd-sniped. Lets see... y[n] = (x[n] + x[n-1]) / 2 y = prefixSum(y); I'm pretty sure prefixsum over Y would solve the problem (see: https://en.wikipedia.org/wiki/Prefix_sum#Algorithm_1:_Shorter_span,_more_parallel https://en.wikipedia.org/wiki/Prefix_sum#Algorithm_1:_Shorte... for the Vectorized solution to prefix sum). Without the addition of y[n-1], the first line is now trivially parallelizable by most auto-vectorizers. The prefix-sum probably needs to be written with intrinsics by hand, but once written it generalizes to many operations. (Prefix-max, prefix-min, prefix-XOR). "Prefix-sum" is sometimes called "scan" in the literature (https://developer.nvidia.com/gpugems/gpugems3/part-vi-gpu-computing/chapter-39-parallel-prefix-sum-scan-cuda https://developer.nvidia.com/gpugems/gpugems3/part-vi-gpu-co...) ----------- EDIT: A surprising number of "recursive" definitions can become a prefix-operation (or scan). Associativity is the main requirement. (Note that floating-points are not associative: error builds up non-associatively). So you'd need to redo your error analysis if these were floats.

Const-me 6y ago
> the first line is now trivially parallelizable by most auto-vectorizers It’s not that trivial to vectorize, would be something like this for SSE: const __m128 x = _mm_loadu_ps( px ); // Rotate the lanes __m128 xPrev = _mm_shuffle_ps( x, x, _MM_SHUFFLE( 2, 1, 0, 3 ) ); // Insert the highest lane of previousStepX into lowest lane of xPrev // The IMM number means extract lane 3, insert into lane 0, don't zero out any lanes xPrev = _mm_insert_ps( xPrev, previousStepX, 0b11000000 ); const __m128 halfSum = _mm_mul_ps( _mm_add_ps( x, xPrev ), _mm_set1_ps( 0.5f ) ); previousStepX = x; Which compiler do you think would do that automatically?
- dragontamer 6y ago
  Why would you need to rotate the lanes? const __m128 x = _mm_loadu_ps( px ); const __m128 xPrev = _mm_loadu_ps( px + 4 ); _mm_loadu_ps: load-unaligned. There's no requirement for addresses to load/unload. There are probably penalties for loading across a cache line, but its perfectly valid to do this. ----------- > Which compiler do you think would do that automatically? GCC simply does the unaligned load: https://godbolt.org/z/2HwYMi https://godbolt.org/z/2HwYMi Clang does the shuffle successfully: https://godbolt.org/z/Y3hbQn https://godbolt.org/z/Y3hbQn The shuffle is probably optimal... I'm surprised CLang is good enough to recognize the opportunity. But simply doing 2x loads from memory is probably fine, especially because the 2nd load is guaranteed to be from L1 cache. Besides, assuming 128-bit vectors, 2x 128-bit vectorized loads every 4-elements (averaging 1/2 a load per element) is still superior to scalar code (1-load per element). Honestly, Clang surprised me by generating the superior shuffle code. I was expecting both GCC and Clang to generate the suboptimal (but simpler) unaligned loads.
  - Const-me 6y ago
    > Why would you need to rotate the lanes? To half the count of RAM access instructions. Even L1D cache hit is still slower than these two instructions for shuffle and insert, both are quite fast, they deliver the result on the next cycle. > GCC simply does the unaligned load Right, that’s twice as many loads as necessary. > Clang does the shuffle successfully I wouldn’t call that code “successfully”. Clang uses shufps instructions to shuffle integer vectors. On many CPUs, especially older than a couple of years, there’s a non-trivial latency, a few cycles, to pass data between integer and floating-point execution units. The approach is correct with these shuffles, but not the implementation. For integers it should have used at least _mm_shuffle_epi32 = pshufd or _mm_slli_si128 = pslldq, but ideally _mm_alignr_epi8 = palignr from SSSE3.
    - dragontamer 6y ago
      I guess we have different measures of "success". I don't expect any compiler to beat me in intrinsics or hand-crafted assembly. After all, I can use the compiler. But the compiler can't use my brain. So in all circumstances, I will always at least tie the compiler. Its an "auto-vectorization success", because the pseudo-code I've laid out auto-vectorizes to SSE 128-bit, AVX 256-bit, and AVX512 512-bit code. Furthermore, CLang supports many architectures, and probably can convert my simple for-loop into SIMD on other architectures without any intervention (ie: ARM Neon) The code doesn't need to be optimal, it just needs to be faster than the scalar form and portable.
qppo 6y ago
That's an elegant way of doing it. The cool thing about these kinds of algorithms is that there are multiple forms with different numerical and performance behaviors, and it can be tough to pick the best one for the given context. The more classical formal methods are the state-space canonical forms [0] which can be derived for higher order sequences and may be faster than a prefix sum implementation, or might not be, or may be better behaved numerically. [0] https://www.engr.mun.ca/~millan/Eng6825/canonicals.pdf https://www.engr.mun.ca/~millan/Eng6825/canonicals.pdf

cjhanks 6y ago

It is possible to SIMD parallelize multiple linear problems when you can guarantee convergence of each subproblem and compute iterative bounds. But sure, SIMD is no panacea. Especially when the latency between the vector instruction shares transistors with the previous instruction.