sgemm does not multithread when dgemm does - Intel MKL

Question

I am using the ?GEMM functions from Intel MKL to multiply matrices. Consider the following two matrix multiplications:

            cblas_?gemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m,n,k,
                            1.0,
                            Matrix1,k,
                            Matrix2,n,
                            0.0,
                            A,n);

where m=1E5, and n=1E4, k=5. When I use pca_dgemm and pca_sgemm, this users all 12 cores, and executes beautifully.

However, when I do the following matrix multiplication:

    cblas_?gemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m,l,n,
                    1.0,
                    A,n,
                    Ran,l,
                    0.0,
                    Q,l);

where m=1E5, n=1E5, and l=7 (note the order of the parameters passed is differnet though. this is (m,n) * (n,l)). pca_dgemm uses all 12 cores and executes beautifully.

However, pca_sgemm does not. It uses only 1 core, and of course, takes much longer. Of course, for sgemm I am using arrays of floats, whereas for dgemm I am using an arrays of doubles.

Why could this be? They both give accurate results, but sgemm only multithreads on the former, whereas dgemm multithreads and both! How could simply changing the data type make this kind of difference?

Note that all arrays were allocated in using mkl_malloc using an alignment of 64.

Edit 2: Please also note that when l=12, in other words, with a larger matrix, it does thread in the sgemm. In other words, it is clearly that the sgemm version requires larger matrices to parallelize, but dgemm does not have this requirement. Why is this?

Answer 1

The MKL functions do quite a bit of work up-front to try to guess what is going to be the fastest way of executing an operation, so it's no surprise that it comes to a different decision when processing doubles or singles.

When deciding which strategy to take, it has to weigh the cost of doing the operation in a single thread against the overhead of launching threads to do the operation in parallel. One factor that will come into play is that SSE instructions can do operations on single-precision numbers twice as fast as double-precision numbers, so the heuristic might well decide that it's likely quicker to do the operation on singles as SSE SIMD operations on a single core rather than kicking of twelve threads to do it in parallel. Exactly how many it can do in parallel will depend on the details of your CPU architecture; SSE2, for instance, can do an operation on four single operands or two double-operands, while more recent SSE instruction sets support wider data.

I've found in the past that, for small matrices/vectors, it's often faster to roll your own functions than to use MKL. For instance, if all your operations are on 3-vectors and 3x3 matrices, it's quite a bit faster to just write your own BLAS functions in plain C and faster again to optimise them with SSE (if you can meet the alignment constraints). For a mix of 3- and 6-vectors, it's still faster to write your own optimised SSE version. This is because the cost of the MKL version deciding which strategy to use becomes a considerable overhead when the operations are small.