I need to compute a matrix-matrix product, $A^T B$, where $A$ is $n \times r$ sparse, and $B$ is $n \times q$ dense. The number of rows $n$ is far larger than both $r$ and $q$. In fact $n$ is so large I cannot store the entire matrix $B$ in memory, so I build it one row at a time and update the matrix product $A^T B$ accordingly. The way I currently do this is to store $A$ in CSR format and then apply the following algorithm:
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], :))
This algorithm seems to have complexity $O(nkq)$, where $k$ is the number of non-zero entries in a row of $A$ (this is the same for all rows of $A$). This algorithm works fine, but is quite slow so I am looking for ideas on how to speed it up.
Dense BLAS algorithms often take advantage of blocks of data instead of operating on one row at a time, so my thinking was to do the same (instead of building B one row at a time, build a block of rows, where the block size is still small enough to fit in memory, and then apply a matrix multiplication algorithm). For this type of method, I don't think CSR format would work, since we would need the columns of $A$ in sequential memory locations, so CSC format might work better. A possible algorithm (with A in CSC format) is:
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, :))
I think the above algorithm has complexity $O(k r q \times NB)$, where $k$ is the number of non-zero entries in a column of $A$, and $NB$ is the number of blocks (num_blocks).
I haven't yet coded up this method, since I'm not sure it will actually be faster. Does anyone have any suggestions on how to make a fast algorithm to do this? Perhaps there is another sparse matrix storage format which would be better to use?