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

oknwwptz 于 2023-01-28 发布在其他

关注(0)|答案(3)|浏览(224)

我尝试通过Eigen库用一个2000*2000的矩阵乘法实现将Fortran代码重写为C++。我发现Eigen中的for循环比Fortran中的do循环要慢得多（〉3x）。代码如下所示：
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的循环实现要慢得多？

c++

来源：https://stackoverflow.com/questions/65367088/eigen-is-much-slower-than-fortran-in-matrix-multiplication-using-an-explicit-loo

3条答案

按热度按时间

kqlmhetl1#

循环的最大问题是，无论是C++（应该是行为主）还是Fortran（应该是列为主），循环的顺序都不正确，这会给性能带来很大的影响，尤其是对于大型矩阵。
JohnAlexiou的nativemul实现（使用dot_product）也有同样的问题，所以我很惊讶他声称它更快（我发现它不是;见下文。也许他的（英特尔？）编译器重写代码以在内部使用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

使用gfortran10.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对于大型矩阵要快得多。

赞(0）回复(0）举报 2023-01-28

tquggr8v2#

在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

赞(0）回复(0）举报 2023-01-28

z4iuyo4d3#

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;
        }
    }

赞(0）回复(0）举报 2023-01-28

我来回答

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

3条答案

相关问题

热门标签

最新问答