查看 scipy 稀疏产品的 Python 代码。请注意,它分两遍调用了编译后的代码。
看起来 mkl 代码做了同样的事情
https://software.intel.com/en-us/node/468640 https://software.intel.com/en-us/node/468640
如果 request=1,则例程仅计算长度为 m + 1 的数组 ic 的值,必须预先分配该数组的内存。退出时,值 ic(m+1) - 1 是数组 c 和 jc 中元素的实际数量。
如果 request=2,则先前已使用参数 request=1 调用了例程,则输出数组 jc 和 c 在调用程序中分配,并且它们的长度至少为 ic(m+1) - 1。
您首先分配ic
基于行数C
(你从输入中知道这一点),并使用以下命令调用 mkl 代码request=1
.
For request=2
你必须分配c
and jc
数组,基于大小ic(m+1) - 1
。这与数量不一样nnz
在输入数组中。
您正在使用request1 = c_int(0)
,这要求c
数组的大小是正确的——如果不实际执行的话你是不知道的dot
(or a request 1
).
=================
File: /usr/lib/python3/dist-packages/scipy/sparse/compressed.py
Definition: A._mul_sparse_matrix(self, other)
第 1 遍分配indptr
(注意大小),并传递指针(数据在此传递并不重要)
indptr = np.empty(major_axis + 1, dtype=idx_dtype)
fn = getattr(_sparsetools, self.format + '_matmat_pass1')
fn(M, N,
np.asarray(self.indptr, dtype=idx_dtype),
np.asarray(self.indices, dtype=idx_dtype),
np.asarray(other.indptr, dtype=idx_dtype),
np.asarray(other.indices, dtype=idx_dtype),
indptr)
nnz = indptr[-1]
第 2 遍分配indptr
(不同尺寸),并基于nnz
indices
and data
.
indptr = np.asarray(indptr, dtype=idx_dtype)
indices = np.empty(nnz, dtype=idx_dtype)
data = np.empty(nnz, dtype=upcast(self.dtype, other.dtype))
fn = getattr(_sparsetools, self.format + '_matmat_pass2')
fn(M, N, np.asarray(self.indptr, dtype=idx_dtype),
np.asarray(self.indices, dtype=idx_dtype),
self.data,
np.asarray(other.indptr, dtype=idx_dtype),
np.asarray(other.indices, dtype=idx_dtype),
other.data,
indptr, indices, data)
最后使用这些数组创建一个新数组。
return self.__class__((data,indices,indptr),shape=(M,N))
The mkl
库应该以同样的方式使用。
===================
https://github.com/scipy/scipy/blob/master/scipy/sparse/sparsetools/csr.h https://github.com/scipy/scipy/blob/master/scipy/sparse/sparsetools/csr.h
有c代码csr_matmat_pass1
and csr_matmat_pass2
===================
如果有帮助,这里是这些过程的纯 Python 实现。直译,不尝试利用任何数组操作。
def pass1(A, B):
nrow,ncol=A.shape
Aptr=A.indptr
Bptr=B.indptr
Cp=np.zeros(nrow+1,int)
mask=np.zeros(ncol,int)-1
nnz=0
for i in range(nrow):
row_nnz=0
for jj in range(Aptr[i],Aptr[i+1]):
j=A.indices[jj]
for kk in range(Bptr[j],Bptr[j+1]):
k=B.indices[kk]
if mask[k]!=i:
mask[k]=i
row_nnz += 1
nnz += row_nnz
Cp[i+1]=nnz
return Cp
def pass2(A,B,Cnnz):
nrow,ncol=A.shape
Ap,Aj,Ax=A.indptr,A.indices,A.data
Bp,Bj,Bx=B.indptr,B.indices,B.data
next = np.zeros(ncol,int)-1
sums = np.zeros(ncol,A.dtype)
Cp=np.zeros(nrow+1,int)
Cj=np.zeros(Cnnz,int)
Cx=np.zeros(Cnnz,A.dtype)
nnz = 0
for i in range(nrow):
head = -2
length = 0
for jj in range(Ap[i], Ap[i+1]):
j, v = Aj[jj], Ax[jj]
for kk in range(Bp[j], Bp[j+1]):
k = Bj[kk]
sums[k] += v*Bx[kk]
if (next[k]==-1):
next[k], head=head, k
length += 1
print(i,sums, next)
for _ in range(length):
if sums[head] !=0:
Cj[nnz], Cx[nnz] = head, sums[head]
nnz += 1
temp = head
head = next[head]
next[temp], sums[temp] = -1, 0
Cp[i+1] = nnz
return Cp, Cj, Cx
def pass0(A,B):
Cp = pass1(A,B)
nnz = Cp[-1]
Cp,Cj,Cx=pass2(A,B,nnz)
N,M=A.shape[0], B.shape[1]
C=sparse.csr_matrix((Cx, Cj, Cp), shape=(N,M))
return C
Both A
and B
不得不csr
。它使用转置,需要转换,例如B = A.T.tocsr()
.