Czy obsługa wielowątkowości scipy dla mnożenia macierzy rzadkich podczas korzystania z MKL BLAS?

Question

Zgodnie z dokumentacją MKL BLAS "Wszystkie operacje macierzy (poziom 3) są gwintowane zarówno dla gęstego, jak i rzadkiego BLASa." http://software.intel.com/en-us/articles/parallelism-in-the-intel-math-kernel-library

Zbudowałem Scipy z MKL BLAS. Używając poniższego kodu testowego, widzę oczekiwane przyspieszenie wielowątkowe dla gęstego, ale nie rzadkiego, mnożenia macierzy. Czy są jakieś zmiany w Scipy, aby włączyć wielowątkowe operacje rozrzedzone?

# 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

O ile mi wiadomo, odpowiedź brzmi nie. Ale możesz zbudować własne opakowanie wokół rzadkich procedur mnożenia MKL. Pytałeś o pomnożenie dwóch rzadkich matryc. Poniżej znajduje się pewien kod opakowania użyty do pomnożenia jednej rzadkiej macierzy razy gęstszy wektor, więc nie powinno być trudno przystosować się (spójrz na referencję Intel MKL dla mkl_cspblas_dcsrgemm). Należy również pamiętać o przechowywaniu tablic scipy: domyślnie jest to gruch, ale csr (lub csc) może być lepszym wyborem. Wybrałem csr, ale MKL obsługuje większość typów (po prostu wywołaj odpowiednią procedurę).

Z tego, co mogłem powiedzieć, zarówno domyślne scipy, jak i MKL są wielowątkowe. Zmieniając OMP_NUM_THREADS mogłem zauważyć różnicę w wydajności.

Aby skorzystać z poniższej funkcji, jeśli masz najnowszą wersję MKL, po prostu upewnij się, że masz LD_LIBRARY_PATHS zestaw, aby dołączyć odpowiednie katalogi MKL. W przypadku starszych wersji musisz zbudować określone biblioteki. Mam informacje od 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