乘以对角矩阵的最佳方法是什么（在fortran中）

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

2021-12-09 14:20:08

什么是最好的计算方法：

Y = D X

$Y = D X$ 在哪里

D \in R^{m \times m}

$D \in \mathbb{R}^{m\times m}$ 是对角线并且

X \in C^{m \times n}

$X \in \mathbb{C}^{m \times n}$ 很一般。我最感兴趣的是这两种情况：

$m >> n$ , $m > 10^7$
$n >> m$ , $m < 10^4$

选项

我可以想到四种没有明显缺陷的方法：循环、forall、zgbmv上的循环、 zdscal上的循环。

环形

do i = 1,n
  do j = 1,m
    Y(j,i) = D(j) * X(j,i)
  enddo
enddo

优点：易于阅读，按顺序读取 D、X、Y
缺点：不重复使用 D

对所有人

forall (i = 1:n, j = 1:m) Y(j,i) = D(j) * X(j,i)

优点：简洁，给编译器自由？
缺点：给编译器自由？
注意：这些评论和这篇文章中对 forall 的先前讨论

zgbmv

 Dz = cmplx(D)
 do i = 1,n
   call zgbmv('N', m, m, 0, 0, one, D, 1, X(1,i), 1, zero, Y(1,i), 1)
 enddo

优点：类似于循环，但可能包含 BLAS 魔法
缺点：不重复使用 D，通过强制转换为复杂的 D 大小加倍

zdscal

 Y = X
 do i = 1,m
   call zdscal(n,D(i),Y(i,1),m)
 enddo

优点：重复使用 D，可能包含 BLAS 魔法
缺点：跨步读取 Y，如果不是就地需要复制

想法

两个主要的权衡似乎是顺序读取 X 与重用 D 和使用 fortran 库与将 D 视为真实而不是转换为复杂。在这两种情况下，自定义实现都可以两全其美，但我对特定于架构的参数持怀疑态度。最好的情况是本机表达操作（例如循环或forall）并告诉编译器完成其余操作。

1个回答

tl;dr使用循环

我的数字表明，ifort 足够聪明，可以识别循环、forall和do concurrent相同，并且在每种情况下都达到了我期望的“峰值”。另一方面，gfortran 用 and 做得不好（慢 10 倍或更多）forall，do concurrent尤其是在N变大时。ifort 和 gfortran 似乎对forall和产生相同的结果do concurrent。

我将 MKL 用于 BLASifort和gfortran. gbmv 比“峰值”慢 3 倍。对于小问题，'scal' 接近峰值，尤其是 small N，但很快就落后了。不过，它永远不会比gfortran's forall 和 do concurrent 更糟糕。

N在像我这样的系统（标准工作站配置）上，看起来循环是所有M. do concurrent没有优势forall：两者都不好。

您可以在此处找到结果表和代码。IBM 或 PGI 的结果是否相似？

有趣的是，作为复制时间倍数的性能Y = X并不取决于N. 我原以为重用D会提高性能N，类似于 GEMM 优于 GEMV。

笔记

$ ifort --version
ifort (IFORT) 14.0.1 20131008

$ gfortran --version
GNU Fortran (Debian 4.7.2-5) 4.7.2

其它你可能感兴趣的问题

上一篇对于非线性偏微分方程，源项是否应离散化为你j你j或平均超过(你j + 1+你j - 1) / 2(你j+1+你j-1)/2? 下一篇用 mkl LAPACK 求解矩阵的零空间