I have to perform a [complex] basis transformation on a large number of [real] diagonal matrices: $$ \langle b_i | A | b_j \rangle = \sum_k \langle b_i | \bar{b}_k\rangle \langle\bar{b}_k | A | \bar{b}_k \rangle \langle\bar{b}_k | b_j \rangle $$
What is the most efficient way to perform this operation, ideally relying on BLAS/LAPACK or similar libraries?
It may be important that the $\bar{B}$ basis is much larger than the $B$ basis (that's why we're doing it, after all). Also, only the upper or lower triangular of the matrix $\langle b_i | A | b_j \rangle$ is actually needed, as its will end up in an eigenvalue solver. $\langle\bar{b}_k | b_j \rangle$ is stored with $ k $ as the fast index and $ j $ as the slow one. On distributed systems, the $k$ index is distributed very close to evenly.
I'll put my current implementation as an answer below. I don't think it is bad, but this accounts for a large enough fraction of run-time that any improvement would be awesome.