3
$\begingroup$

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.

|cite|improve this question
$\endgroup$
  • $\begingroup$ Since you're using BLAS/LAPACK you may find this useful. A surprisingly large fraction of time can be spent in mundane ways you wouldn't have guessed, that can be easily fixed. $\endgroup$ – Mike Dunlavey Sep 1 '12 at 14:46
5
$\begingroup$

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.

|cite|improve this answer
$\endgroup$
  • $\begingroup$ I've found that this solution works quite well. If $D$ is positive, I scale $B$ by $\sqrt(D)$ and then call zherk. If $D$ has negative elements or is complex, I scale $B^H$ by $D$ and then just call ZGEMM, blocking along the diagonal if $B$ is large enough. I still wish this was just in BLAS... $\endgroup$ – Max Hutchinson Jan 31 '13 at 6:04
0
$\begingroup$

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.

|cite|improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.