With the deprecation of sparsetools in 1.4, no supported method exists for taking the product of a sparse matrix and a dense array that accepts an output array (like the functions sparsetools.*_matvec). Also, the numpy.dot function produces an unexpected result.

import numpy as np
from scipy.sparse import csr_matrix
S = csr_matrix(np.eye(2))
d = np.arange(S.shape[1])
o = np.zeros(S.shape[0])
x = np.dot(S, d, out=o)  # returns len(d) scalar multiplied copies of S and ignores o.
y = np.dot(S, csr_matrix(d).T, out=o)  # returns correct csr_matrix, but ignores o.
z = S.dot(d)  # returns correct dense array, but has no out parameter.

This might be solvable by moving the spmatrix.__mul__ code to spmatrix.dot and adding an optional output array. Of course handling this optional parameter appropriately in all cases presents some challenges, but it might be acceptable to simply raise an exception in most cases.

Calls to various implementations of _mul_vector and _mul_multivector can easily be modified to accommodate the optional output array.