Efficiently changing basis on many diagonal matrices

Question

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.

Answer 1

Unfortunately the desired routine is not part of BLAS but is very similar to zherk, which performs a Hermitian rank-k update. In particular, zherk supports operations which look like

$$A := \alpha B B^H + \beta A,$$

where only one triangle of $A$ is updated (the traditional variable names are different in the zherk prototype, where usually $C$ is the matrix to be updated). What you want is the generalization

$$A := \alpha B D B^H + \beta A,$$

where $D$ is a real diagonal matrix. I needed something similar for computing Hermitian matrix functions and wrote my own version which computes

$$A := \alpha B (D B)^H + \beta A,$$

where the user can specify that, even though the left and right matrices for the outer product are different, the result will still be Hermitian. You can see my distributed-memory implementation here. I call the generalization a Triangular Rank K update (Trrk). It is also useful for $LDL$ factorizations.

Answer 2

I do the right-side of the transformation $b_j$ by $b_j$ and then the left-side through a matrix-vector call. After accumulating in a condensed vector, and summing on distributed systems, I re-expand into the upper triangular.

  ! trans(k,j) is $<\bar{b}_k | b_j>$
  ij = 1
  do j = 1, N1
    ! right-side of transformation
    jtmp(1:N2) = A(1:N2) * conjg(trans(1:N2, j))
    ! left-size of transformation
    call ZGEMV('T', N2, j, one, &
                trans(:,:), N2, &
                jtmp, 1, zero, &
                res(ij), 1)
    ij = ij + j
  enddo

  ! global sum if needed
  [global sum(res)]

  ! expand to upper triangular
  ij=0
  do j = 1, N2; do i = 1,j
    ij = ij + 1
    A_new(i,j) = conjg(res(ij))
  enddo; enddo

I am unhappy with:

• BLAS not knowing that the right and left side of the transformations are the same (though the cache might, on machines with cache).

• The element-wise multiply not being a library call. I've noticed that BLAS's ZDOTC noticeably outperforms Fortran's intrinsic dot_product (mkl, Intel compilers), so I don't have a ton of confidence in that line being optimally vectorized.