Scipy supporta il multithreading per la moltiplicazione con matrice sparsa quando si utilizza MKL BLAS?

Question

Secondo la documentazione MKL BLAS "Tutte le operazioni matrice matrice (livello 3) sono filettate per BLAS sia denso che sparse." http://software.intel.com/en-us/articles/parallelism-in-the-intel-math-kernel-library

Ho costruito Scipy con MKL BLAS. Usando il codice di prova qui sotto, vedo l'attesa velocità multithread per la moltiplicazione densa, ma non sparsa, della matrice. Ci sono delle modifiche a Scipy per abilitare le operazioni sparse multithread?

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

# test sparse matrix multiplication 
from scipy import sparse 
x = sparse.rand(10000,10000) 
t1 = time.time() 
foo = dot(x.T, x) 
print time.time() - t1 

Answer 1

Per quanto ne so, la risposta è no. Ma puoi costruire il tuo wrapper attorno alle routine multiple di MKL. Hai chiesto della moltiplicazione di due matrici sparse. Di seguito è riportato un codice wrapper che ho usato per moltiplicare una matrice sparsa per un vettore denso, quindi non dovrebbe essere difficile adattarlo (guarda il riferimento Intel MKL per mkl_cspblas_dcsrgemm). Inoltre, sii consapevole di come vengono memorizzati gli array scipy: il valore predefinito è coo, ma csr (o csc) potrebbe essere una scelta migliore. Ho scelto CSR, ma MKL supporta la maggior parte dei tipi (basta chiamare la routine appropriata).

Da quello che ho potuto vedere, sia il default di scipy che MKL sono multithreaded. Cambiando OMP_NUM_THREADS ho notato una differenza nelle prestazioni.

Per utilizzare la funzione seguente, se si dispone di una versione recente di MKL, è sufficiente assicurarsi di avere LD_LIBRARY_PATHS impostato per includere le relative directory MKL. Per le versioni precedenti, è necessario creare alcune librerie specifiche. Ho ricevuto le mie informazioni da IntelMKL in python

def SpMV_viaMKL(A, x): 
""" 
Wrapper to Intel's SpMV 
(Sparse Matrix-Vector multiply) 
For medium-sized matrices, this is 4x faster 
than scipy's default implementation 
Stephen Becker, April 24 2014 
stephen.beckr@gmail.com 
""" 

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", e.g., 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