为什么犰狳矩阵计算比Fortran慢得多

Posted

tags:

篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了为什么犰狳矩阵计算比Fortran慢得多相关的知识，希望对你有一定的参考价值。

我尝试通过Armadillo库使用矩阵实现从Fortran到C ++重写代码。两个代码的结果相同，但C ++代码比Fortran慢（> 10x）。代码涉及小矩阵（2x2,4x4）逆，乘法和加法。我在这里放了一部分相似的代码进行测试。

============================

clang++ cplusplus.cc -o cplusplus --std=c++14 -larmadillo -O2

ifort fort.f90 -o fort -O2

C ++代码时间：0.39404s

Fortran代码时间：0.068秒

============================

C ++代码：

#include <armadillo>
#include <iostream>

int main()
{
  const int niter = 1580000;
  const int ns = 3;
  arma::cx_cube m1(2, 2, ns), m2(2, 2, ns), m3(2, 2, ns);
  arma::wall_clock timer;
  timer.tic();
  for (auto i=0; i<niter; ++i) {
    for (auto j=0; j<ns; ++j)
      m1.slice(j) += m2.slice(j) * m3.slice(j);
  }
  double n = timer.toc();
  std::cout << "time: " << n << "s" << std::endl;
  return 0;
}

Fortran代码：

program main
  implicit none
  integer, parameter :: ns = 3, niter = 1580000
  complex*16 m1(2, 2, ns), m2(2, 2, ns), m3(2, 2, ns)
  integer i, j
  real :: start, finish
  call cpu_time(start)
  do i = 1, niter
     do j = 1, ns
        m1(1, 1, j) = m1(1, 1, j) + m2(1, 1, j) * m3(1, 1, j) + m2(1, 2, j) * m3(2, 1, j)
        m1(1, 2, j) = m1(1, 2, j) + m2(1, 1, j) * m3(1, 2, j) + m2(1, 2, j) * m3(2, 2, j)
        m1(2, 1, j) = m1(2, 1, j) + m2(2, 1, j) * m3(1, 1, j) + m2(2, 2, j) * m3(2, 1, j)
        m1(2, 2, j) = m1(2, 2, j) + m2(2, 1, j) * m3(1, 2, j) + m2(2, 2, j) * m3(2, 2, j)
     end do
  end do
  call cpu_time(finish)
  print *, "time: ", finish-start, " s"

end program main

====================================================================

关注@ewcz @ user5713492建议

============================

clang++ cplusplus.cc -o cplusplus --std=c++14 -larmadillo -O2

ifort fort.f90 -o fort -O2

ifort fort2.f90 -o fort2 -O2

C ++代码（cplusplus.cc）时间：0.39650s

Fortran代码（fort.f90）（显式操作）时间：0.020s

Fortran代码（fort2.f90）（matmul）时间：0.064s

============================

C ++代码（cplusplus.cc）：

#include <armadillo>
#include <iostream>
#include <complex>

int main()
{
  const int niter = 1580000;
  const int ns = 3;
  arma::cx_cube m1(2, 2, ns, arma::fill::ones),
    m2(2, 2, ns, arma::fill::ones),
    m3(2, 2, ns,arma::fill::ones);
  std::complex<double> result;
  arma::wall_clock timer;
  timer.tic();
  for (auto i=0; i<niter; ++i) {
    for (auto j=0; j<ns; ++j)
      m1.slice(j) += m2.slice(j) * m3.slice(j);
  }

  double n = timer.toc();
  std::cout << "time: " << n << "s" << std::endl;
  result = arma::accu(m1);
  std::cout << result << std::endl;
  return 0;
}

Fortran代码（fort.f90）：

program main
  implicit none
  integer, parameter :: ns = 3, niter = 1580000
  complex*16 m1(2, 2, ns), m2(2, 2, ns), m3(2, 2, ns)
  integer i, j
  complex*16 result
  real :: start, finish
  m1 = 1
  m2 = 1
  m3 = 1
  call cpu_time(start)
  do i = 1, niter
     do j = 1, ns
        m1(1, 1, j) = m1(1, 1, j) + m2(1, 1, j) * m3(1, 1, j) + m2(1, 2, j) * m3(2, 1, j)
        m1(1, 2, j) = m1(1, 2, j) + m2(1, 1, j) * m3(1, 2, j) + m2(1, 2, j) * m3(2, 2, j)
        m1(2, 1, j) = m1(2, 1, j) + m2(2, 1, j) * m3(1, 1, j) + m2(2, 2, j) * m3(2, 1, j)
        m1(2, 2, j) = m1(2, 2, j) + m2(2, 1, j) * m3(1, 2, j) + m2(2, 2, j) * m3(2, 2, j)
     end do
  end do
  call cpu_time(finish)
  result = sum(m1)
  print *, "time: ", finish-start, " s"
  print *, result

end program main

Fortran代码（fort2.f90）：

program main
  implicit none
  integer, parameter :: ns = 3, niter = 1580000
  complex*16 m1(2, 2, ns), m2(2, 2, ns), m3(2, 2, ns)
  integer i, j
  complex*16 result
  real :: start, finish
  m1 = 1
  m2 = 1
  m3 = 1
  call cpu_time(start)
  do i = 1, niter
     do j = 1, ns
        m1(:,:,j) = m1(:,:,j)+matmul(m2(:,:,j),m3(:,:,j))
     end do
  end do
  call cpu_time(finish)
  result = sum(m1)
  print *, "time: ", finish-start, " s"
  print *, result

end program main

======================================================================

复数可能是犰狳如此缓慢的原因之一。如果我在C ++中使用arma::cube而不是arma::cx_cube并在Fortran中使用real*8，那么时间是：

C ++代码时间：0.08s

Fortran代码（fort.f90）（显式操作）时间：0.012s

Fortran代码（fort2.f90）（matmul）时间：0.028s

但是，我的计算需要复数。奇怪的是，犰狳图书馆的计算时间增长非常大，但对于Fortran而言则略有增加。

答案

你没有在gfortran中计算任何东西。它可以在-O2级别看到您不使用m1的值，因此它完全跳过计算。同样在Fortran中，您的阵列未初始化，因此您可以使用NaN进行计算，这可能会大大减慢速度。

因此，您应该初始化数组并使用某种输入，如命令行，用户输入或文件内容，以便编译器无法预先计算结果。

然后您可以考虑将Fortran中的循环内容更改为

m1(:,:,j) = m1(:,:,j)+matmul(m2(:,:,j),m3(:,:,j))

这样才能与C ++的东西保持一致。 （gfortran在做这件事时似乎放慢了很多但是ifort对它非常满意。）

然后你必须在最后打印出你的数组，这样编译器就不会断定你正在计时的循环可以像gfortran那样被跳过。编辑修复程序，让我们了解新结果。

另一答案

我会说你的Fortran版本在这个特定的例子中从显式扩展到基本操作中获得了显着的利润。为了证明这一点，我们假设以下修改：

  implicit none
  integer, parameter :: ns = 3, niter = 1580000
  complex*16 m1(2, 2, ns), m2(2, 2, ns), m3(2, 2, ns)
  integer i, j
  real :: start, finish
  call cpu_time(start)
  m2 = 1
  m3 = 1
  do i = 1, niter
     do j = 1, ns
        !m1(1, 1, j) = m1(1, 1, j) + m2(1, 1, j) * m3(1, 1, j) + m2(1, 2, j) * m3(2, 1, j)
        !m1(1, 2, j) = m1(1, 2, j) + m2(1, 1, j) * m3(1, 2, j) + m2(1, 2, j) * m3(2, 2, j)
        !m1(2, 1, j) = m1(2, 1, j) + m2(2, 1, j) * m3(1, 1, j) + m2(2, 2, j) * m3(2, 1, j)
        !m1(2, 2, j) = m1(2, 2, j) + m2(2, 1, j) * m3(1, 2, j) + m2(2, 2, j) * m3(2, 2, j)
        m1(:, :, j) = m1(:, :, j) + MATMUL(m2(:, :, j), m3(:, :, j))
     end do
  end do
  WRITE(*, *) SUM(m1)
  call cpu_time(finish)
  print *, "time: ", finish-start, " s"

这里，最后，程序打印m1的总和，以便至少部分确定整个循环没有被消除。使用显式乘法（和-O2），我得到大约0.05s的运行时间，而一般MATMUL它大约是0.2s，即类似于犰狳方法......

此外，尽管Armadillo基于模板很多，因此通过slice()创建子多维数据集视图的许多函数调用可能会被消除，但是在使用Fortran时，你仍然原则上有一些开销，你直接操作连续的内存块。

以上是关于为什么犰狳矩阵计算比Fortran慢得多的主要内容，如果未能解决你的问题，请参考以下文章

为啥以下简单的并行化代码比 Python 中的简单循环慢得多？

为啥 map.values().stream() 比 Array.stream(array) 慢得多

pandas 比 numpy 慢得多？

为啥 gcc 的输出比 Visual Studio 慢得多（对于此代码）？

为啥 TensorFlow matmul() 比 NumPy multiply() 慢得多？

犰狳 inplace_plus 明显慢于“正常”加操作