[SciPy-Dev] Improving performance of sparse matrix multiplication

[marc]

Dear Scipy developers,

We are developing a program that perform a large number of sparse matrix 
multiplications. We recently wrote a Python version of this program for 
several reasons (the original code is in Fortran).

We are trying now to improve the performance of the Python version and 
we noticed that one of the bottlenecks are the sparse matrix 
multiplications, as example,

import numpy as np
from scipy.sparse import csr_matrix

row = np.array([0, 0, 1, 2, 2, 2])
col = np.array([0, 2, 2, 0, 1, 2])
data = np.array([1, 2, 3, 4, 5, 6], dtype=np.float32)

csr = csr_matrix((data, (row, col)), shape=(3, 3))
print(csr.toarray())

A = np.array([1, 2, 3], dtype=np.float32)

print(csr*A)


I started to look at the Scipy code to see how this functions were 
implemented, and realized that there is no openmp parallelization over 
the for loops. Like in function  csr_matvec in sparse/sparsetools/csr.h 
(line 1120). Is it possible to parallelize this loops with openmp? Do 
you have maybe better ideas to improve the performances for this kind of 
operations?

Best regards,
Marc Barbry