NOTE: this performance problem has been mostly resolved in MATLAB 7.6

BLAS performance bug and its effect on sparse "bench" in MATLAB 7.4

CHOLMOD uses mostly transposed BLAS calls (C=A*B' in DGEMM, C=A*A' in DSYRK, and X=B/L' in DTRSM). For DTRSM, the side='R' option is the culprit (for Intel and Sun, but not AMD). For DGEMM, the problem is made worse because the matrices tend to be tall and thin. These kernels are not intrinsically slow, but all BLAS authors have ignored their performance when creating multithreaded versions. With multithreading enabled, the performance of CHOLMOD, and the MATLAB sparse bench, will be slow until the BLAS authors fix the performance bug. In the meantime, if you do a lot of sparse matrix work in MATLAB, be sure to disable multithreading.

The following plots were provided by Pat Quillen, at The MathWorks, Inc. He compared the single-threaded and multi-threaded performance of DTRSM in the Intel MKL 9.0 library, with various inputs. CHOLMOD uses DTRSM applied to the right, lower triangular, transposed (X=B/L' in MATLAB notation, where L is lower triangular with non-unit diagonal).

TRSM performance results, x=L\b

In these plots, lower is better. The y-axis plots the time for each number of threads, divided by the time for one thread. Here, the four plots are for X=L\B (Lower, NoTrans), X=L'\B (Lower, Trans), X=U\B (Upper, NoTrans), and X=U'\B (Upper, Trans). The x-axis is the size of the matrix L or U (called "A" in the figures). These plots show that using 2 threads leads to nearly perfect speedup when the matrix A is large enough, with relative times close to 0.5. An ideal speedup of 4 would correspond a ratio of 0.25; these results show that DTRSM gets just under 0.3 with 4 threads. Thus, multi-threading is better, leading to speedups even for tiny matrices. Unfortunately, CHOLMOD does not use DTRSM this way. It applies the matrix to the right, not left.

TRSM performance results, x=b/L

In these plots, lower is still better; the x and y axes are the same as above. Here, the four plots are for X=B/L (Lower, NoTrans), X=B/L' (Lower, Trans), X=B/U (Upper, NoTrans), and X=B/U' (Upper, Trans). CHOLMOD uses X=B/L' (Lower, Trans), with many tiny matrices (5-by-5 and 10-by-10 are common). In sparse bench, CHOLMOD calls DTRSM about 8000 times, with many small matrices. These plots show that with 2 threads and a 5-by-5 matrix, DTRSM as used by CHOLMOD (Lower, Trans) is almost SIX times slower than when using just one thread. DTRSM is over NINE times slower with four threads.

Thus CHOLMOD is slow with multi-threading enabled

There are also performance problems in DGEMM. The AMD ACML BLAS does not have the same problem in DTRSM as the Intel MKL 9.0 BLAS does. However, it still has problems in DGEMM and perhaps DSRYK, when multi-threading is enabled. In the following plot, from the June 2007 Cleve's Corner, higher is better. The y-axis is the time for sparse bench with one thread, divided by the time with 1, 2, or 4 threads. The x-axis is the number of MATLAB "labs" running on a platform with 16 computers, each with two AMD processors, and each processor has 2 cores (64 cores total). Each "lab" is an independent instance of MATLAB, which can communicate via message passing with other labs. No message passing is done here; each lab is running its own independent "bench". Multi-threading is always worse than single-threading for this matrix, even when running just one lab. For larger matrices, multi-threading can help, but the matrix in sparse "bench" is not large enough (it's from a 300-by-300 2D mesh of the L-shaped membrane).

Thus, if you use sparse Cholesky in MATLAB 7.4, be sure to leave multi-threading turned off for best performance.