据我所知,答案是否定的。但是,您可以围绕MKL稀疏乘法例程构建自己的包装器。您询问了两个稀疏矩阵相乘的问题。以下是我用于将一个稀疏矩阵乘以一个密集向量的包装器代码,因此很容易进行适应(查看Intel MKL参考mkl_cspblas_dcsrgemm)。此外,请注意您的scipy数组存储方式:默认为coo,但csr(或csc)可能是更好的选择。我选择了csr,但MKL支持大多数类型(只需调用适当的例程)。
据我所知,scipy的默认设置和MKL都是多线程的。通过更改
OMP_NUM_THREADS
,我可以看到性能上的差异。
要使用下面的函数,如果您有最新版本的MKL,请确保已将
LD_LIBRARY_PATHS
设置为包括相关的MKL目录。对于旧版本,您需要构建一些特定的库。我从
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
if not sparse.isspmatrix_csr(A):
raise Exception("Matrix must be in csr format")
(m,n) = A.shape
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))
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))
if x.dtype.type is not np.double:
x = x.astype(np.double,copy=True)
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))
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
c_char("N")
更改为c_char(b'N'))
(一个字节字符串),但除此之外,这个代码非常好用。 - Nelewout