Search code examples
fortranfunction-pointersopenblas

Pointer to OpenBLAS subroutines in fortran

In my code, I want to set pointers to OpenBLAS subroutines, to compile the code either in single or double precisions.

To do so, I define two modules and I define function interfaces for single(sgemv) or double precision(dgemv) OpenBLAS function:


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

and

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

In the main program, I take a flag in the compile command, and I use the appropriate module:

program test_MatrixMultiplication


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


end program test_MatrixMultiplication

To compile, I use the following command:

! 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

I suppose that now I can use MatrixVectorMultiplication in my code, and the MatrixVectorMultiplication in the later development becomes independent of its dependecy to either sgemv, or dgemv.

But the compiler raise the following error:

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

I would appreciate anyone who can help me to find a solution to this problem.

Solution

  • As noted in the comments

    1. You don't need value. The blas interfaces are designed to be Fortran 77 callable, and thus do not use value
    2. I wouldn't do it this way at all. I would use overloading. Then all you need to change is the kind parameter to choose the precision
    3. Don't use literal constants for kind values - always use symbolic constants (aka parameters) set by an appropriate inquiry routine. Literal constants for kinds are inherently non-portable, not guaranteed to work, and even if it does work it might not do what you intend. Further the BLAS interfaces are defined using default integer kind, so why not use default integer kind?

    Here is how I would do it. If you look carefully at the two versions of the code you will see is all that has changed is setting the working precision kind wp to single ( sp ) or double ( 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$