Does scipy support multithreading for sparse matrix multiplication when using MKL BLAS?

Question

According to MKL BLAS documentation "All matrix matrix operations (level 3) have threads for tight and sparse BLAS." http://software.intel.com/en-us/articles/parallelism-in-the-intel-math-kernel-library

I created Scipy with MKL BLAS. Using the test code below, I see the expected multi-threaded acceleration for dense but not sparse matrix multiplication. Are there any changes to Scipy to enable multi-threaded sparse operations?

# test dense matrix multiplication from numpy import * import time x = random.random((10000,10000)) t1 = time.time() foo = dot(xT, x) print time.time() - t1 # test sparse matrix multiplication from scipy import sparse x = sparse.rand(10000,10000) t1 = time.time() foo = dot(xT, x) print time.time() - t1 

Answer 1

As far as I know, the answer is no. But you can create your own wrapper around sparse multiple MKL routines. You asked about the multiplication of two sparse matrices. The following is the shell code that I used to multiply one sparse matrix by a different dense vector, so it’s hard to adapt (see Intel MKL Link for mkl_cspblas_dcsrgemm). Also, remember how your scipy arrays are stored: coo is used by default, but csr (or csc) might be a better choice. I chose csr, but MKL supports most types (just call the appropriate procedure).

From what I could say, both scipy default and MKL are multi-threaded. By changing OMP_NUM_THREADS , I saw a performance difference.

To use the function below, if you have a recent version of MKL, just make sure you have LD_LIBRARY_PATHS to include the appropriate MKL directories. For older versions, you need to build some specific libraries. I got info from IntelMKL in python

 def SpMV_viaMKL( A, x ): """ Wrapper to Intel SpMV (Sparse Matrix-Vector multiply) For medium-sized matrices, this is 4x faster than scipy default implementation Stephen Becker, April 24 2014 [e-mail] """ import numpy as np import scipy.sparse as sparse from ctypes import POINTER,c_void_p,c_int,c_char,c_double,byref,cdll mkl = cdll.LoadLibrary("libmkl_rt.so") SpMV = mkl.mkl_cspblas_dcsrgemv # Dissecting the "cspblas_dcsrgemv" name: # "c" - for "c-blas" like interface (as opposed to fortran) # Also means expects sparse arrays to use 0-based indexing, which python does # "sp" for sparse # "d" for double-precision # "csr" for compressed row format # "ge" for "general", eg, the matrix has no special structure such as symmetry # "mv" for "matrix-vector" multiply if not sparse.isspmatrix_csr(A): raise Exception("Matrix must be in csr format") (m,n) = A.shape # The data of the matrix data = A.data.ctypes.data_as(POINTER(c_double)) indptr = A.indptr.ctypes.data_as(POINTER(c_int)) indices = A.indices.ctypes.data_as(POINTER(c_int)) # Allocate output, using same conventions as input nVectors = 1 if x.ndim is 1: y = np.empty(m,dtype=np.double,order='F') if x.size != n: raise Exception("x must have n entries. x.size is %d, n is %d" % (x.size,n)) elif x.shape[1] is 1: y = np.empty((m,1),dtype=np.double,order='F') if x.shape[0] != n: raise Exception("x must have n entries. x.size is %d, n is %d" % (x.size,n)) else: nVectors = x.shape[1] y = np.empty((m,nVectors),dtype=np.double,order='F') if x.shape[0] != n: raise Exception("x must have n entries. x.size is %d, n is %d" % (x.size,n)) # Check input if x.dtype.type is not np.double: x = x.astype(np.double,copy=True) # Put it in column-major order, otherwise for nVectors > 1 this FAILS completely if x.flags['F_CONTIGUOUS'] is not True: x = x.copy(order='F') if nVectors == 1: np_x = x.ctypes.data_as(POINTER(c_double)) np_y = y.ctypes.data_as(POINTER(c_double)) # now call MKL. This returns the answer in np_y, which links to y SpMV(byref(c_char("N")), byref(c_int(m)),data ,indptr, indices, np_x, np_y ) else: for columns in xrange(nVectors): xx = x[:,columns] yy = y[:,columns] np_x = xx.ctypes.data_as(POINTER(c_double)) np_y = yy.ctypes.data_as(POINTER(c_double)) SpMV(byref(c_char("N")), byref(c_int(m)),data,indptr, indices, np_x, np_y ) return y