# Block matrix and DSYRK

I want to compute the matrix $$A = \sum_{i=1}^N v_i v_i^T$$ where each $$v_i$$ is a given vector of length $$2500$$, so that $$A$$ is $$2500 \times 2500$$, and my $$N$$ is about 2 million. Rather than call DSYR 2 million times, I found I can speed up the calculation significantly, by storing the $$v_i$$ in block matrices, and then calling DSYRK when the block matrix is full. I define a $$\textrm{block_size} \times 2500$$ matrix $$W_k$$ as: $$W_k = \pmatrix{ v_1^T \\ \vdots \\ v_b^T }$$ so that $$A = \sum_{k=1}^{\textrm{nblocks}} W_k^T W_k$$ The reason I set the rows of $$W_k$$ equal to the vectors $$v_i$$ and not the columns is because I use C and this way the incX of the rows will be 1 which I think helps to speed up the BLAS operations.

I am finding this procedure is still slower than I would like, and my question is how do I choose the $$\textrm{block_size}$$ (i.e. number of rows of $$W_k$$) to get the best performance?

Currently, I am choosing $$\textrm{block_size} = 50000$$, which I only chose so that the total size of the $$W_k$$ storage matrix is 1GB ($$50000 \times 2500 \times 8 = 1 GB$$). Is there a better choice of $$\textrm{block_size}$$ (number of rows of $$W_k$$) to get optimal performance?

In case it matters I am using the ATLAS CBLAS library on a machine with 24 cores.

• You’re going to have to access the matrix both by rows and columns no matter whether your start with your $W$ matrix or its transpose. It’s possible that your BLAS might be more efficient one way or the other. You could also try other BLAS implementations (MKL, OpenBLAS, etc.) Jan 11, 2019 at 19:38
• Of these ~2 million vectors $v_i$, at most 2500 can be linearly independent. Can't you use that somehow? Jan 14, 2019 at 7:24