Pointer to BLAS functions

5

@certik’s solution works fine: at the end the following code coule be compiled in either single precision (-Dsingle) or double precision (-Ddouble) and MatrixVectorMultiplication() will do both.

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

   
    interface

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

    end subroutine dgemv

   end interface 

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

end module DoublePrecision

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

    interface

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

    end subroutine sgemv

   end interface 

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

end module SinglePrecision

program test_MatrixMultiplication


#ifdef single
   use SinglePrecision, only : precision, MatrixVectorMultiplication => sgemv
#else
   use DoublePrecision, only : precision, MatrixVectorMultiplication => dgemv
#endif

  Integer, Parameter :: n = 3
  Real( precision ), Dimension( 1:n, 1:n ) :: a
  Real( precision ), Dimension(      1:n ) :: x
  Real( precision ), Dimension(      1:n ) :: y_blas
  Call Random_number( a )
  Call Random_number( x )
  Call MatrixVectorMultiplication( 'N', n, n, 1.0_precision, a, n, x, 1, 0.0_precision, y_blas, 1 )

end program test_MatrixMultiplication

To compile the code, one needs to pass -ldl flag for dynamic linking:

gfortran -Dsingle -cpp -L/path/to/openblas/lib/ -I/path/to/openblasinclude/ -lopenblas -ldl  -ffree-line-length-1024 test_MatrixMultiplication.f90

Thanks for this solution

1 Like