logo标签列表

使用显式循环进行矩阵乘法，Eigen比Fortran慢得多。

c++fortraneigen
3

我尝试使用Eigen库将一个2000x2000的矩阵乘法从Fortran重写为C++，但是我发现Eigen中的for循环比Fortran中的do循环要慢得多（> 3倍）。以下是代码：

test.f90

program main
implicit none
integer :: n,i,j,k
integer :: tic,toc
real(8),ALLOCATABLE ::a(:,:),b(:,:),c(:,:)
real(8) :: s

n = 2000
allocate(a(n,n),b(n,n),c(n,n))
do i=1,n
    do j =1,n
        a(j,i) = i * 1.0
        b(j,i) = i * 1.0
    enddo
enddo

call system_clock(tic)
do j=1,n
    do i=1,n
        s = 0.0
        do k=1,n
            s = s + a(i,k) * b(k,j)
        enddo
        c(i,j) = s
    enddo
enddo
call system_clock(toc)
print*,'Fortran with loop:', (toc - tic) / 1000.0

call system_clock(tic)
c = matmul(a,b)
call system_clock(toc)
print*,'Fortran with matmul:', (toc - tic) / 1000.0


DEALLOCATE(a,b,c)
end

test.cpp

#include<Eigen/Core>
#include<time.h>
#include<iostream>
using Eigen::MatrixXd;

int main(){
    int n = 2000;
    MatrixXd a(n,n),b(n,n),c(n,n);
    for(int i=0;i<n;i++){
    for(int j=0;j<n;j++){
            a(i,j) = i * 1.0;
            b(i,j) = j * 1.0;
        }
    }
    clock_t tic,toc;
    tic = clock();
    for(int j=0;j<n;j++){
        for(int i=0;i<n;i++){
            double s= 0.0;
            for(int k=0;k<n;k++){
                s += a(i,k) * b(k,j);
            }
            c(i,j) = s;
        }
    }
    toc = clock();
    std::cout << (double)((toc - tic)) / CLOCKS_PER_SEC << std::endl;

    tic = clock();
    c=  a * b;
    toc = clock();
    std::cout << (double)((toc - tic)) / CLOCKS_PER_SEC << std::endl;
}

使用gcc-8.4在Ubuntu-18.04编译生成

gfortran test.f90 -O3 -march=native -o testf
g++ test.cpp -O3 -march=native -I/path/to/eigen -o testcpp

然后我得到了结果：

Fortran with loop:   10.9700003
Fortran with matmul:   0.834999979
Eigen with loop: 38.2188
Eigen with *: 0.40625

内部实现的速度相当，但为什么Eigen在循环实现方面要慢得多？

2
我认为在Fortran中，operator()是一种本机方式来访问内存，但在C++中，它将是operator[]。Eigen中的operator()似乎返回引用对象以访问一个具有重载的operator=和type()的内存，并修改它。这种开销可能是Eigen中循环实现比Fortran慢的原因。 - Semyon Burov
4
"-ffast-math" 命名不准确，应该被命名为 "-funsafe-math"。至少在 Fortran 中，它会破坏符合规范的 Fortran 程序。 - evets
2
@SemyonBurov 可能是可能的（附注：在Fortran中，索引不是运算符），但首先我不明白为什么有人会从Eigen的朴素定义循环计算矩阵乘法。 - Vladimir F Героям слава
3
标题有些误导性，如果Eigen的"Matrix * Matrix"比Fortran的"matmul"更快。如果你认为需要手动访问大量系数，则应考虑使用"-DNDEBUG"来构建（在确保代码按预期运行后）。 - chtz
2
请注意，system_clock不一定以毫秒为单位测量。请查看rate参数。 - Ian Bush
显示剩余7条评论
3个回答
3
循环的最大问题是它们没有按照 C++（应该是行优先）或 Fortran（应该是列优先）的正确顺序进行。这会给你带来很大的性能损失，尤其是在处理大矩阵时。
John Alexiou 的 nativemul 实现（使用 dot_product）也存在同样的问题，所以我非常惊讶他声称它更快。 （而且我发现它不是，见下文。也许他（intel？）编译器重写了代码以内部使用 matmul。）
这是 Fortran 的正确循环顺序：
    c = 0
    do j=1,n
        do k=1,n
            do i=1,n
                c(i,j) = c(i,j) + a(i,k) * b(k,j)
            enddo
        enddo
    enddo

使用 gfortran 版本 10.2.0 编译时，加上 -O3 参数，我得到了以下结果：

 Fortran with original OP's loop:        53.5190010
 Fortran with John Alexiou's nativemul:  53.4309998
 Fortran with correct loop:              11.0679998
 Fortran with matmul:                     2.3699999

在C++中正确的循环应该给你类似的性能。

显然，对于大矩阵，matmul/BLAS要快得多。

Eigen矩阵默认使用列主序存储。您可以通过模板参数选择存储顺序：https://eigen.tuxfamily.org/dox/group__TopicStorageOrders.html。这意味着Fortran naive matmul建议的顺序应该在naive Eigen C++版本中使用。 - IPribec
0
在Fortran代码中，我遇到了同样的问题，但是后来我将矩阵乘法移入了一个子程序，结果速度几乎与matmul一样好。我还与BLAS Level 3函数进行了比较。
Fortran with loop:   9.220000
Fortran with matmul:   8.450000
Fortran with blas3:   2.050000

以及生成它的代码

program ConsoleMatMul
use BLAS95
implicit none
integer :: n,i,j
integer :: tic,toc
real(8),ALLOCATABLE :: a(:,:),b(:,:),c(:,:),xe(:,:)

n = 2000
allocate(a(n,n),b(n,n),c(n,n),xe(n,n))
do i=1,n
    do j =1,n
        a(j,i) = i * 1.0
        b(j,i) = i * 1.0
    enddo
enddo

call system_clock(tic)
call nativemul(a,b,c)
call system_clock(toc)
print*,'Fortran with loop:', (toc - tic) / 1000.0

call system_clock(tic)
c = matmul(a,b)
call system_clock(toc)
print*,'Fortran with matmul:', (toc - tic) / 1000.0
c = b
xe = 0d0
call system_clock(tic)
call gemm(a,c,xe) ! BLAS MATRIX/MATRIX MUL
call system_clock(toc)
print*,'Fortran with blas3:', (toc - tic) / 1000.0

DEALLOCATE(a,b,c)

contains

pure subroutine nativemul(a,b,c)
real(8), intent(in), allocatable :: a(:,:), b(:,:)
real(8), intent(out), allocatable :: c(:,:)
real(8) :: s
integer :: n, i,j,k
    n = size(a,1)
    if (.not. allocated(c)) allocate(c(n,n))
    do j=1,n
        do i=1,n
            s = 0.0d0
            do k=1,n
                s = s + a(i,k) * b(k,j)
            end do
            c(i,j) = s
        end do
    end do
end subroutine    

end program ConsoleMatMul

在我将代码移入子程序之前，我得到了以下结果

Fortran with loop:   85.450000

更新当内部循环被dot_product（）替换时，本地乘法达到或超过matmul水平。

pure subroutine nativemul(a,b,c)
real(8), intent(in) :: a(:,:), b(:,:)
real(8), intent(out) :: c(:,:)
integer :: n, i,j
    n = size(a,1)
    do j=1,n
        do i=1,n
            c(i,j) = dot_product(a(i,:),b(:,j))
            ! or  = sum(a(i,:)*b(:,j))
        end do
    end do
end subroutine
3
为什么你想要nativemul的a和b是可分配的？ - francescalus
1
请注意，您可以使用 BLAS 完成 matmul。 - Vladimir F Героям слава
1
矩阵c在nativemul开始时总是未分配的，因为它被声明为intent(out), allocatable。因此，“Fortran with loop”中的nativemul计时包括重新分配的时间。 - knia
@A.Hennink - 很好的观点。real(8), intent(inout), allocatable :: c(:,:) 可以解决这个问题吗？ - John Alexiou
是的...但你需要的是intent(out)，不需要allocatable。 "实际参数"(在调用过程中放入的内容)是可分配的，但虚拟参数(在程序中声明的参数)不需要是可分配的。 - knia
-1

C++中的前置自增比后置自增更快...

for(int j=0;j<n;++j){
        for(int i=0;i<n;++i){
            double s= 0.0;
            for(int k=0;k<n;++k){
                s += a(i,k) * b(k,j);
            }
            c(i,j) = s;
        }
    }
1
那么，这意味着什么？为什么这个事实很重要？ - Vladimir F Героям слава
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接