# Computing sparse matrix products into a dense result

I need to assemble a matrix (in dense form, of moderate size, say dimension 1000) which is most easily expressed as the product of several (4) sparse matrices. These matrices are most easily expressed by their action (a matrix-vector multiplication routine), although I suppose it's not too hard to write them out explicitly.

My question is what is the best way for computing the product, keeping in mind that the result will wind up being stored as a dense matrix anyways. Clearly, expressing everything as a dense matrix and computing the products seems ridiculously inefficient, since the intermediates are still quite sparse.

More details: The actual form of the matrix is like $$A = DBD^HC$$ where $D$ is a mesh adjacency matrix (only two nonzeros per row), and $B$ and $C$ are banded matrices with bandwidth under 10.

If you really want all of the entries of the dense matrix product, you could compute it column by column, by computing the action on columns of the identity matrix.

The first question is how the assembled matrix $A$ will be used later, as the construction of the matrix itself is rarely a goal. In the question, it is clearly defined that one needs to have a dense matrix in the end of the day. However, for some purposes, having a dense matrix instead of its factorization into several matrices might be detrimental and should be avoided. The simplest example would be a solution of the system of linear equations using an iterative solver.
• As clipper pointed out, the entries of the dense matrix $A$ can be manually computed column-by-column by applying matrix-vector products to the columns of the identity matrix: $D(B(D^H(CI_n)))$, where $I_n$ is the $n$th column of the identity matrix. In such a way, you will be able to avoid forming dense intermediate results and take advantage of the sparsity of $D$, $B$, and $C$.
• Explicitly use some SparseBLAS implementation. For example, Intel MKL now supports subroutines mkl_?csrmultcsr and mkl_?csrmultd that compute the product of two sparse matrices and store the results in a sparse CSR or dense format, respectively.