计算科学 - 在许多对角矩阵上有效地改变基 - 吾爱随笔录

在许多对角矩阵上有效地改变基

计算科学线性代数表现布拉斯

2021-12-09 02:29:06

我必须对大量 [real] 对角矩阵执行 [complex] 基变换：

⟨ b_{i} | A | b_{j} ⟩ = \sum_{k} ⟨ b_{i} | {\bar{b}}_{k} ⟩ ⟨ {\bar{b}}_{k} | A | {\bar{b}}_{k} ⟩ ⟨ {\bar{b}}_{k} | b_{j} ⟩

$\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$

执行此操作的最有效方法是什么，理想情况下依赖 BLAS/LAPACK 或类似库？

基大得多可能很重要（毕竟这就是我们这样做的原因）。的上三角或下三角实际上是需要的，因为它最终会出现在特征值求解器中。存储作为快速索引，作为慢速索引。在分布式系统上，索引的分布非常接近均匀。 $\bar{B}$ $B$ $\langle b_i | A | b_j \rangle$ $\langle\bar{b}_k | b_j \rangle$ $k$ $j$ $k$

我将把我当前的实现作为答案放在下面。我不认为这很糟糕，但这占运行时间的很大一部分，任何改进都会很棒。

2个回答

不幸的是，所需的例程不是 BLAS 的一部分，但与 zherk 非常相似，后者执行 Hermitian rank-k 更新。特别是，zherk 支持的操作看起来像

A := α B B^{H} + β A,

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

其中只更新了的一个三角形（传统的变量名称在 zherk 原型中有所不同，其中通常是要更新的矩阵）。你想要的是概括 $A$ $C$

A := α B D B^{H} + β A,

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

其中是实对角矩阵。我需要类似的东西来计算 Hermitian 矩阵函数并编写了我自己的版本来计算 $D$

A := α B (D B)^{H} + β A,

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

用户可以指定，即使外积的左右矩阵不同，结果仍然是 Hermitian。您可以在此处查看我的分布式内存实现。我将泛化称为三角秩 K 更新 (Trrk)。它也可用于分解。 $LDL$

我通过的右侧，然后通过矩阵向量调用进行左侧。在一个压缩向量中累积，并在分布式系统上求和后，我重新展开到上三角。 $b_j$ $b_j$

  ! 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

我不满意：

BLAS 不知道转换的右侧和左侧是相同的（尽管缓存可能在具有缓存的机器上）。
逐元素乘法不是库调用。我注意到 BLAS 的 ZDOTC 明显优于 Fortran 的内在 dot_product（mkl，Intel 编译器），所以我对优化矢量化的那条线没有很大的信心。

其它你可能感兴趣的问题

上一篇算子的半有界性及其对稳定性的影响下一篇如何检测行驶道路上的关键转折点？