Поддерживает ли scipy многопоточность для разреженного умножения матриц при использовании MKL BLAS?
Согласно документации МКЛ БЛАС »Все матрично-матричные операции (уровень 3) ориентированы как на плотную, так и на разреженную BLAS. " http://software.intel.com/en-us/articles/parallelism-in-the-intel-math-kernel-library
Я построил Сципи с MKL BLAS. Используя тестовый код ниже, я вижу ожидаемое многопоточное ускорение для плотного, но не разреженного, умножения матриц. Есть ли какие-либо изменения в Scipy для включения многопоточных разреженных операций?
# 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
Ответы на вопрос(1)
ответ - нет. Но вы можете создать свою собственную обертку вокруг подпрограмм разреженного умножения MKL. Вы спрашивали о умножении двух разреженных матриц. Ниже приведен код-обертка, который я использовал для умножения одной разреженной матрицы на плотный вектор, поэтому он не долженТяжело адаптироваться (см. справочник Intel MKL для mkl_cspblas_dcsrgemm). Кроме того, помните, как хранятся ваши scipy-массивы: по умолчанию это coo, но csr (или csc) может быть лучшим выбором. Я выбрал csr, но MKL поддерживает большинство типов (просто вызовите соответствующую процедуру).
Из того, что я мог сказать, обаs по умолчанию и MKL являются многопоточными. Путем измененияOMP_NUM_THREADS Я мог видеть разницу в производительности.
Чтобы использовать функцию ниже, если у вас есть последняя версия MKL, просто убедитесь, что у вас естьLD_LIBRARY_PATHS установить для включения соответствующих каталогов MKL. Для более старых версий вам нужно собрать несколько конкретных библиотек. Я получил информацию отIntelMKL в питоне
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
[email protected]
"""
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