Listable built-in functions - Why are they so fast on numerical vectors?

Question

Consider the following small example, where I calculate Exp on a large vector by employing the listability which many mathematical built-in functions have:

data = RandomReal[1, 10^7];
RepeatedTiming[Exp[data];]
(* {0.0115, Null} *)

The list is packed but even unpacking it still has a good performance (it's most likely re-packed)

data2 = Developer`FromPackedArray[data];
RepeatedTiming[Exp[data2];]
(* {0.0705, Null} *)

As soon as we append a non-numerical element or a number of a different type, the performance drops drastically

data3 = Append[data, 1];
RepeatedTiming[Exp[data3];]
(* {3.11, Null} *)

Now, I'm wondering if I missed something in the documentation or if there is a place, where it is explained how exactly these calls are made so performant. Specifically, I want to know

1. Does such a call use vectorization, is it parallelized or does it just call fast, but sequential code?
2. How does this compare to Compile with Parallelization turned on?

Answer 1

I'm collecting what we have found out and I'm concentrating for the moment on the call Exp[data] where data is a vector of reals in machine precision. For this, the following two statements seem to apply in Mathematica 11.3 with a modern CPU:

1. Single instruction, multiple data (SIMD) instructions are used provided by the MMX technology and Streaming SIMD Extensions
2. Mathematica employs the Intel Integrated Performance Primitives library under the hood. For the Exp case, the function ippsExp_64fc_A53 is called.

Using a slightly more complex example, like

Exp[x]*Log[x + 1.0]

shows that all parts of the expression are relaid to the IPP libraries

• libippvm.so (ippsLn_64f_A53)
• libippvm.so (ippsExp_64f_A53)
• libipps.so (ippsAddC_64f)
• libippvm.so (ippsMul_64f_A53)

Evidently, not all functions are available in the IPP and then it gets complicated and there is no general rule. For instance

data = RandomReal[1, 10^7];
Do[LogGamma[data], {200}] // AbsoluteTiming
Do[Log[Gamma[data]], {200}] // AbsoluteTiming
(* {10.2816, Null} *)
(* {9.09973, Null} *)

It appears that the Gamma function is implemented in libWolframEngine and LogGamma is a combination of Wolfram code for Gamma and the log function from libm on my machine. The second call however uses IPP for the logarithm and appears to be a bit faster here.

As soon as we set a specific precision, different implementations are used that partly employ the The GNU Multiple Precision Arithmetic Library. As expected, the runtime drops significantly several orders of magnitude (100x - 1000x).

My second question was, how this relates to compiled code. I want to differentiate two cases

1. A parallelized compiled function that gets one number and is made listable.
2. A compiled function that receives a tensor of numbers and calls Exp[list] in its body

In the first case, it appears that no vectorization is used at all. The parallelism comes from spawning multiple threads. The second case is more interesting. My tests showed that even in the compiled code, IPP library functions are called. For simple cases, however, the runtime of the compiled code cannot compete with simply using high-level calls.

Answer 2

Here is one that I could find in Wolfram Support:

"The kernel uses highly optimized multithreaded libraries such as Intel MKL and IPP, which are tuned for optimal performance and take advantage of advanced CPU features when available. This is important for vectorized machine arithmetic and numerical linear algebra (BLAS, …) routines, which are fundamental building blocks for many computational tasks."

I wonder if further details can be found in the Wolfram documentation.