$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.