# Optimized parallel routine for $X' W X$ with $W$ diagonal

$X$ is a dense matrix of real doubles, typically of size 20 million rows and 500 columns, and $W$ is a diagonal matrix of real, non-negative doubles stored as a vector. I'm working in C and have looked at the documentation for BLAS but couldn't find a routine that seemed like a perfect fit for this.

The machine has over 20 CPU cores and I would like to make use of them if possible. There is plenty of capacity for $X$ and $W$ to be fully in memory but I would like to avoid making copies of data unless really necessary as it is a shared computing resource.

I compute $X'WX$ hundreds of times. At each iteration $X$ is constant but $W$ changes.

I would be happy to write my own code for this in C with OpenMP but ideally would use an existing optimized routine if that would be fastest.

You haven't said anything about how much storage you have available to use for this computation. A 20 million by 500 array of double precision floats would require 80 gigabytes of RAM to store.

If you don't have that much RAM then you'll want to organize the computation to work with smaller strips of the array that you bring in from disk one at a time. For example, you could work with strips of 25 columns at a time and you'd need only 4 gigabytes of RAM to store each strip. You'd compute the results for each set of 25 columns and add the resulting 500 by 500 matrices together to get the final result.

If your weights $W$ are nonnegative, then you can simply take their square roots and write $W$ as $W=D^{2}$. Then your product is $(DX)^{T}(DX)$. You can scale by $D$ in one pass and then call the BLAS routine DGEMM to do the actual matrix multiplications.

If your weights $W$ are mixed positive and negative, then I would make a copy of $X$ and scale it by $W$ before performing the matrix multiplication by DGEMM. Again, if the matrices are too big to fit in memory you can block the computations by hand.

Libraries like Intel's MKL and AMD's ACML provide highly optimized implementations of the DGEMM routine. You should definitely take advantage of the efforts that have been put into optimizing DGEMM in such libraries.

• Thanks for the detailed answer Brian. I've edited my question to provide more detail. – user2179977 Nov 3 '15 at 13:04
• Brian, it's not as bad as you make out. Only 74.5 GB (80 billion bytes) of RAM would be needed to store the double precision array. :) – Mark L. Stone Nov 3 '15 at 14:26

Just an improvement on Brian's answer. Instead of using DGEMM, you can use DSYRK to multiply a matrix by its transpose.

• True, although you might find that DSYRK hasn't been optimized in your BLAS library as well as DGEMM. – Brian Borchers Nov 4 '15 at 5:24