matrix - Fortran中多线程矩阵的转置
问题描述
我正在使用 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 例程允许您指定要使用矩阵的转置形式,因此您实际上不必进行转置。当然你能不能做取决于你正在做什么的确切细节,但不做几乎可以肯定是最快的方法!
接下来我将看看单线程性能。矩阵转置是现代内存子系统可以做的最糟糕的事情之一,因此投入更多内核不太可能改善事情,因为限制因素是内存的速度,而不是你可以多快进行计算。现在对于小于缓存大小的小矩阵,这可能不是问题 - 缓存会拯救你。但是对于比缓存更大的大型矩阵,您会遇到问题。解决这个问题的方法是将操作分成更小的块并依次转置每个块 - 并选择块大小使其小于缓存。下面是一个以两种略有不同的方式执行此操作的代码,以及我的笔记本电脑上的一些单线程结果
! 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+ 的改进 [是的,我知道我对这里是固有的,但重点仍然存在];这相当于大量的线程!此外,如果您查看应该最好地反映缓存使用的读取绑定版本,您可以看到最大速度是 608 的阻塞因子。给定每个转置需要两个大小的块(阻塞因子)*(阻塞因子)这意味着对于 bfac=608,它需要 (608*608*8*2)/(1024*1024)=5.64 MB,略低于我机器的 L3 缓存大小(8MB)。因此,在这种情况下,考虑如何使用内存比如何使用内核更重要。
当然你也可以多线程。下面是我的笔记本电脑上的结果,同样是相对于内在的加速。移动到 2 个线程可以让你多一点,但不如上面的改进那么多。
推荐阅读
- playwright - 如何在禁用网络安全的情况下执行剧作家代码生成
- php - Facebook 总是重定向到 Web App 的主页 - PHP
- mongodb - 在将在 Jaspersoft Studio 中执行的 $match 中使用 $if 或 $ifNull 条件
- git - Terraform 会根据运行时环境不断更改多行 heredoc 的行尾
- arrays - 我正在尝试将列表中的元素加在一起,但不断得到一个奇怪的数字
- popup - 您可以禁用 Microsoft Edge-chromium 上弹出的还原页面吗?
- java - oracle.jdbc.autoCommitSpecCompliant VS setAutoCommit(false)
- swiftui - 如何在 SwiftUI watchOS 上创建闹钟
- linux - 删除/清除活动的 cron 作业
- salesforce - dropbox-api get_thumbnail_v2 & get_thumbnail 在 rhombs () 中返回问号 (?)