6 ms·

At a conference presentation the other day a speaker made the claim that well-written Julia code was as fast as well-written C code while as expressive as Pytho

by dalke 6y ago

At a conference presentation the other day a speaker made the claim that well-written Julia code was as fast as well-written C code while as expressive as Python, so it solved the two language problem. My Python package uses a C extension for performance. That C extension uses AVX2 intrinsics. How well does Julia support intrinsics? Even supporting something like __builtin_popcountll() would be nice, but http://rosettacode.org/wiki/Population_count#Julia http://rosettacode.org/wiki/Population_count#Julia suggests that it's not supported in Julia. I've no experience with Julia and my attempt at finding this out on my own failed - does anyone here know about Julia's support for intrinsics?

bacr 6y ago
My understanding is that the Julia community is quite interested in having SIMD via e.g. AVX “just work”. I recall reading this post on it a while back: https://juliacomputing.com/blog/2017/09/27/auto-vectorization-in-julia.html https://juliacomputing.com/blog/2017/09/27/auto-vectorizatio...
- dalke 6y ago
  Sure, but you can only do that if either that's a way to express what you want directly or, in the example you gave, there's a common idiomatic style that the compiler can recognize and handle. What is the idiomatic way to write the popcount of the intersection of two 256-byte byte strings? My C code is: static int byte_intersect_256(const unsigned char *fp1, const unsigned char *fp2) { int num_words = 2048 / 64; int intersect_popcount = 0; /* Interpret as 64-bit integers and assume possible mis-alignment is okay. */ uint64_t *fp1_64 = (uint64_t *) fp1, *fp2_64 = (uint64_t *) fp2; for (int i=0; i<num_words; i++) { intersect_popcount += __builtin_popcountll(fp1_64[i] & fp2_64[i]); } return intersect_popcount; } I haven't figured out the Julia way to write it so it would use the POPCNT instruction (if available), the AVX2 popcount technique (if available), or the VPOPCNTDQ AVX-512 instruction (if available) - falling back, I suppose, to the SSSE3 and Lauradoux implementations - the last being the fastest generic C implementation I found. (See https://jcheminf.biomedcentral.com/articles/10.1186/s13321-019-0398-8/tables/1 https://jcheminf.biomedcentral.com/articles/10.1186/s13321-0... ).
timholy 6y ago
Check out LoopVectorization: https://github.com/chriselrod/LoopVectorization.jl https://github.com/chriselrod/LoopVectorization.jl From its benchmarks (https://chriselrod.github.io/LoopVectorization.jl/latest/examples/matrix_multiplication/#Matrix-Multiplication-1 https://chriselrod.github.io/LoopVectorization.jl/latest/exa...), a 9-line naive matrix multiplication routine in Julia + LV slightly edges out Intel's MKL.
- dalke 6y ago
  It looks like that @avx macro just tells Julia that it's okay to use the AVX instructions? My specific question is, how do I tell Julia that I want to compute the popcount of the intersection of two byte strings of length 256 bytes? A reference C code is at http://www.dalkescientific.com/writings/diary/archive/2020/10/07/intersection_popcount.html http://www.dalkescientific.com/writings/diary/archive/2020/1... in byte_intersect_256() and threshold_bin_tanimoto_search(); my blog posts shows the important parts - I link to the full definition for the C code.
  - eigenspace 6y ago
    @avx does a ton more than just use AVX instructions. It'll reorder and unroll loops when advantageous, swap out some functions for more vectorizable version of those functions and a few other tricks. Julia uses avx instructions by default if your code is amenable to it.
    - dalke 6y ago
      Sorry, I see how what I wrote could be interpreted that way, and part of what I wrote was out of ignorance. While I didn't write it, I assumed the macro was doing some sort of equivalent to the template metaprogramming I've heard about in C++ to do similar things. What I did't see was how I could use the AVX2 instructions myself. Checking now, since Julia's count_ones() maps to the LLVM popcount instruction, and recent clang versions know how to optimize that fixed-length sequence in C even for AVX-512, the Julia equivalent to the code I wrote should have good performance. There are a few optimizations (keeping one AVX register loaded with a constant byte string, and using prefetch instructions) which might be missing. I'll be talking with the conference participant who brought up Julia to work this out in more detail. Thanks for the comment!
Oreb 6y ago
Isn’t __builtin_pocountll() exactly the same as Julia’s count_ones() function?
- dalke 6y ago
  Thanks! Yes! I kept looking for "popcount" and "population count", not "Number of ones in the binary representation of x." A Julia person should update the Rosetta Code implementation to use that function. That function maps to ctpop_int which is implemented as uu_iintrinsic_slow(LLVMCountPopulation, ctpop_int, u). I ... don't know what the "slow" means. There are follow-on problems, like my AVX2 implementation for 128 bytes stores one of the byte strings in AVX registers, and I get extra performance by prefetching. Once I get to the AVX-512 implementation, the same will happen. I don't know how well the optimizer can handle those details. ("My" = implemented by HN user nkurz.) But this is enough of a clue for a first decently-high-performance implementation - thanks again! I talked with the person from the conference who promoted Julia. He's going to work on this/we'll work on it. I'm looking forward to seeing how the timings work out.

by dalke 6y ago

bacr 6y ago

My understanding is that the Julia community is quite interested in having SIMD via e.g. AVX “just work”. I recall reading this post on it a while back: https://juliacomputing.com/blog/2017/09/27/auto-vectorization-in-julia.html https://juliacomputing.com/blog/2017/09/27/auto-vectorizatio...

dalke 6y ago
Sure, but you can only do that if either that's a way to express what you want directly or, in the example you gave, there's a common idiomatic style that the compiler can recognize and handle. What is the idiomatic way to write the popcount of the intersection of two 256-byte byte strings? My C code is: static int byte_intersect_256(const unsigned char *fp1, const unsigned char *fp2) { int num_words = 2048 / 64; int intersect_popcount = 0; /* Interpret as 64-bit integers and assume possible mis-alignment is okay. */ uint64_t *fp1_64 = (uint64_t *) fp1, *fp2_64 = (uint64_t *) fp2; for (int i=0; i<num_words; i++) { intersect_popcount += __builtin_popcountll(fp1_64[i] & fp2_64[i]); } return intersect_popcount; } I haven't figured out the Julia way to write it so it would use the POPCNT instruction (if available), the AVX2 popcount technique (if available), or the VPOPCNTDQ AVX-512 instruction (if available) - falling back, I suppose, to the SSSE3 and Lauradoux implementations - the last being the fastest generic C implementation I found. (See https://jcheminf.biomedcentral.com/articles/10.1186/s13321-019-0398-8/tables/1 https://jcheminf.biomedcentral.com/articles/10.1186/s13321-0... ).

timholy 6y ago

Check out LoopVectorization: https://github.com/chriselrod/LoopVectorization.jl https://github.com/chriselrod/LoopVectorization.jl From its benchmarks (https://chriselrod.github.io/LoopVectorization.jl/latest/examples/matrix_multiplication/#Matrix-Multiplication-1 https://chriselrod.github.io/LoopVectorization.jl/latest/exa...), a 9-line naive matrix multiplication routine in Julia + LV slightly edges out Intel's MKL.

dalke 6y ago
It looks like that @avx macro just tells Julia that it's okay to use the AVX instructions? My specific question is, how do I tell Julia that I want to compute the popcount of the intersection of two byte strings of length 256 bytes? A reference C code is at http://www.dalkescientific.com/writings/diary/archive/2020/10/07/intersection_popcount.html http://www.dalkescientific.com/writings/diary/archive/2020/1... in byte_intersect_256() and threshold_bin_tanimoto_search(); my blog posts shows the important parts - I link to the full definition for the C code.
- eigenspace 6y ago
  @avx does a ton more than just use AVX instructions. It'll reorder and unroll loops when advantageous, swap out some functions for more vectorizable version of those functions and a few other tricks. Julia uses avx instructions by default if your code is amenable to it.
  - dalke 6y ago
    Sorry, I see how what I wrote could be interpreted that way, and part of what I wrote was out of ignorance. While I didn't write it, I assumed the macro was doing some sort of equivalent to the template metaprogramming I've heard about in C++ to do similar things. What I did't see was how I could use the AVX2 instructions myself. Checking now, since Julia's count_ones() maps to the LLVM popcount instruction, and recent clang versions know how to optimize that fixed-length sequence in C even for AVX-512, the Julia equivalent to the code I wrote should have good performance. There are a few optimizations (keeping one AVX register loaded with a constant byte string, and using prefetch instructions) which might be missing. I'll be talking with the conference participant who brought up Julia to work this out in more detail. Thanks for the comment!

Oreb 6y ago

Isn’t __builtin_pocountll() exactly the same as Julia’s count_ones() function?

dalke 6y ago
Thanks! Yes! I kept looking for "popcount" and "population count", not "Number of ones in the binary representation of x." A Julia person should update the Rosetta Code implementation to use that function. That function maps to ctpop_int which is implemented as uu_iintrinsic_slow(LLVMCountPopulation, ctpop_int, u). I ... don't know what the "slow" means. There are follow-on problems, like my AVX2 implementation for 128 bytes stores one of the byte strings in AVX registers, and I get extra performance by prefetching. Once I get to the AVX-512 implementation, the same will happen. I don't know how well the optimizer can handle those details. ("My" = implemented by HN user nkurz.) But this is enough of a clue for a first decently-high-performance implementation - thanks again! I talked with the person from the conference who promoted Julia. He's going to work on this/we'll work on it. I'm looking forward to seeing how the timings work out.