title: Multiplication of two huge dense matrices Hadamard-multiplied by a sparse matrix
body: <p>I have two dense matrices <code>A</code> and <code>B</code>, and each of them has a size fo <code>3e5x100</code>. Another sparse binary matrix, <code>C</code>, with size <code>3e5x3e5</code>. I want to find the following quantity: <code>C ∘ (AB')</code>, where <code>∘</code> is Hadamard product (i.e., element wise) and <code>B'</code> is the transpose of <code>B</code>. Explicitly calculating <code>AB'</code> will ask for crazy amount of memory (~500GB). Since the end result won't need the whole <code>AB'</code>, it is sufficient to only calculate the multiplication <code>A_iB_j'</code> where <code>C_ij != 0</code>, where <code>A_i</code> is the column <code>i</code> of matrix <code>A</code> and <code>C_ij</code> is the element at location <code>(i,j)</code> of the matrix <code>C</code>. A suggested approach would be like the algorithm below:</p>
  <pre class="lang-py prettyprint-override"><code>result = numpy.initalize_sparse_matrix(shape = C.shape)
  while True:
   (i,j) = C_ij.pop_nonzero_index() #prototype function returns the nonzero index and then points to the next nonzero index
   if (i,j) is empty:
     break
   result(i,j) = A_iB_j'
  </code></pre>
  <p>This algorithm however takes too much time. Is there anyway to improve it using <code>LAPACK/BLAS</code> algorithms? I  am coding in Python so I think <code>numpy</code> can be more human friendly wrapper for <code>LAPACK/BLAS</code>.</p>