指向 fortran 中 OpenBLAS 子例程的指针

Pointer to OpenBLAS subroutines in fortran

提问人：Navid Khairdast 提问时间：11/17/2023 更新时间：11/17/2023 访问量：45

问：

在我的代码中，我想设置指向 OpenBLAS 子例程的指针，以单精度或双精度编译代码。

为此，我定义了两个模块，并为单精度 （sgemv） 或双精度 （dgemv） OpenBLAS 函数定义了函数接口：


module DoublePrecision  
   implicit none 
   integer, parameter                                :: precision = selected_real_kind(15, 307)

   external dgemv

   abstract interface

    subroutine gemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy) 
      
      import                                         :: precision
      character                                      :: trans
      integer(4), value                              :: m, n, lda, incx, incy
      real(precision), value                         :: alpha, beta
      real(precision), dimension(lda, *), intent(in) :: a
      real(precision), dimension(*), intent(in)      :: x
      real(precision), dimension(*), intent(inout)   :: y

    end subroutine gemv

   end interface 

   procedure(gemv), pointer                          :: MatrixVectorMultiplication => null()

end module DoublePrecision

和

module SinglePrecision  
   implicit none 
   integer, parameter                                :: precision = selected_real_kind(6, 37)
   
   external sgemv

   abstract interface

    subroutine gemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
      
      import                                         :: precision
      character                                      :: trans
      integer(4), value                              :: m, n, lda, incx, incy
      real(precision), value                         :: alpha, beta
      real(precision), dimension(lda, *), intent(in) :: a
      real(precision), dimension(*), intent(in)      :: x
      real(precision), dimension(*), intent(inout)   :: y

    end subroutine gemv

   end interface 

   procedure(gemv), pointer                          :: MatrixVectorMultiplication => null()

end module SinglePrecision

在主程序中，我在编译命令中取了一个标志，并使用了相应的模块：

program test_MatrixMultiplication


#ifdef single
   use SinglePrecision
   MatrixVectorMultiplication => sgemv
#else
   use DoublePrecision
   MatrixVectorMultiplication => dgemv
#endif


end program test_MatrixMultiplication

为了编译，我使用以下命令：

! for single precision
gfortran -Dsingle -cpp -L/usr/local/opt/openblas/lib/ -I/usr/local/opt/openblas/include/ -lopenblas -ffree-line-length-1024 SingleDoublePrecisionMatrixMultiplication.f90

! or for double precision
gfortran -Ddouble -cpp -L/usr/local/opt/openblas/lib/ -I/usr/local/opt/openblas/include/ -lopenblas -ffree-line-length-1024 SingleDoublePrecisionMatrixMultiplication.f90

我想现在我可以在我的代码中使用 MatrixVectorMultiplication，并且在以后的开发中，MatrixVectorMultiplication 将独立于它对 sgemv 或 dgemv 的依赖关系。

但是编译器会引发以下错误：

Error: Explicit interface required for ‘sgemv’ at (1): value argument

我将不胜感激任何可以帮助我找到解决这个问题的人。

Fortran Function-Pointers OpenBlas

评论

0赞 Ian Bush 11/17/2023
首先，你为什么认为你需要这个属性？标准的 blas 接口被设计为从 Fortran 77 调用，而 Fortran 77 早于它valuevalue
0赞 Ian Bush 11/17/2023
其次，为什么这么复杂？简单的重载将是我的选择
0赞 Ian Bush 11/17/2023
第三，不要像在这里对整数那样对种类使用显式常量。它是不可移植的，不能保证得到支持，并且可能无法执行您认为的功能。此外，BLAS 接口显示默认整数，那么为什么不使用默认整数呢？
0赞 Navid Khairdast 11/17/2023
谢谢@IanBush，您的解决方案运行良好。实际上，我对这个解决方案并不强硬。非常感谢！

答：

2赞 Ian Bush 11/17/2023 #1

正如评论中所指出的

  1. 你不需要.blas 接口被设计为 Fortran 77 可调用，因此不使用 valuevalue
  2. 我根本不会这样做。我会使用重载。然后，您需要更改的只是选择精度的种类参数
  3. 不要对种类值使用文本常量 - 始终使用由适当的查询例程设置的符号常量（也称为参数）。种类的文字常量本质上是不可移植的，不能保证工作，即使它确实工作，也可能无法达到您的预期。此外，BLAS 接口是使用默认整数类型定义的，那么为什么不使用默认整数类型呢？

这是我的做法。如果你仔细看一下代码的两个版本，你会发现所有的变化都是将工作精度类型设置为单精度 （ ） 或双精度 （wpspdp )

ijb@ijb-Latitude-5410:~/work/stack$ cat blas.f90
Module blas_overload_module

  Implicit None

  integer, parameter :: sp = selected_real_kind(  6, 37  )
  integer, parameter :: dp = selected_real_kind( 15, 307 )

  Interface gemv

     subroutine sgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy) 
       import                                  :: sp
       character                               :: trans
       integer                                 :: m, n, lda, incx, incy
       real(sp)                                :: alpha, beta
       real(sp), dimension(lda, *), intent(in) :: a
       real(sp), dimension(*), intent(in)      :: x
       real(sp), dimension(*), intent(inout)   :: y
     end subroutine sgemv

     subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy) 
       import                                  :: dp
       character                               :: trans
       integer                                 :: m, n, lda, incx, incy
       real(dp)                                :: alpha, beta
       real(dp), dimension(lda, *), intent(in) :: a
       real(dp), dimension(*), intent(in)      :: x
       real(dp), dimension(*), intent(inout)   :: y
     end subroutine dgemv

  End Interface gemv


End Module blas_overload_module

Program testit

  Use blas_overload_module, Only : sp, dp, gemv

  Implicit None

  Integer, Parameter :: n = 3

  Real( sp ), Dimension( 1:n, 1:n ) :: as
  Real( sp ), Dimension(      1:n ) :: xs
  Real( sp ), Dimension(      1:n ) :: ys_mm, ys_blas

  Real( dp ), Dimension( 1:n, 1:n ) :: ad
  Real( dp ), Dimension(      1:n ) :: xd
  Real( dp ), Dimension(      1:n ) :: yd_mm, yd_blas

  ! Change this line to set the precision - best put in a module
  Integer, Parameter :: wp = sp

  Real( wp ), Dimension( 1:n, 1:n ) :: aw
  Real( wp ), Dimension(      1:n ) :: xw
  Real( wp ), Dimension(      1:n ) :: yw_mm, yw_blas

  Call Random_number( as )
  Call Random_number( xs )
  Call gemv( 'N', n, n, 1.0_sp, as, n, xs, 1, 0.0_sp, ys_blas, 1 )
  ys_mm = Matmul( as, xs )
  Write( *, * ) 'Single: ', ys_mm - ys_blas

  Call Random_number( ad )
  Call Random_number( xd )
  Call gemv( 'N', n, n, 1.0_dp, ad, n, xd, 1, 0.0_dp, yd_blas, 1 )
  yd_mm = Matmul( ad, xd )
  Write( *, * ) 'Double: ', yd_mm - yd_blas

  Call Random_number( aw )
  Call Random_number( xw )
  Call gemv( 'N', n, n, 1.0_wp, aw, n, xw, 1, 0.0_wp, yw_blas, 1 )
  yw_mm = Matmul( aw, xw )
  Write( *, * ) 'Working: ', yw_mm - yw_blas

End Program testit
ijb@ijb-Latitude-5410:~/work/stack$ gfortran -Wall -Wextra -fcheck=all -std=f2018 -O -g blas.f90 -lopenblas
ijb@ijb-Latitude-5410:~/work/stack$ ./a.out
 Single:    0.00000000       0.00000000       0.00000000    
 Double:    0.0000000000000000        0.0000000000000000        0.0000000000000000     
 Working:    0.00000000       0.00000000       0.00000000    
ijb@ijb-Latitude-5410:~/work/stack$ gfortran -Wall -Wextra -fcheck=all -std=f2018 -O -g blas.f90 -lopenblas
ijb@ijb-Latitude-5410:~/work/stack$ cat blas.f90
Module blas_overload_module

  Implicit None

  integer, parameter :: sp = selected_real_kind(  6, 37  )
  integer, parameter :: dp = selected_real_kind( 15, 307 )

  Interface gemv

     subroutine sgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy) 
       import                                  :: sp
       character                               :: trans
       integer                                 :: m, n, lda, incx, incy
       real(sp)                                :: alpha, beta
       real(sp), dimension(lda, *), intent(in) :: a
       real(sp), dimension(*), intent(in)      :: x
       real(sp), dimension(*), intent(inout)   :: y
     end subroutine sgemv

     subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy) 
       import                                  :: dp
       character                               :: trans
       integer                                 :: m, n, lda, incx, incy
       real(dp)                                :: alpha, beta
       real(dp), dimension(lda, *), intent(in) :: a
       real(dp), dimension(*), intent(in)      :: x
       real(dp), dimension(*), intent(inout)   :: y
     end subroutine dgemv

  End Interface gemv


End Module blas_overload_module

Program testit

  Use blas_overload_module, Only : sp, dp, gemv

  Implicit None

  Integer, Parameter :: n = 3

  Real( sp ), Dimension( 1:n, 1:n ) :: as
  Real( sp ), Dimension(      1:n ) :: xs
  Real( sp ), Dimension(      1:n ) :: ys_mm, ys_blas

  Real( dp ), Dimension( 1:n, 1:n ) :: ad
  Real( dp ), Dimension(      1:n ) :: xd
  Real( dp ), Dimension(      1:n ) :: yd_mm, yd_blas

  ! Change this line to set the precision - best put in a module
  Integer, Parameter :: wp = dp

  Real( wp ), Dimension( 1:n, 1:n ) :: aw
  Real( wp ), Dimension(      1:n ) :: xw
  Real( wp ), Dimension(      1:n ) :: yw_mm, yw_blas

  Call Random_number( as )
  Call Random_number( xs )
  Call gemv( 'N', n, n, 1.0_sp, as, n, xs, 1, 0.0_sp, ys_blas, 1 )
  ys_mm = Matmul( as, xs )
  Write( *, * ) 'Single: ', ys_mm - ys_blas

  Call Random_number( ad )
  Call Random_number( xd )
  Call gemv( 'N', n, n, 1.0_dp, ad, n, xd, 1, 0.0_dp, yd_blas, 1 )
  yd_mm = Matmul( ad, xd )
  Write( *, * ) 'Double: ', yd_mm - yd_blas

  Call Random_number( aw )
  Call Random_number( xw )
  Call gemv( 'N', n, n, 1.0_wp, aw, n, xw, 1, 0.0_wp, yw_blas, 1 )
  yw_mm = Matmul( aw, xw )
  Write( *, * ) 'Working: ', yw_mm - yw_blas

End Program testit
ijb@ijb-Latitude-5410:~/work/stack$ gfortran -Wall -Wextra -fcheck=all -std=f2018 -O -g blas.f90 -lopenblas
ijb@ijb-Latitude-5410:~/work/stack$ ./a.out
 Single:    0.00000000       0.00000000       0.00000000    
 Double:    0.0000000000000000        0.0000000000000000        0.0000000000000000     
 Working:    0.0000000000000000        0.0000000000000000        0.0000000000000000     
ijb@ijb-Latitude-5410:~/work/stack$

评论

0赞 Navid Khairdast 11/20/2023
这篇文章中提出了另一个安静类似的解决方案：fortran-lang.discourse.group/t/pointer-to-blas-functions/6806/...
0赞 Ian Bush 11/20/2023
好吧，但我不喜欢激烈的预处理，我喜欢查看和编辑实际使用的内容。

上一个：为什么需要PLT存根开头的“jmp”？

下一个：迭代循环不一致，完全迭代