I have a program, structured in two parts, $A$ and $B$. Both parts are capable of running as standalone units, and written in C++. $A$ is written for cluster systems, running entirely on CPU-nodes, connected via MPI (both on the same CPU for multiple cores and for different nodes), and $B$ is targeted to run on a single GPU. $A$ is doing FEM-calculations, while $B$ is only doing matrix-vector-multiplications, with a static dense matrix and varying vectors (matrix-size is typically 5kx5k-15kx15k complex double elements). The result of the matrix-vector-multiplication must be accessible from all threads in $A$.
In order to keep the memory load on the GPU as low as possible my strategy until now is to create the static matrix only in the first MPI thread in $A$, but the result matrix in all threads (giving me a large data chunk in GPU memory for the first thread, and small data chunks for all other threads) by creating as many instances of $B$ as there are threads in $A$, but only generating the matrix for one instance. After doing the multiplications in thread 0, the result is distributed using MPI to all the other threads. $A$ does calculations based on that result, and afterwards restarts $B$ with different parameters.
Now I have to port $B$ to the same nodes as $A$, without a GPU. I was looking at PBLAS-functions (PZGEMV) or the PLASMA-library, and at the Trilinos-library (which I already use in $A$, thus integration in $B$ should be easy), but there I am not sure if dense distributed matrix-vector multiplications are available. What would be the best strategy for porting, and to make the resulting program as efficient as possible?