我需要计算一个矩阵-矩阵乘积，其中是稀疏的，是密集的。行数远大于和。实际上太大了，我无法将整个矩阵存储在内存中，所以我一次构建一行并相应地更新矩阵乘积。我目前这样做的方式是将存储为 CSR 格式，然后应用以下算法：$A^T B$$A$$n \times r$$B$$n \times q$$n$$r$$q$$n$$B$$A^T B$$A$

loop i = 1, ..., n     // loop over rows of B
  bi = build_row(B, i) // construct row i of B
  loop p = A->row_ptr[i], ..., A->row_ptr[i+1] - 1 // loop over non-zero entries in row i of A
    DAXPY(A->data[p], bi, ATB(A->col_ind[p], :))

该算法似乎具有复杂性，其中的一行中非零条目的数量（这对于 A 的所有行都是相同的。该算法运行良好，但速度很慢，因此我正在寻找有关如何加快速度的想法。$O(nkq)$$k$$A$$A$

密集的 BLAS 算法经常利用数据块而不是一次操作一行，所以我的想法是做同样的事情（而不是一次构建 B 一行，而是构建一个行块，其中块大小为仍然足够小以适合内存，然后应用矩阵乘法算法）。对于这种类型的方法，我认为 CSR 格式行不通，因为我们需要的列在顺序内存位置，因此 CSC 格式可能会更好。一种可能的算法（使用 CSC 格式的 A）是：$A$

loop i = 1, ..., num_blocks
  row0 = (i-1) * block_size           // row of start of block in B
  row1 = row0 + block_size            // row of end of block in B
  blockB = build_block(B, row0, row1) // build B(row0:row1,:)
  loop j = 1, ..., r // loop over rows of A^T
    loop p = A->col_ptr[j], ..., A->col_ptr[j+1]-1 // loop over non-zero entries in row j of A^T
      if (row0 <= A->row_ind[p]) AND (A->row_ind[p] <= row1) // if this row of A is in the right range, perform DAXPY operation
        DAXPY(A->data[p], blockB(A->row_ind[p]-row0, :), ATB(j, :))

我认为上述算法具有复杂度，其中列中非零条目的数量，而是块数（num_blocks）。$O(k r q \times NB)$$k$$A$$NB$

我还没有编写这种方法，因为我不确定它实际上会更快。有没有人对如何制作一个快速算法来做到这一点有任何建议？也许还有另一种更好用的稀疏矩阵存储格式？