LAPACK/BLAS sgemm() 比自定义矩阵乘法慢

Question

我刚刚使用 LAPACK/BLAS 运行了一次矩阵 * 矩阵乘法，并使用自定义循环优化（平铺）运行了一次。我有点恼火，因为简单的循环平铺方法比 BLAS 算法快大约 43%。基本上，我的问题是我在应用 BLAS 例程时是否犯了错误。这是我的代码：

program test 
  implicit none

  integer, parameter :: N = 1000, tile = 2
  real*4, dimension(N,N) :: a,b,c,temp 
  integer :: i,j,k,x,y,z
  double precision ::  E,S
  real :: alpha = 1.0, beta = 0.0

  call random_seed()
  call random_number(a)
  call random_number(b)

  call cpu_time(S)

  ! call sgemm('n','n',N, N, N, alpha,a,N,b,N, beta,c,N)

  do j = 1,N,tile
     do k = 1,N,tile
        do i = 1,N,tile
           do y = j, min( j+tile-1,N)
              do x = i, min( i+tile-1,N)
                 do z = k, min( k+tile-1,N)
                    c(x,y) = c(x,y) + a(x,z) * b(z,y)
                 enddo
              enddo
           enddo
        enddo
     enddo
  enddo

  call cpu_time(E)
  print*,(E-S)
end program test

我在配备 4gb DRAM 和 3096kb 缓存的 Intel Dual Core2 机器上运行此计算。该程序编译为：

$gfortran -O3 test.f03 -o test
0.9359

循环和：

$gfortran test.f03 -lblas -O3 -o test
1.3399

那么我是不是没有得到关于 BLAS 的一些信息，我是否遗漏了一些东西（编译器优化，或者我只是不知道什么）？我在有和没有 Eigen::Matrix 的情况下使用 C++ 运行了类似的代码，并从使用 Eigen 库的 MMM 获得了可观的收益，这就是为什么我的期望与 BLAS 库相似的原因。

Answer 1

正确使用了 BLAS 例程。 唯一的区别是 BLAS 正在执行

C = 0.0*C + 1.0*A*B

和你的循环

C = C + A*B

在您的循环中，您正在尝试提高 cpu 缓存内存的使用率。 有执行类似操作的 BLAS 变体。 我建议您尝试使用 openblas、atlas 或 mkl（英特尔编译器）库。您将获得极大的时间改进。