JuliaLang/julia — pull request #5991 by timholy (2014-02-28 21:39:18+00:00)
Title: RFC: automated BLAS nthreads tuning

Lately there has been a certain amount of distress about how OpenBLAS uses threads. This PR is an attempt to do two things: (1) give Julia a mechanism to set its own policies, and (2) help OpenBLAS developers collect the performance data they need to make things better. CC @xianyi, @wernsaar.

The way this works is the following:

- Measure the performance of a function for different settings of `blas_set_num_threads` and different problem sizes. For an eigenvalue problem, there's just 1 parameter, `n`, as in the `nxn` matrix. For matrix multiply, there are 3 parameters, `m`, `k`, and `n` (for `mxk` times `kxn`).
- Store the performance data to disk, since it can take a while to acquire and is interesting in its own right
- Analyze the performance data, in particular with respect to a model
```
nthreads = c0 + c1*log(p1) + c2*log(p2) + ...
```
and optimizing the parameters to give the best overall performance, as defined by the average (over all problem sizes) of the ratio `time_for_chosen_number_of_threads/time_for_best_choice_of_threads`. In other words, a `2x2x2` task that takes `1e-5` seconds will contribute to the objective function just as much as a `512x1024x512` task that takes much longer.
- Enshrine the model parameters in the file `base/linalg/blas_tune.jl`, which provides a function to calculate the number of threads for each problem analyzed.

This PR enables just two functions, used by `eigfact!` and matrix multiplication. Perhaps we'd want to expand this to more functions. To add a new function, you have to:

1. Set up the test to [run](https://github.com/JuliaLang/julia/blob/teh/blastune/etc/blas_calibrate.jl#L59-L60) and be [analyzed](https://github.com/JuliaLang/julia/blob/teh/blastune/etc/blas_calibrate.jl#L59-L60)
2. To exploit the results, add calls to the auto-generated `nthreads` function inside the [relevant LAPACK function](https://github.com/JuliaLang/julia/blob/teh/blastune/etc/blas_calibrate.jl#L59-L60) (that `ccall` was already there, of course, and was not changed).

To run the performance testing on your own machine, navigate to the `etc` directory and do this:
```
include("blas_calibrate.jl")
perf_measure()
tune()
```
The timing results are stored in a JLD/HDF5 file. If you're mostly interested in the raw results, here's an example:
```
using HDF5, JLD, Winston
T, sz = jldopen("blas_calibrate.jld") do fid
g = fid["A_mul_B!"]
T = read(g, "T")
sz = read(g, "P")
close(g)
T, sz
end
rng(sz) = (log2(sz[1]), log2(sz[end]))
r = slice(T, 1, :,19,:)./slice(T,4,:,19,:)
p = imagesc(rng(sz[1]), rng(sz[3]), log10(r), (-1,1))
# Work around a bug in Winston: flip the y axis labels
yr = getattr(p, "yrange")
setattr(p, "yrange", (yr[2],yr[1])); display(p)
xlabel("log2(m)")
ylabel("log2(n)")
title("log10(T_{1thread}/T_{4threads})")
```

This has some significant limitations. Most notably, `blas_set_num_threads` does not control the actual number of threads used; instead, it seems to set an upper limit on the number of threads, but OpenBLAS might use fewer. That means we're not really taking advantage of the full range of possibilities---if OpenBLAS sometimes uses fewer threads than it should be, we can't measure that.

This fact helps explain some features of the performance data, which may be of interest to the OpenBLAS folks. For example, consider matrix multiplication:
![image](https://f.cloud.github.com/assets/1525481/2298427/fde9ec02-a0bc-11e3-932b-45df24a5f5a1.png)
Here, `x` and `y` axes are the size of `m` and `n`, and `k=1024`. You can see that the ratio is 1 when `m` and `n` are small, suggesting that OpenBLAS chooses to use 1 thread. But note two oddities:

- Just above the cutoff, it's obvious that using 4 threads is the better solution. Presumably OpenBLAS could do even better for such problems by using more threads. (The large value of `k` presumably accounts for the desirability of using multiple threads here.)
- Curiously, the performance advantage of threads tails off for larger `m` and `n` (see the upper right quadrant). This makes no sense to me (I thought each CPU had its own L1 cache).

Also on my machine, it's always best to use just 1 thread for `eigfact!`. This is consistent with some user reports, and seems to suggest a problem with `geevx!`.

Some things that could be done:
- If the OpenBLAS developers want to exploit this, consider adding an option to turn off all choices about the number of threads to use, and just let Julia set this number. That would give you the most useful performance data. It's also possible that the model-fit I've done here might be fairly flexible and robust mechanism for choosing the proper number of threads for different problem sizes internally within OpenBLAS.
- On the Julia side, we could run these tests on more functions
- Someone probably needs to think about how to handle and store tuning data for multiple architectures.