使用 OMP 加速对称矩阵的计算

Question

我的矩阵计算是：C=C-A*B

这里 C 是一个对称矩阵，所以我想通过只考虑上三角然后取相反的元素来加快这个计算。我使用了 OMP，发现我的实现比整个矩阵 C 的正常计算要慢。

我还看到 C=C-AxB 的计算比 C=C+AxB 慢。

附上我的程序。请多多指教！

    Program testspeed
implicit none
integer nstate,nmeas,i,j,l
integer(kind=8) :: tclock1, tclock2, clock_rate
real(kind=8) :: elapsed_time
double precision, allocatable, dimension(:,:):: B,C,A
nstate =20000
nmeas=10000
allocate (B(nmeas,nstate),C(nstate,nstate),A(nstate,nmeas))
A=1d0
B=1d0
call system_clock(tclock1)
write(*,*) "1"
!$omp parallel do
do j = 1, nstate
    do l = 1,nmeas
        do i = 1, j
            C(j,i) = C(j,i) - A(j,l)*B(l,i)
            C(i,j)=C(j,i)
        end do
    end do
end do
!$omp end parallel do
write(*,*) "2"
call system_clock(tclock2, clock_rate)
elapsed_time = float(tclock2 - tclock1) / float(clock_rate)
write(*,*) elapsed_time
end Program testspeed

Answer 1

我教给我的学生的基本规则之一是，现在没有人应该编写自己的密集矩阵乘法运算——而且 30 多年来都不应该这样做。您应该改用 BLAS 库。下面我将使用 BLAS 库与您的循环排序和更好的循环排序进行比较，并与我用作参考的 Fortran 内部函数 matmul 进行比较，以检查结果是否正确。 BLAS 和 matmul 没有利用 C 的对称性，但它们仍然是最快的例程 - BLAS 比您编写的循环顺序快大约 200-300 倍。请注意，我还减少了一些矩阵大小，因为我厌倦了等待原件运行对于更大的情况：

ijb@ijb-Latitude-5410:~/work/stack$ cat mm.f90
Program testspeed

  Use, Intrinsic :: iso_fortran_env, Only : wp => real64, li => int64

  Use omp_lib, Only : omp_get_max_threads
  
  Implicit None
  Integer nstate,nmeas,i,j,l
  Integer(li) :: tclock1, tclock2, clock_rate
  Real(wp) :: elapsed_time
  Real( wp ), Allocatable, Dimension(:,:):: B,C,A
  Real( wp ), Allocatable, Dimension(:,:):: C_test
  Real( wp ), Allocatable, Dimension(:,:):: C_start

  Write( *, * ) 'Using ', omp_get_max_threads(), ' threads'
  
!!$  nstate =2000
!!$  nmeas=1000
  nstate = 5000
  nmeas  = 2500
  Allocate (B(nmeas,nstate),C(nstate,nstate),A(nstate,nmeas))
  Allocate( C_test, Mold = C )
  Allocate( C_start, Mold = C )
!!$  A=1.0_wp
!!$  B=1.0_wp
  ! Random numbers are a much better test
  Call Random_number( A )
  B = Transpose( A ) ! make sure result is symmetric
  Call Random_number( C_start )
  ! Make Initial C Symmetric
  C_start = 0.5_wp * ( C_start + Transpose( C_start ) )

  Write( *, * ) 'Matix sizes ', Shape( A ), Shape( B ), Shape( C )
  
  C_test = C_start
  Call system_Clock(tclock1)
  C_test = C_test -  Matmul( A, B )
  Call system_Clock(tclock2, clock_rate)
  elapsed_time = Real(tclock2 - tclock1,wp) / Real(clock_rate,wp)
  Write( *,'( a, t20, f8.3 )' ) 'Matmul', elapsed_time

  C = C_start
  Call system_Clock(tclock1)
  !$omp parallel do
  Do j = 1, nstate
     Do l = 1,nmeas
        Do i = 1, j
           C(j,i) = C(j,i) - A(j,l)*B(l,i)
           C(i,j)=C(j,i)
        End Do
     End Do
  End Do
  !$omp end parallel do
  Call system_Clock(tclock2, clock_rate)
  elapsed_time = Real(tclock2 - tclock1,wp) / Real(clock_rate,wp)
  Write(*,'( a, t20, f8.3, t30, "Max error ", g20.12 )' ) &
       'Orig loops', elapsed_time, Maxval( Abs( C_test - C ) )

  C = C_start
  Call system_Clock(tclock1)
  !$omp parallel default( none ) shared ( nstate, nmeas, A, B, C ), private( i, j, l )
  !$omp do 
  Do i = 1, nstate
     Do l = 1,nmeas
        Do j = 1, i
           C(j,i) = C(j,i) - A(j,l)*B(l,i)
        End Do
     End Do
  End Do
  !$omp end do
  !$omp do
  Do i = 1, nstate
     Do j = 1, i
        C( i, j ) = C( j, i )
     End Do
  End Do
  !$omp end do
  !$omp end parallel
  Call system_Clock(tclock2, clock_rate)
  elapsed_time = Real(tclock2 - tclock1,wp) / Real(clock_rate,wp)
  Write(*,'( a, t20, f8.3, t30, "Max error ", g20.12 )' ) &
       'Sensible loops', elapsed_time, Maxval( Abs( C_test - C ) )

  C = C_start
  Call system_Clock(tclock1)
  Call dgemm( 'N', 'N', nstate, nstate, nmeas, -1.0_wp, A, Size( A, Dim = 1 ), &
                                                        B, Size( B, Dim = 1 ), &
                                               +1.0_wp, C, Size( C, Dim = 1 ) )
  Call system_Clock(tclock2, clock_rate)
  elapsed_time = Real(tclock2 - tclock1,wp) / Real(clock_rate,wp)
  Write(*,'( a, t20, f8.3, t30, "Max error ", g20.12 )' ) &
       'BLAS ', elapsed_time, Maxval( Abs( C_test - C ) )
  
End Program testspeed
ijb@ijb-Latitude-5410:~/work/stack$ gfortran --version
GNU Fortran (Ubuntu 9.3.0-17ubuntu1~20.04) 9.3.0
Copyright (C) 2019 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

ijb@ijb-Latitude-5410:~/work/stack$ gfortran -fopenmp -Wall -Wextra -std=f2018 -O3  mm.f90 -lopenblas
ijb@ijb-Latitude-5410:~/work/stack$ export OMP_NUM_THREADS=1
ijb@ijb-Latitude-5410:~/work/stack$ ./a.out
 Using            1  threads
 Matix sizes         5000        2500        2500        5000        5000        5000
Matmul                4.793
Orig loops          421.564  Max error   0.488853402203E-11
Sensible loops       20.742  Max error   0.488853402203E-11
BLAS                  2.185  Max error   0.682121026330E-12
ijb@ijb-Latitude-5410:~/work/stack$ export OMP_NUM_THREADS=2
ijb@ijb-Latitude-5410:~/work/stack$ ./a.out
 Using            2  threads
 Matix sizes         5000        2500        2500        5000        5000        5000
Matmul                4.968
Orig loops          324.319  Max error   0.466116034659E-11
Sensible loops       17.656  Max error   0.466116034659E-11
BLAS                  1.161  Max error   0.682121026330E-12
ijb@ijb-Latitude-5410:~/work/stack$ export OMP_NUM_THREADS=3
ijb@ijb-Latitude-5410:~/work/stack$ ./a.out
 Using            3  threads
 Matix sizes         5000        2500        2500        5000        5000        5000
Matmul                4.852
Orig loops          243.268  Max error   0.500222085975E-11
Sensible loops       15.802  Max error   0.500222085975E-11
BLAS                  0.852  Max error   0.682121026330E-12
ijb@ijb-Latitude-5410:~/work/stack$ export OMP_NUM_THREADS=4
ijb@ijb-Latitude-5410:~/work/stack$ ./a.out
 Using            4  threads
 Matix sizes         5000        2500        2500        5000        5000        5000
Matmul                4.994
Orig loops          189.189  Max error   0.477484718431E-11
Sensible loops       14.245  Max error   0.477484718431E-11
BLAS                  0.707  Max error   0.682121026330E-12

对于 BLAS，我使用了 openblas - 它是免费提供的。在 Linux 系统上，一个简单的 apt get 或类似的东西就足够了。

另请注意

如果您必须编写自己的循环，您的最内层循环应该尽可能地遍历数组的第一个索引。这是因为 Fortran 将其数组排序为主要列。这就是我在“明智的”循环排序中所做的
Real( 8 ) 不可移植，不能保证你的编译器支持，也不能保证做你期望的，不应该被使用。 Integer( 8 ) 类似。请查看我为获得适用于所有编译器的更好方法所做的工作。
Float 不是标准内在函数 - 使用 Real 正如我所做的那样
如果结果不正确，基准测试就毫无意义，因此您应该始终包含一种测试结果的方法。这里我使用 Fortran intrinsic matmul 提供一个参考版本。您的原始代码没有初始化 C，因此结果不可信 - 但由于您不检查您是否获得了 C 的正确值，您无法知道这一点。
我个人非常不喜欢 !$omp parallel do，我认为 OpenMP 中有这样的捷径是个错误。而是将它们分成 !$omp parallel 和 !$omp do - 了解线程创建和工作共享是不同的东西非常重要，将它们卷成一行并不是学习 OpenMP 的好方法。

使用 OMP 加速对称矩阵的计算

Speedup of calculation for symmetric matrix using OMP

fortran

openmp

transpose

matrix-multiplication

intel-fortran