5 ms·

How does Pythran fare against the code in this link: https://discourse.julialang.org/t/i-just-decided-to-migrate-from-python-fortran-to-julia-as-julia-was-fast

How does Pythran fare against the code in this link: https://discourse.julialang.org/t/i-just-decided-to-migrate-from-python-fortran-to-julia-as-julia-was-faster-in-my-test/61188/9 https://discourse.julialang.org/t/i-just-decided-to-migrate-... The author made some julia performance mistakes in his benchmarks. The code in this version: https://discourse.julialang.org/t/i-just-decided-to-migrate-from-python-fortran-to-julia-as-julia-was-faster-in-my-test/61188/12 https://discourse.julialang.org/t/i-just-decided-to-migrate-... should be very tricky to beat, and I suspect that the only way is doing a lot of customized assembly.

cycomanic 5y ago
That is the code I'm comparing against. So the fastest Julia version is actually the one using @tturbo in your link. I'm not surprised, pythran gets me the best performance across numba, customized cython and (for some comparisons) julia in most of my tests.
- eigenspace 5y ago
  Ah very interesting, this is great to know! Sorry I didn't properly read your initial comment. I'm now quite interested to see what we can do to match or beat Pythran here. What CPU are you using? I know my CPU, a Zen+ Ryzen 5 2600, is one of the worst modern CPUs for this sort of thing and there's much better scaling on Intel CPUs or the Zen 3 CPUs, but those things are probably also effecting the Pythran code as well. Is Pythran any good at producing BLAS microkernels?
  - celrod 5y ago
    Would be good to confirm the number of threads being used by both, but it could also be that pythran is using a better `sincos` implementation. On very old hardware (which it would have to be if multithreaded), the exp implementation might not be great either. That'd especially be the case if it has no AVX, as the exp implementation relies on a table (which requires gather instructions). For 2000x2000, their time was 14ms, vs 6ms for the 2600.
    - eigenspace 5y ago
      Yeah, I tried to get Pythran working on my computer to compare against, but it wouldn't build for me.
    - cycomanic 5y ago
      Just an additional comment, using sincos in the comparison is not really fair, because it is a change to the algorithm (it assumes that a is real), because in general exp(x) = exp(x.real())+ sincos(x.imag()) (for simplification sincos(x) returns cos(x) + jsin(x), but you get the point). The other code can't easily optimise for that (although it would be a nice optimisation). Interestingly (and I can't explain why, maybe someone knowing more about julia), if I change the code of that third method using tturbo to use the cis function (or exp) it becomes a factor of 20 slower and is actually the slowest of all the different algorithms)
      - celrod 5y ago
        exp(x) = exp(x.real()) * sincos(x.imag()) Anyway, it is the same algorithm as long as `x` is real. The reason for manually inlining exp(::Complex) and manually decomposing the complex number into its real and imaginary parts is that `@tturbo` currently only supports real inputs. Given types that it does not support, it will silently run the original loop single-threaded and optimized simply with `@inbounds @fastmath`. I imagine in changing to cis or exp, you also switched to using complex numbers? If so, that would explain the dramatic slowdown. If not, mind sharing your example, here or in an issue at LoopVectorization.jl? It's planned that "some day" LoopVectorization will support complex numbers through the AbstractInterpreter interface, but this is probably still a ways off.
        cycomanic 5y ago
        I agree, that it's the same algorithm as long as it x is real. My point was that with this optimisation we give more information to this function than the comparable python function (or the other julia implementations), so it gains an advantage for optimising. Yes I was using complex inputs when switching to to cis/exp. LoopVectorization not supporting complex does explain this. As a side note, I generally dislike silent fall-backs, for optimisation libraries/modules. That's why I always do python=False when using numba (and I wish it was default). I'm using the library/module because I want speed up, the problem with the fall-backs is they are often slower than the default code (e.g. unwrapped loops in numba, and in this example), so I would like to know if it doesn't speed up the code, because rather not use it then. A failure would save me benchmarking the slow-down.
- celrod 5y ago
  Did you start Julia with multiple threads? Julia unfortunately starts with only a single thread by default, and the slow time you reported (14 ms) makes it appear that this may have been the case. You can start `julia -t4` for 4 threads, for example. Or you can set the environmental variable `JULIA_NUM_THREADS=4`. If you were using the same number of threads for both, mind letting me know which CPU you're using?
  - cycomanic 5y ago
    You were correct, I just tested with julia -t4 and that improved times by about a factor of 3. N = 2000 27.137 ms (24 allocations: 61.04 MiB) 49.496 ms (3 allocations: 61.05 MiB) 4.597 ms (2 allocations: 61.04 MiB) The CPU is a ryzen 3600 btw.
    - celrod 5y ago
      -t6 would give 6 threads. FWIW, the Ryzen 2600 with 6 threads took 6 ms, so the 3600's floating point and vector improvement is quite evident from your timing. (The 2600 emulates 256-bit operations with 2x 128-bit, while the 3600 has full support.)
      - cycomanic 5y ago
        Yes, I've also run with 6 threads and get an additional minor improvement on the loopvectorized code (3.9 vs 4.5ms). I have to say LoopVectorization.jl is quite an impressive piece of work. I have to try out some more julia at some point, unfortunately pretty much everything I do is using complex values.

cycomanic 5y ago

That is the code I'm comparing against. So the fastest Julia version is actually the one using @tturbo in your link. I'm not surprised, pythran gets me the best performance across numba, customized cython and (for some comparisons) julia in most of my tests.

eigenspace 5y ago
Ah very interesting, this is great to know! Sorry I didn't properly read your initial comment. I'm now quite interested to see what we can do to match or beat Pythran here. What CPU are you using? I know my CPU, a Zen+ Ryzen 5 2600, is one of the worst modern CPUs for this sort of thing and there's much better scaling on Intel CPUs or the Zen 3 CPUs, but those things are probably also effecting the Pythran code as well. Is Pythran any good at producing BLAS microkernels?
- celrod 5y ago
  Would be good to confirm the number of threads being used by both, but it could also be that pythran is using a better `sincos` implementation. On very old hardware (which it would have to be if multithreaded), the exp implementation might not be great either. That'd especially be the case if it has no AVX, as the exp implementation relies on a table (which requires gather instructions). For 2000x2000, their time was 14ms, vs 6ms for the 2600.
  - eigenspace 5y ago
    Yeah, I tried to get Pythran working on my computer to compare against, but it wouldn't build for me.
  - cycomanic 5y ago
    Just an additional comment, using sincos in the comparison is not really fair, because it is a change to the algorithm (it assumes that a is real), because in general exp(x) = exp(x.real())+ sincos(x.imag()) (for simplification sincos(x) returns cos(x) + jsin(x), but you get the point). The other code can't easily optimise for that (although it would be a nice optimisation). Interestingly (and I can't explain why, maybe someone knowing more about julia), if I change the code of that third method using tturbo to use the cis function (or exp) it becomes a factor of 20 slower and is actually the slowest of all the different algorithms)
    - celrod 5y ago
      exp(x) = exp(x.real()) * sincos(x.imag()) Anyway, it is the same algorithm as long as `x` is real. The reason for manually inlining exp(::Complex) and manually decomposing the complex number into its real and imaginary parts is that `@tturbo` currently only supports real inputs. Given types that it does not support, it will silently run the original loop single-threaded and optimized simply with `@inbounds @fastmath`. I imagine in changing to cis or exp, you also switched to using complex numbers? If so, that would explain the dramatic slowdown. If not, mind sharing your example, here or in an issue at LoopVectorization.jl? It's planned that "some day" LoopVectorization will support complex numbers through the AbstractInterpreter interface, but this is probably still a ways off.
      - cycomanic 5y ago
        I agree, that it's the same algorithm as long as it x is real. My point was that with this optimisation we give more information to this function than the comparable python function (or the other julia implementations), so it gains an advantage for optimising. Yes I was using complex inputs when switching to to cis/exp. LoopVectorization not supporting complex does explain this. As a side note, I generally dislike silent fall-backs, for optimisation libraries/modules. That's why I always do python=False when using numba (and I wish it was default). I'm using the library/module because I want speed up, the problem with the fall-backs is they are often slower than the default code (e.g. unwrapped loops in numba, and in this example), so I would like to know if it doesn't speed up the code, because rather not use it then. A failure would save me benchmarking the slow-down.
celrod 5y ago
Did you start Julia with multiple threads? Julia unfortunately starts with only a single thread by default, and the slow time you reported (14 ms) makes it appear that this may have been the case. You can start `julia -t4` for 4 threads, for example. Or you can set the environmental variable `JULIA_NUM_THREADS=4`. If you were using the same number of threads for both, mind letting me know which CPU you're using?
- cycomanic 5y ago
  You were correct, I just tested with julia -t4 and that improved times by about a factor of 3. N = 2000 27.137 ms (24 allocations: 61.04 MiB) 49.496 ms (3 allocations: 61.05 MiB) 4.597 ms (2 allocations: 61.04 MiB) The CPU is a ryzen 3600 btw.
  - celrod 5y ago
    -t6 would give 6 threads. FWIW, the Ryzen 2600 with 6 threads took 6 ms, so the 3600's floating point and vector improvement is quite evident from your timing. (The 2600 emulates 256-bit operations with 2x 128-bit, while the 3600 has full support.)
    - cycomanic 5y ago
      Yes, I've also run with 6 threads and get an additional minor improvement on the loopvectorized code (3.9 vs 4.5ms). I have to say LoopVectorization.jl is quite an impressive piece of work. I have to try out some more julia at some point, unfortunately pretty much everything I do is using complex values.