From its benchmarks (https://chriselrod.github.io/LoopVectorization.jl/latest/exa...), a 9-line naive matrix multiplication routine in Julia + LV slightly edges out Intel's MKL.
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/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.
Julia uses avx instructions by default if your code is amenable to it.
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!
> What I did't see was how I could use the AVX2 instructions myself.
If you ever find yourself in a situation where you want manual control over vectorization, the package SIMD.jl [1] is pretty good for manual, handwritten vectorization. There's also VectorizationBase.jl [2] which LoopVectorization.jl uses. Which one of these two packages are most appropriate just kinda depends on what sort of interface you prefer.
However! The second shows an example with Base.llvmcall, which lets me write LLVM assembly. That should let me do anything I want.
However however - the stuff I actually do doesn't need all that, and could likely be implemented with LoopVectorization.
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-0... ).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.
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 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?