fortran - 带有 BLAS 的 OpenMP-6ren

fortran - 带有 BLAS 的 OpenMP

转载作者：行者123 更新时间：2023-12-04 08:04:29

related question
我试图在上述链接的答案中扩展代码，以包括交叉检查和 openmp。

Program reshape_for_blas

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

  Implicit None

  Real( wp ), Dimension( :, :    ), Allocatable :: a
  Real( wp ), Dimension( :, :, : ), Allocatable :: b
  Real( wp ), Dimension( :, :, : ), Allocatable :: c1, c2, c3, c4, c5
  Real( wp ), Dimension( :, :    ), Allocatable :: d
  Real( wp ), Dimension( :, :    ), Allocatable :: e

  Integer :: na, nb, nc, nd, ne
  Integer :: la, lb, lc, ld  
  Integer( li ) :: start, finish, rate, numthreads

  numthreads = 2

  call omp_set_num_threads(numthreads)  
  
  
  Write( *, * ) 'na, nb, nc, nd ?'
  Read( *, * ) na, nb, nc, nd
  ne = nc * nd
  Allocate( a ( 1:na, 1:nb ) ) 
  Allocate( b ( 1:nb, 1:nc, 1:nd ) ) 
  Allocate( c1( 1:na, 1:nc, 1:nd ) ) 
  Allocate( c2( 1:na, 1:nc, 1:nd ) ) 
  Allocate( c3( 1:na, 1:nc, 1:nd ) )
  Allocate( c4( 1:na, 1:nc, 1:nd ) )
  Allocate( c5( 1:na, 1:nc, 1:nd ) )  
  Allocate( d ( 1:nb, 1:ne ) ) 
  Allocate( e ( 1:na, 1:ne ) ) 

  ! Set up some data
  Call Random_number( a )
  Call Random_number( b )

  ! With reshapes
  Call System_clock( start, rate )
  
  !write (*,*) 'clock', start, rate
  
  
  d = Reshape( b, Shape( d ) )
  Call dgemm( 'N', 'N', na, ne, nb, 1.0_wp, a, Size( a, Dim = 1 ), &
                                            d, Size( d, Dim = 1 ), &
                                    0.0_wp, e, Size( e, Dim = 1 ) )
  c1 = Reshape( e, Shape( c1 ) )
  Call System_clock( finish, rate )
  
  !write (*,*) 'clock', finish, rate
  
  
  
  Write( *, * ) 'Time for reshaping method ', Real( finish - start, wp ) / rate
  Write( *, * ) 'Difference between result matrices ', Maxval( Abs( c1 - c2 ) )
  
  
  ! Direct
  Call System_clock( start, rate )
  Call dgemm( 'N', 'N', na, ne, nb, 1.0_wp, a , Size( a , Dim = 1 ), &
                                            b , Size( b , Dim = 1 ), &
                                            0.0_wp, c2, Size( c2, Dim = 1 ) )
  Call System_clock( finish, rate )
  Write( *, * ) 'Time for straight  method ', Real( finish - start, wp ) / rate
    
    
  
  Call System_clock( start, rate )
  
  !$omp parallel
  ! Direct
  Call dgemm( 'N', 'N', na, ne, nb, 1.0_wp, a , Size( a , Dim = 1 ), &
                                            b , Size( b , Dim = 1 ), &
                                            0.0_wp, c4, Size( c4, Dim = 1 ) )
  !$omp end parallel
  Call System_clock( finish, rate )  
  Write( *, * ) 'Time for straight  method omp', Real( finish - start, wp ) / rate
    
  
 
  !naive
  Call System_clock( start, rate )

  do la = 1, na 
    do lc = 1, nc
      do ld = 1, nd
         c3(la,lc,ld) = 0.0_wp
      enddo
    enddo
  enddo
  
  do la = 1, na 
    do lb = 1, nb
      do lc = 1, nc
        do ld = 1, nd
          c3(la,lc,ld) = c3(la,lc,ld)  + a(la,lb) * b(lb, lc, ld)
        enddo  
      enddo
    enddo
  enddo  
 
  
  Call System_clock( finish, rate )
  Write( *, * ) 'Time for loop', Real( finish - start, wp ) / rate  
 
   

  !naive omp
  Call System_clock( start, rate )  
  !$omp parallel

  do la = 1, na 
    do lc = 1, nc
      do ld = 1, nd
         c5(la,lc,ld) = 0.0_wp
      enddo
    enddo
  enddo
  
  !$omp do private(la, lb, lc, ld) schedule(dynamic) reduction(+: c5)    
  do la = 1, na 
    do lb = 1, nb
      do lc = 1, nc
        do ld = 1, nd
          c5(la,lc,ld) = c5(la,lc,ld)  + a(la,lb) * b(lb, lc, ld)
        enddo  
      enddo
    enddo
  enddo  
  !$omp end do  
  !$omp end parallel
 
  
  Call System_clock( finish, rate )
  Write( *, * ) 'Time for loop omp', Real( finish - start, wp ) / rate  
 
  
  
  
  
  do la = 1, na 
    do lc = 1, nc
      do ld = 1, nd
         
         if ( dabs(c3(la,lc,ld) - c1(la,lc,ld))  > 1.e-6 ) then 
           write (*,*) '!!! c1', c3(la,lc,ld) - c1(la,lc,ld)
         endif         
         
         
         if ( dabs(c3(la,lc,ld) - c2(la,lc,ld))  > 1.e-6 ) then 
           write (*,*) '!!! c2', c3(la,lc,ld) - c2(la,lc,ld)
         endif
         
         if ( dabs(c3(la,lc,ld) - c4(la,lc,ld))  > 1.e-6 ) then 
           write (*,*) '!!! c4', la,lc,ld, c3(la,lc,ld) - c4(la,lc,ld)
         endif
         
         if ( dabs(c3(la,lc,ld) - c5(la,lc,ld))  > 1.e-6 ) then 
           write (*,*) '!!! c5', la,lc,ld, c3(la,lc,ld) - c5(la,lc,ld)
         endif         
      enddo
    enddo
  enddo  
  
End Program reshape_for_blas

我有两个问题:

BLAS 或幼稚循环都没有显着的加速。例如，通过 gfortran -std=f2008 -Wall -Wextra -fcheck=all reshape.f90 -lblas -fopenmp , 并输入 30 30 30 60 , 我得到了

30 30 30 60
 Time for reshaping method    2.9443999999999998E-003
 Difference between result matrices    12.380937791257775
 Time for straight  method    1.0016000000000001E-003
 Time for straight  method omp   2.4878000000000001E-003
 Time for loop   6.6072500000000006E-002
 Time for loop omp  0.100242600000000002

当维度变大时，例如 60 60 60 60在输入中，openmp BLAS 结果可以得到与天真的循环不同的值，似乎我错过了一些控制选项。

OpenMP 在这里会出现什么问题？
编辑
我在 c5 的初始化中添加了 omp 行。部分并注释掉两个打印行，


Program reshape_for_blas

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

  Implicit None

  Real( wp ), Dimension( :, :    ), Allocatable :: a
  Real( wp ), Dimension( :, :, : ), Allocatable :: b
  Real( wp ), Dimension( :, :, : ), Allocatable :: c1, c2, c3, c4, c5
  Real( wp ), Dimension( :, :    ), Allocatable :: d
  Real( wp ), Dimension( :, :    ), Allocatable :: e

  Integer :: na, nb, nc, nd, ne
  Integer :: la, lb, lc, ld  
  Integer( li ) :: start, finish, rate, numthreads

  numthreads = 2

  call omp_set_num_threads(numthreads)  
  
  
  Write( *, * ) 'na, nb, nc, nd ?'
  Read( *, * ) na, nb, nc, nd
  ne = nc * nd
  Allocate( a ( 1:na, 1:nb ) ) 
  Allocate( b ( 1:nb, 1:nc, 1:nd ) ) 
  Allocate( c1( 1:na, 1:nc, 1:nd ) ) 
  Allocate( c2( 1:na, 1:nc, 1:nd ) ) 
  Allocate( c3( 1:na, 1:nc, 1:nd ) )
  Allocate( c4( 1:na, 1:nc, 1:nd ) )
  Allocate( c5( 1:na, 1:nc, 1:nd ) )  
  Allocate( d ( 1:nb, 1:ne ) ) 
  Allocate( e ( 1:na, 1:ne ) ) 

  ! Set up some data
  Call Random_number( a )
  Call Random_number( b )

  ! With reshapes
  Call System_clock( start, rate )
  
  !write (*,*) 'clock', start, rate
  
  
  d = Reshape( b, Shape( d ) )
  Call dgemm( 'N', 'N', na, ne, nb, 1.0_wp, a, Size( a, Dim = 1 ), &
                                            d, Size( d, Dim = 1 ), &
                                    0.0_wp, e, Size( e, Dim = 1 ) )
  c1 = Reshape( e, Shape( c1 ) )
  Call System_clock( finish, rate )
  
  !write (*,*) 'clock', finish, rate
  
  
  
  Write( *, * ) 'Time for reshaping method ', Real( finish - start, wp ) / rate
  Write( *, * ) 'Difference between result matrices ', Maxval( Abs( c1 - c2 ) )
  
  
  ! Direct
  Call System_clock( start, rate )
  Call dgemm( 'N', 'N', na, ne, nb, 1.0_wp, a , Size( a , Dim = 1 ), &
                                            b , Size( b , Dim = 1 ), &
                                            0.0_wp, c2, Size( c2, Dim = 1 ) )
  Call System_clock( finish, rate )
  Write( *, * ) 'Time for straight  method ', Real( finish - start, wp ) / rate
    
    

  !naive loop
  Call System_clock( start, rate )

  do la = 1, na 
    do lc = 1, nc
      do ld = 1, nd
         c3(la,lc,ld) = 0.0_wp
      enddo
    enddo
  enddo
  
  do la = 1, na 
    do lb = 1, nb
      do lc = 1, nc
        do ld = 1, nd
          c3(la,lc,ld) = c3(la,lc,ld)  + a(la,lb) * b(lb, lc, ld)
        enddo  
      enddo
    enddo
  enddo  
 
  
  Call System_clock( finish, rate )
  Write( *, * ) 'Time for loop', Real( finish - start, wp ) / rate  
   



  !dgemm omp 
  Call System_clock( start, rate )
  
  !$omp parallel
  ! Direct
  Call dgemm( 'N', 'N', na, ne, nb, 1.0_wp, a , Size( a , Dim = 1 ), &
                                            b , Size( b , Dim = 1 ), &
                                            0.0_wp, c4, Size( c4, Dim = 1 ) )
  !$omp end parallel
  Call System_clock( finish, rate )  
  Write( *, * ) 'Time for straight  method omp', Real( finish - start, wp ) / rate
    
  
 

  !loop omp
  Call System_clock( start, rate )  
  !$omp parallel

  do la = 1, na 
    do lc = 1, nc
      do ld = 1, nd
         c5(la,lc,ld) = 0.0_wp
      enddo
    enddo
  enddo
  
  !$omp do private(la, lb, lc, ld) schedule(dynamic) reduction(+: c5)    
  do la = 1, na 
    do lb = 1, nb
      do lc = 1, nc
        do ld = 1, nd
          c5(la,lc,ld) = c5(la,lc,ld)  + a(la,lb) * b(lb, lc, ld)
        enddo  
      enddo
    enddo
  enddo  
  !$omp end do  
  !$omp end parallel
 
  
  Call System_clock( finish, rate )
  Write( *, * ) 'Time for loop omp', Real( finish - start, wp ) / rate  
 
  


!single core: c1 c2 c3
! c1 reshape blas
! c2 blas
! c3 naive loop (reference)
! parallel: c4 c5
! c4 dgemm parallel
! c5 naive loop parallel


  do la = 1, na 
    do lc = 1, nc
      do ld = 1, nd
         
         if ( dabs(c3(la,lc,ld) - c1(la,lc,ld))  > 1.e-6 ) then 
           write (*,*) '!!! c1', c3(la,lc,ld) - c1(la,lc,ld)
         endif         
         
         
         if ( dabs(c3(la,lc,ld) - c2(la,lc,ld))  > 1.e-6 ) then 
           write (*,*) '!!! c2', c3(la,lc,ld) - c2(la,lc,ld)
         endif
         
         if ( dabs(c3(la,lc,ld) - c4(la,lc,ld))  > 1.e-6 ) then 
           write (*,*) '!!! c4', la,lc,ld, c3(la,lc,ld) - c4(la,lc,ld)
         endif
         
         if ( dabs(c3(la,lc,ld) - c5(la,lc,ld))  > 1.e-6 ) then 
           write (*,*) '!!! c5', la,lc,ld, c3(la,lc,ld) - c5(la,lc,ld)
         endif         
      enddo
    enddo
  enddo  
  
End Program reshape_for_blas

然后 gfortran reshape.f90 -lblas -fopenmp , 和 30 30 30 30输入导致

 Time for reshaping method    1.3519000000000001E-003
 Difference between result matrices    12.380937791257775
 Time for straight  method    6.2549999999999997E-004
 Time for straight  method omp   1.2600000000000001E-003
 Time for naive loop   1.0008599999999999E-002
 Time for naive loop omp   1.6678999999999999E-002

虽然速度不好。

最佳答案

您调用DGEMM并行使用相同的变量集(因为并行区域中的变量在 Fortran 中默认共享)。由于数据竞争，这不起作用并产生奇怪的结果。你有两个选择:

找到一个并行的 BLAS 实现，其中 DGEMM已经线程了。英特尔 MKL 和 OpenBLAS 是主要候选者。英特尔 MKL 使用 OpenMP，更具体地说，它是使用英特尔 OpenMP 运行时构建的，因此它可能无法很好地处理使用 GCC 编译的 OpenMP 代码，但它可以完美地处理非线程代码。

调用 DGEMM并行但不具有相同的参数集。相反，执行一个或两个张量的 block 分解，并让每个线程为单独的 block 进行收缩。由于 Fortran 使用列优先存储，分解第二个张量可能是合适的:

C[i,k,l=1..L] = A[i,j] * B[j,k,l=1..L]

变成两个线程:

thread 0: C[i,k,l=1..L/2] = A[i,j] * B[j,k,l=1..L/2]
thread 1: C[i,k,l=L/2+1..L] = A[i,j] * B[j,k,l=L/2+1..L]

对于任意数量的线程，归结为计算 l 的开始和结束值。每个线程中的索引并调整 DGEMM 的参数因此。

就个人而言，我会选择并行的 BLAS 实现。使用 Intel MKL，您只需要链接并行驱动程序，它就会自动使用所有可用的 CPU。
下面是 block 分解的示例实现。仅显示对原始代码的添加和更改:

  ! ADD: Use the OpenMP module
  Use :: omp_lib

  ! ADD: Variables used for the decomposition
  Integer :: ithr, istart, iend

  ! CHANGE: OpenMP with block decomposition
  !$omp parallel private(ithr, istart, iend)
    ithr = omp_get_thread_num()

    ! First index (plane) in B for the current thread
    istart = ithr * nd / omp_get_num_threads()
    ! First index (plane) in B for the next thread
    iend = (ithr + 1) * nd / opm_get_num_threads()

    Call dgemm('N', 'N', na, nc * (iend - istart), nb, 1.0_wp, a, nd, &
               b(1, 1, 1 + istart), Size(b, Dim = 1), &
               0.0_wp, c4(1, 1, 1 + istart), Size(c4, Dim = 1))
  !$omp end parallel

istart是 B 的第一个平面的索引每个单独的线程都在工作。 iend是下一个线程的第一个平面，所以 iend - istart是当前线程的平面数。 b(1, 1, 1 + istart)是 B 中的平面 block 开始的地方。 c4(1, 1, 1 + istart)是结果张量中的 block 开始的地方。
确保你做其中之一，但不要同时做。即，如果您的 BLAS 实现是线程化的，但您决定使用 block 分解，请禁用 BLAS 库中的线程化。相反，如果您在 BLAS 实现中使用线程，请不要在代码中执行 block 分解。

关于fortran - 带有 BLAS 的 OpenMP，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/66296334/

文章推荐： javascript - Asp.Net core Vue C# + JS动态添加输入文本框

文章推荐： C# 从 HTML Body 解析 XML 并保存到文件

文章推荐： java - 禁用错误消息 : AccessDeniedException: Access is denied

blas - 使用 100M 行转置 CSR 时出现 MKL 稀疏 BLAS 段错误
我正在尝试将 MKL 稀疏 BLAS 用于 CSR 矩阵，行数/列数约为 100M。当我将其增加到 100M 时，我的源代码似乎适用于 10M 行/列，但由于段错误而失败。我将问题隔离到以下代码片段
scala - 为什么 spark blas 在 1 级例程中使用 f2jBLAS 而不是原生 BLAS？
我在 BLAS.scala 中找到了以下代码: // For level-1 routines, we use Java implementation. private def f2jBLAS: Ne
c++ - Mac OS 中的 BLAS/LAPACK 和原来的 BLAS/LAPACK 有什么区别(如果有的话)？
我最近从 Linux 切换到 Mac OS。我需要 BLAS 和 LAPACK 来做一些计算。通过查看 BLAS 的维基百科，我了解到这两个库已经在 Mac OS 中实现。不过，据说 Apple's
c++ - Armadillo C++ 和 BLAS 和 ATLAS 在 mingw32 下找不到 lapack blas
全部，总结... 我正在尝试编译 example.cpp arma.sourceforge.net/docs.html#example_prog 我尝试用 lapack 和 blas 编译和链接 A
r - 在没有 root 访问权限的情况下，当 R 与引用 BLAS 链接时，使用调整后的 BLAS 运行
谁能告诉我为什么我可以不成功测试 OpenBLAS dgemm通过以下方式在 R 中性能(在 GFLOP 中)？将 R 与“引用 BLAS”链接起来 libblas.so 编译我的 C 程序 mmp
blas - BLAS例程中的incx
有某些BLAS例程将向量X的增量即incX作为参数。我找不到增量，以及增量如何影响计算结果。任何人都可以提供示例或任何其他类型的信息吗？更新: 我在这里找到了最好的信息: Intel HPC m
blas - 为什么加减法没有BLAS例程
在 BLAS 中有这样的例程 dscal scale a vector by a constant dinit initialize a vector with given value
blas - blas中incx和incxy是什么意思？
对于 BLAS 函数 sdot (n, x, incx, y, incy)。 incx 指定 x 元素的增量。参数incx和incy是什么意思？最佳答案 sdot (n, x, incx, y,
matrix - BLAS 矩阵乘矩阵转置乘法
我必须以 A'A 的形式计算一些产品或更一般的 A'DA ，其中 A是将军mxn矩阵和 D是对角线 mxm矩阵。他们都是满级；即 rank(A)=min(m,n) . 我知道你可以节省大量时间是这样的
benchmarking - BLAS 库基准
是否有比较不同 BLAS(基本线性代数子程序)库的基准？我对单核和多核系统的稀疏矩阵乘法特别感兴趣？最佳答案 BLAS 性能在很大程度上取决于系统，因此您最好在要使用的机器上自己进行基准测试。由于只
optimization - 优化的汇编代码 (BLAS)
因此，关于通过汇编代码提高性能的问题的答案通常是“别打扰，编译器比你更聪明”。我明白了。但是，我注意到优化的线性代数库(例如 ACML)可以实现比标准编译库高 2 到 5 倍的性能改进。例如，在我的
java - BLAS.dgemm方法多线程计算误差
当我使用spark mllib多层感知器模型来预测 vector 时，我发现同一 vector 在多线程中有时会给出不同的结果。我阅读了源代码，发现它是基于BLAS lib的。我为BLAS在多线程中编
c - 缺少头文件时链接到 BLAS
我正在尝试用 C 语言编译一个程序，该程序使用线性代数的 BLAS 接口(interface)。该系统在 /usr/lib64/libblas.* 中具有 BLAS 库(.a 和 .so 文件)但没有
c++ - Blas 看起来很慢
我在我的机器上比较矩阵乘法，似乎 c++ blas 非常慢。一个1000x1000的矩阵相乘大约需要4秒，而在python中同样需要1.5秒左右。我认为链接可能有问题，但我真的不知道如何解决这些问题。
c++ - BLAS 中矩阵之间的元素明智乘法？
我开始使用 C++(特别是英特尔 MKL)中的 BLAS 函数来创建我的一些旧 Matlab 代码的更快版本。到目前为止它运行良好，但我无法弄清楚如何对 2 个矩阵(Matlab 中的 A.* B)
c++ - BLAS 是如何获得如此极致的性能的？
出于好奇，我决定对我自己的矩阵乘法函数与 BLAS 实现进行基准测试......我对结果最不感到惊讶: Custom Implementation, 10 trials of 1000x1000 ma
fortran - 在 BLAS 中转置还是先自己做？
我在 Fortran 77 中整理了一些科学代码，我正在争论什么会更快。基本上，我有一个 MxN 矩阵，我们称之为 A。M 大于 N。稍后在代码中，我需要将 transpose(A) 乘以一堆向量。
language-agnostic - BLAS 如何结合矩阵链乘法优化
BLAS(基本线性代数子程序)提供了许多其他编程语言，比如我使用的 Matlab，以及快速例程来执行矩阵乘法等操作。然而，当将多个矩阵相乘时，有一个最佳顺序来“括号”矩阵。取自 wikipedia
boost - BLAS 和 CUBLAS
我想知道 NVIDIA 的 cuBLAS 库。有没有人有这方面的经验？例如，如果我使用 BLAS 编写一个 C 程序，我是否能够用对 cuBLAS 的调用替换对 BLAS 的调用？或者甚至更好地实现一
matrix - 如何使用 BLAS 执行向量矩阵乘法？
BLAS 定义了 GEMV(矩阵向量乘法)2 级运算。如何使用 BLAS 库执行向量矩阵乘法？这可能很明显，但我不知道如何使用 BLAS 运算进行乘法运算。我本来希望进行 GEVM 操作。最佳答案

行者123

个人简介

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

作者热门文章

滴滴打车优惠券免费领取

全站热门文章

首页

博学

6Ren·AI

商城

fortran - 带有 BLAS 的 OpenMP