add SLEEF in SLiM for SIMD vectorization of lots of math functions

[bhaller]

@andrewkern mentioned SLEEF in #578:

wrt further SIMDizing Eidos math functions-- transcendental functions aren't natively handled by SIMD instructions, but there is a quite portable library that most people use called sleef that we could consider bringing in: https://sleef.org/

If you were open to that direction, sleef could be included as a git submodule and we could build it as a static library with all its bells and whistles turned off. Alternatively, sleef has a way to build inline headers where we could include inline implementations really easily with a simple #include sleef.h

Aha, interesting. SLEEF looks quite remarkable and worth pursuing. I'm not sure which way of including it would be best. I would love the inline header option; but you have to build the inline header yourself, it is not provided. Do you know whether that has to be done on the target machine? I.e., does the generated header's code depends upon the target architecture/compiler? If so, that might make the benefits of that approach somewhat moot. I've never worked with git submodules, but I've heard a lot of people saying that it's painful (e.g., https://www.reddit.com/r/ExperiencedDevs/comments/1mttw8i/to_git_submodule_or_not_to/), and it would complicate the basic workflow of checking out and building SLiM, so I'm not crazy about that. But I'm not crazy about external link dependencies either. This is an interesting option: https://sleef.org/compile.xhtml#import but it would require a network connection to build, I guess (the first time, or every time?), and maybe we would not be able to control which version of SLEEF got imported and built against. So: a lot of options I'm not crazy about. If it is possible to build an inline header that is cross-platform and can just be included in SLiM's repo, that's clearly the best solution, to me, but I don't know if it's possible. The discussion of it here doesn't say: https://sleef.org/additional.xhtml#inline. It does note that SLEEF's "dispatcher" mechanism is not compatible with the inline header, which I think means that you'd have to have #ifdefs to select the correct version of a function to call, I think? Is that an acceptable sacrifice?

[bhaller]

The platforms SLEEF works on are limited: https://sleef.org/index.xhtml#environment. On macOS it has to be Clang (which is the standard compiler on macOS so that's OK); on Linux it has to be GCC (ditto); for Windows it has to be MSVC (which I think is not the compiler we target for SLiM on Windows at the moment). Everything else is "experimental", which I don't think would be acceptable for SLiM. So this is somewhat limiting.

[bhaller]

This is starting to feel pretty complicated and I'm wondering if it's worth it. It would speed up trig, power, and log functions (and related stuff like inverse trig, hyperbolic functions, etc.), and error and gamma functions. That stuff is present in Eidos, and doubtless some models do use it; but I don't think any of it is used in the SLiM core, or in a "typical" SLiM model. It also does min, max, abs, and mod, which are probably somewhat more commonly used, but still not really core functionality.

The one thing that it does that seems hella useful is exp(), which is core functionality because Gaussian functions are so commonly used in simulation. That would be great to have. But this is a lot of complication for that one thing. @andrewkern do you see another path here? Maybe we can see how SLEEF handles the exp() case specifically, and just implement a clone of that ourselves in SLiM to get the same thing for ourselves?

@petrelharp What are your thoughts here? Do you find yourself doing floating-point math like this on large vectors in SLiM, for some reason? Note that @andrewkern already sped up sqrt, abs, floor, ceil, trunc, round, sum, product in PR #578; this issue is about getting SIMD optimizations for more advanced math functions. @andrewkern What about you – do you actually use these kinds of functions on big vectors in your models?

It has always puzzled me that chip designers, even nowadays when chips are SO huge and complicated, don't include built-in instructions for (a) exp(), and (b) pseudo-random number generation. Doing those two things in-chip would make such a huge difference. Oh well.

[andrewkern]

The platforms SLEEF works on are limited: https://sleef.org/index.xhtml#environment. On macOS it has to be Clang (which is the standard compiler on macOS so that's OK); on Linux it has to be GCC (ditto); for Windows it has to be MSVC (which I think is not the compiler we target for SLiM on Windows at the moment). Everything else is "experimental", which I don't think would be acceptable for SLiM. So this is somewhat limiting.

yeah, annoying... what's the windows compiler we want?

[bhaller]

The platforms SLEEF works on are limited: https://sleef.org/index.xhtml#environment. On macOS it has to be Clang (which is the standard compiler on macOS so that's OK); on Linux it has to be GCC (ditto); for Windows it has to be MSVC (which I think is not the compiler we target for SLiM on Windows at the moment). Everything else is "experimental", which I don't think would be acceptable for SLiM. So this is somewhat limiting.

yeah, annoying... what's the windows compiler we want?

I believe at the moment we're set up for mingw64 and ucrt64. I'm not sure why; I don't know anything about Windows, the Windows support for SLiM was mostly set up by @rdinnager. But from what I've seen online, MSVC is a pretty weird compiler/toolchain that does a lot of things in a non-standard way, so I wouldn't be surprised if it took some work to get the SLiM build working with it.

[andrewkern]

okay a quick search reveals that MSYS2 already has these as packages:
https://packages.msys2.org/packages/mingw-w64-x86_64-sleef
https://packages.msys2.org/packages/mingw-w64-ucrt-x86_64-sleef

so i'm reasonably confident we can use pacman to grab these if building on windows with sleef

[bhaller]

okay a quick search reveals that MSYS2 already has these as packages: https://packages.msys2.org/packages/mingw-w64-x86_64-sleef https://packages.msys2.org/packages/mingw-w64-ucrt-x86_64-sleef

so i'm reasonably confident we can use pacman to grab these if building on windows with sleef

Yeah, it might work. The external linking path sounds like it requires LTO to be performant for SLEEF, which is presently not enabled because of issues users have encountered with it; but probably it is time to bite the bullet on that and say that people with broken toolchains need to figure out how to fix them, LTO is pretty basic. External link dependencies just always bring up weird issues, though. I don't know. What do you think about the bigger question: is this actually worth it? Is it going to make a real difference to performance for a significant fraction of SLiM users? I'm not sure I'm seeing it, except for exp(), which maybe we can solve in some other way. Who actually calls these functions, and why?

[andrewkern]

yeah good question. let me do some reading on what it would take to roll our own SIMD exp()-- it might be possible without too much pain

[andrewkern]

okay so i did a bunch of work last night and made good progress on incorporating SLEEF into SLiM. I have a draft PR up at #587 which now builds / tests cleanly across all the platforms targeted in CI.

This PR went the inline header route. Basically I built SLEEF headers for ARM (NEON) & X86_64 (AVX2) and include those headers directly in their own directory at eidos/sleef. This presents a portable way to build code off of SLEEF for use in SLiM.

The PR has exp, log, log10, and log2 implemented currently, but there are others we might want to include (e.g., pow, sin, cos, etc). These functions show good speedups on our AVX2 hardware:

1M elements, 10 repetitions (x86_64 AVX2):

Function	Without SLEEF	With SLEEF	Speedup
exp()	8.30 ms	4.05 ms	2.1x
log()	6.17 ms	3.37 ms	1.8x
log10()	10.79 ms	3.66 ms	2.9x
log2()	5.81 ms	3.99 ms	1.5x

For the moment I've left SSE4.2 out of SLEEF optimizations. It's doable, but would require another header file build and SSE4.2 machines are at a minimum 10 years old at this point-- I'm not sure we need to be working to get full optimizations out of those, but let me know if it's of interest.

I tried to document this all as well in markdown files that I've deposited in simd_benchmarks/*.md. There were a few patches to the header files that needed to be made, and the scripts to do that are included as well.

[andrewkern]

oh also note that I plopped this optimization into EidosValue_SP SpatialMap::ExecuteMethod_exp(...) to make it do something!

[bhaller]

Hi @andrewkern! Wow, you are amazing. :-> I had a quick look at the PR and it looks awesome. Thanks for documenting things thoroughly, I'm sure that will prove useful down the road.

I'll be happy to merge this when you think it's ready. Do you think these improvements could be applied to accelerating Eidos_ExecuteFunction_dnorm()? It's defined in eidos_functions_distributions.cpp. Accelerating that would make a real difference to a whole category of quantitative-trait SLiM models. Right now that calls gsl_ran_gaussian_pdf() but that function is very simple:

double
gsl_ran_gaussian_pdf (const double x, const double sigma)
{
  double u = x / fabs (sigma);
  double p = (1 / (sqrt (2 * M_PI) * fabs (sigma))) * exp (-u * u / 2);
  return p;
}

So if you could provide a SIMD version of that, using SLEEF's exp() and such, that would rock. :-> I think you could also speed up rnorm(), which calls gsl_ran_gaussian(). It'd be good to keep the original code in place for those, commented out of #ifdef'ed out.

I'm at the airport now, about to be offline for a day or two. Great stuff! :->

[andrewkern]

Glad you like this 😄

Yeah Eidos_ExecuteFunction_dnorm() seems like low hanging fruit for SIMD acceleration. That would be easy.

Any others I should look at? I was thinking the spatial interaction kernels might have SIMD opportunities, but haven't looked.

[bhaller]

PR #587 has been merged, but I'll keep this issue open since @andrewkern has further plans. :->

[andrewkern]

okay @bhaller --- #588 has the dnorm() addition.

In my benchmarking it leads to a ~2.7x speedup if there is a single mu/var, and ~1.5x if mu/var are vector valued.

This is a small PR and should be ready to merge once it passes tests.