# 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.) – Brian Borchers Jan 11 '19 at 19:38
• Of these ~2 million vectors $v_i$, at most 2500 can be linearly independent. Can't you use that somehow? – Christoph Jan 14 '19 at 7:24