Fortran 中 OpenBLAS 子例程的指针

Question

在我的代码中，我想设置指向 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

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

Answer 1

如评论中所述

1. 你不需要

   value

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

我会这样做。如果仔细查看代码的两个版本，您会发现所有更改都是将工作精度类型

wp

设置为单精度 (

sp

) 或双精度 (

dp

)

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$