Q:

Parallel computation of big covariance matrices

We need to compute covariance matrices with sizes ranging from $10000\times10000$ to $100000\times100000$. We have access to GPUs and clusters, we wonder what is the best parallel approach for speeding up these computations.

A:

The first thing is to recognize that you can do this using BLAS.  If you data matrix is $X = [x_1 x_2 x_3 ...] \in \mathbb{R}^{m\times n}$ (each $x$ is a column vector corresponding to one measurement; rows are trials), then you can write the covariance as:
$$ C_{ij} = E[x_i,x_j] - E[x_i] E[x_j] = \frac{1}{n} \sum_k x_{ik} x_{jk} - \frac{1}{n^2} \left(\sum_k x_{ik} \right) \left(\sum_k x_{jk}\right) $$
We can write this as:
$$C = \frac{1}{n} X^T X - \frac{1}{n^2} (1^T X)^T (1^T X) $$
where $(1^T)$ is the row-vector with all elements 1 so $(1^T X)$ is a row vector of the column sums of $X$. This can be written completely as BLAS, where $X^T X$ is either a GEMM or, better yet, a SYRK/HERK and you can get the $(1^T X) = b$ with a GEMV, $b^T b$ with GEMM or SYRK/HERK again, and the prefactors with SCAL.
Your data and result matrices can be around 64GB, so you're not going to fit on a single node or a node's worth of GPUs.  For a non-GPU cluster, you might want to look at PBLAS, which feels like scalapack.  For GPUs, multi-node libraries aren't quite there yet.  Magma has some sort of underlying parallel BLAS implementation, but it may not be user friendly.  I don't think CULA does multi-node yet, but its something to keep an eye on.  CUBLAS is single-node.
I'd also suggest that you strongly consider implementing the parallelism yourself , especially if you are familiar with MPI and have to hook this into an existing code-base.  That way, you can switch between CPU and GPU BLAS easily and begin and end with the data exactly where you want it.  You shouldn't need more than a few MPI_ALLREDUCE calls.

