我正在使用 Fortran 计算一个非常大的维度矩阵的转置，P=TRANSPOSE(PP)。我看到 Fortran 中的内置函数 TRANSPOSE 很慢。我想通过并行化代码来加快速度:
subroutine TP(nstate,P,PP)
integer i,j,nstate
double precision P(nstate,nstate),PP(nstate,nstate)
!$omp parallel shared ( P, PP,nstate) private (i, j)
!$omp do
do i=1,nstate
do j=1,nstate
P(j,i) = PP(i,j)
end do
end do
!$omp end do
!$omp end parallel do
end subroutine TP(nstate,P,PP)
不幸的是，我的代码只使用了 1 个线程并且没有任何改进。
最佳答案
要回答您关于线程的问题，最可能的原因是您必须设置编译过程才能使用 OpenMP。我所知道的所有编译器不默认使用 OpenMP，您必须明确打开它。
但是，您真正的问题是关于加快矩阵转置的速度，如果我正在执行此多线程处理，那么我将关注的第三件事。第一个是，如评论中所述，我可以在不明确形成转置的情况下重写我的算法吗？根据我的经验，您很少需要明确地形成转置。例如，许多 BLAS 和 LAPACK 例程允许您指定要使用矩阵的转置形式，因此您实际上不必进行转置。当然你能不能做取决于你在做什么的具体细节，但不做几乎肯定是最快的方法!
接下来我会看看单线程性能。矩阵转置是现代内存子系统可以做的最糟糕的事情之一，因此投入更多内核不太可能改善事情，因为限制因素是内存的速度，而不是你可以多快进行计算。现在对于小于缓存大小的小矩阵，这可能不是问题 - 缓存会拯救你。但是对于比缓存更大的大型矩阵，您会遇到问题。解决这个问题的方法是将操作分成更小的 block 并依次转置每个 block - 并选择 block 大小使其小于缓存。下面是一个以两种略有不同的方式执行此操作的代码，以及我的笔记本电脑上的一些单线程结果
! transpose.f90
Program tran
Use, Intrinsic :: iso_fortran_env, Only : wp => real64, li => int64
Use omp_lib, Only : omp_get_max_threads
Implicit None ( type, external )
Real( wp ), Dimension( :, : ), Allocatable :: a, aT
Real( wp ) :: t, t_in, t_in1, t_in2
Real( wp ) :: error
Integer :: start, finish, rate
Integer :: bfac
Integer :: n
Write( *, * ) 'n ?'
Read ( *, * ) n
Allocate( a ( 1:n, 1:n ) )
Allocate( aT( 1:n, 1:n ) )
Call Random_number( a )
Write( *, * ) 'Size of matrix ', n
Write( *, * ) 'Using ', omp_get_max_threads(), ' threads'
Call System_clock( start, rate )
aT = Transpose( a )
Call System_clock( finish, rate )
error = Maxval( Abs( Transpose( a ) - aT ) )
t_in1 = Real( finish - start, wp ) / rate
Write( *, * ) 'Intrinsic 1: ', t_in1, ' Error: ', error
Call System_clock( start, rate )
aT = Transpose( a )
Call System_clock( finish, rate )
error = Maxval( Abs( Transpose( a ) - aT ) )
t_in2 = Real( finish - start, wp ) / rate
Write( *, * ) 'Intrinsic 2: ', t_in2, ' Error: ', error
t_in = 0.5_wp * ( t_in1 + t_in2 )
Write( *, * ) 'Read bound'
bfac = 10
Do While( bfac <= n )
aT = -2.0_wp
Call System_clock( start, rate )
Call blocked_transpose_rb( a, bfac, aT )
Call System_clock( finish, rate )
error = Maxval( Abs( Transpose( a ) - aT ) )
t = Real( finish - start, wp ) / rate
Write( *, * ) bfac, ' Intrinsic: ', t, ' Error: ', error, ' Speed up ', t_in / t
bfac = Nint( bfac * 1.5_wp )
End Do
Write( *, * ) 'Write bound'
bfac = 10
Do While( bfac <= n )
aT = -2.0_wp
Call System_clock( start, rate )
Call blocked_transpose_wb( a, bfac, aT )
Call System_clock( finish, rate )
error = Maxval( Abs( Transpose( a ) - aT ) )
t = Real( finish - start, wp ) / rate
Write( *, * ) bfac, ' Intrinsic: ', t, ' Error: ', error, ' Speed up: ', t_in / t
bfac = Nint( bfac * 1.5_wp )
End Do
Contains
Subroutine blocked_transpose_rb( a, bfac, aT )
Use, Intrinsic :: iso_fortran_env, Only : wp => real64, li => int64
Real( wp ), Dimension( :, : ), Intent( In ) :: a
Integer , Intent( In ) :: bfac
Real( wp ), Dimension( :, : ), Intent( Out ) :: aT
Integer :: n
Integer :: ib, jb
Integer :: i , j
n = Ubound( a, Dim = 1 )
!$omp parallel default( none ) shared( n, bfac, a, AT ), private( ib, jb, i, j )
!$omp do collapse( 2 )
Do ib = 1, n, bfac
Do jb = 1, n, bfac
Do i = ib, Min( n, ib + bfac - 1 )
Do j = jb, Min( n, jb + bfac - 1 )
aT( j, i ) = a( i, j )
End Do
End Do
End Do
End Do
!$omp end do
!$omp end parallel
End Subroutine blocked_transpose_rb
Subroutine blocked_transpose_wb( a, bfac, aT )
Use, Intrinsic :: iso_fortran_env, Only : wp => real64, li => int64
Real( wp ), Dimension( :, : ), Intent( In ) :: a
Integer , Intent( In ) :: bfac
Real( wp ), Dimension( :, : ), Intent( Out ) :: aT
Integer :: n
Integer :: ib, jb
Integer :: i , j
n = Ubound( a, Dim = 1 )
!$omp parallel default( none ) shared( n, bfac, a, AT ), private( ib, jb, i, j )
!$omp do collapse( 2 )
Do ib = 1, n, bfac
Do jb = 1, n, bfac
Do i = ib, Min( n, ib + bfac - 1 )
Do j = jb, Min( n, jb + bfac - 1 )
aT( i, j ) = a( j, i )
End Do
End Do
End Do
End Do
!$omp end do
!$omp end parallel
End Subroutine blocked_transpose_wb
End Program tran
编译
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 -O3 -fopenmp -Wall -Wextra -g transpose.f90 -o transpose
代码执行
ijb@ijb-Latitude-5410:~/work/stack$ export OMP_NUM_THREADS=1; ./transpose < in
n ?
Size of matrix 20000
Using 1 threads
Intrinsic 1: 9.4619999999999997 Error: 0.0000000000000000
Intrinsic 2: 8.6950000000000003 Error: 0.0000000000000000
Read bound
10 Intrinsic: 2.2450000000000001 Error: 0.0000000000000000 Speed up 4.0438752783964365
15 Intrinsic: 2.0019999999999998 Error: 0.0000000000000000 Speed up 4.5347152847152854
23 Intrinsic: 1.8210000000000000 Error: 0.0000000000000000 Speed up 4.9854475562877543
35 Intrinsic: 1.4319999999999999 Error: 0.0000000000000000 Speed up 6.3397346368715084
53 Intrinsic: 1.2629999999999999 Error: 0.0000000000000000 Speed up 7.1880443388756934
80 Intrinsic: 1.2470000000000001 Error: 0.0000000000000000 Speed up 7.2802726543704885
120 Intrinsic: 1.1870000000000001 Error: 0.0000000000000000 Speed up 7.6482729570345409
180 Intrinsic: 1.1990000000000001 Error: 0.0000000000000000 Speed up 7.5717264386989154
270 Intrinsic: 1.1750000000000000 Error: 0.0000000000000000 Speed up 7.7263829787234037
405 Intrinsic: 1.1510000000000000 Error: 0.0000000000000000 Speed up 7.8874891398783662
608 Intrinsic: 1.1319999999999999 Error: 0.0000000000000000 Speed up 8.0198763250883403
912 Intrinsic: 1.2200000000000000 Error: 0.0000000000000000 Speed up 7.4413934426229513
1368 Intrinsic: 1.4050000000000000 Error: 0.0000000000000000 Speed up 6.4615658362989326
2052 Intrinsic: 3.2240000000000002 Error: 0.0000000000000000 Speed up 2.8159119106699748
3078 Intrinsic: 3.8510000000000000 Error: 0.0000000000000000 Speed up 2.3574396260711503
4617 Intrinsic: 4.0990000000000002 Error: 0.0000000000000000 Speed up 2.2148084898755793
6926 Intrinsic: 4.8529999999999998 Error: 0.0000000000000000 Speed up 1.8706985369874305
10389 Intrinsic: 5.5330000000000004 Error: 0.0000000000000000 Speed up 1.6407916139526477
15584 Intrinsic: 5.9450000000000003 Error: 0.0000000000000000 Speed up 1.5270815811606391
Write bound
10 Intrinsic: 1.5669999999999999 Error: 0.0000000000000000 Speed up: 5.7935545628589660
15 Intrinsic: 1.5389999999999999 Error: 0.0000000000000000 Speed up: 5.8989603638726447
23 Intrinsic: 1.4190000000000000 Error: 0.0000000000000000 Speed up: 6.3978153629316417
35 Intrinsic: 1.2110000000000001 Error: 0.0000000000000000 Speed up: 7.4966969446738227
53 Intrinsic: 1.3109999999999999 Error: 0.0000000000000000 Speed up: 6.9248665141113657
80 Intrinsic: 1.0700000000000001 Error: 0.0000000000000000 Speed up: 8.4845794392523359
120 Intrinsic: 1.0320000000000000 Error: 0.0000000000000000 Speed up: 8.7969961240310077
180 Intrinsic: 1.1960000000000000 Error: 0.0000000000000000 Speed up: 7.5907190635451505
270 Intrinsic: 1.2350000000000001 Error: 0.0000000000000000 Speed up: 7.3510121457489870
405 Intrinsic: 1.2480000000000000 Error: 0.0000000000000000 Speed up: 7.2744391025641022
608 Intrinsic: 1.2849999999999999 Error: 0.0000000000000000 Speed up: 7.0649805447470824
912 Intrinsic: 1.4750000000000001 Error: 0.0000000000000000 Speed up: 6.1549152542372880
1368 Intrinsic: 1.6970000000000001 Error: 0.0000000000000000 Speed up: 5.3497348261638180
2052 Intrinsic: 2.5249999999999999 Error: 0.0000000000000000 Speed up: 3.5954455445544555
3078 Intrinsic: 2.9569999999999999 Error: 0.0000000000000000 Speed up: 3.0701724721001016
4617 Intrinsic: 3.1490000000000000 Error: 0.0000000000000000 Speed up: 2.8829787234042552
6926 Intrinsic: 3.6120000000000001 Error: 0.0000000000000000 Speed up: 2.5134274640088594
10389 Intrinsic: 3.7269999999999999 Error: 0.0000000000000000 Speed up: 2.4358733565870674
15584 Intrinsic: 4.4550000000000001 Error: 0.0000000000000000 Speed up: 2.0378226711560044
引用的加速是与 Transpose() 相比，给定阻塞因子的代码运行速度快了多少，因此加速 2 意味着我的代码运行速度是内在代码的两倍，你可以看到，仅仅通过重组循环，你就可以获得 8+ 的改进 [是的，我知道我对这里的内在函数有点不公平，但重点仍然存在]；这相当于大量的线程!此外，如果您查看应该最好地反射(reflect)缓存使用情况的读取绑定(bind)版本，您可以看到最大速度为 608 的阻塞因子。鉴于每个转置需要两个大小的 block (阻塞因子)*(阻塞因子)，这意味着对于 bfac=608，它需要 (608*608*8*2)/(1024*1024)=5.64 MB，略低于我机器的 L3 缓存大小(8MB)。因此，在这种情况下，考虑如何使用内存比考虑如何使用内核更重要。
当然你也可以多线程。下面是我的笔记本电脑上的结果，同样是相对于内在的加速。改用 2 个线程可以让你多一点，但不如上面的改进那么多。
关于matrix - Fortran中多线程矩阵的转置，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/68344472/