个人中心
gpt4 book ai didi

matrix - Fortran中多线程矩阵的转置

转载 作者：行者123 更新时间：2023-12-05 01:06:51 24 4
gpt4 key购买 nike

我正在使用 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 个线程可以让你多一点，但不如上面的改进那么多。

Transpose Speed up relative to intrinsic function for N=20000

关于matrix - Fortran中多线程矩阵的转置，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/68344472/

24 4 0
文章推荐： JupyterHub Helm 安装 : function "dig" not defined
文章推荐： node.js - AWS S3 Node.js SDK `NotImplemented` Multer 错误
文章推荐： null - 检查变量是否为 NA
文章推荐： github-actions - 仅对非草稿拉取请求运行操作
行者123
个人简介

我是一名优秀的程序员,十分优秀！

滴滴打车优惠券免费领取
滴滴打车优惠券
Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
广告合作:1813099741@qq.com 6ren.com