Pointer to BLAS functions

1

Hello everyone,

I want to do matrix-vector multiplication in either single or double pricision. To do so, I defined two modules for each precision, and at each one I defined pointers to function, dgemv for double precision, and sgemv for single precision.

The code is as follows:

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

   
   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)
   
   

   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

To compile with the single precision, I use

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

and 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

The compiler raise the error below:

SingleDoublePrecisionMatrixMultiplication.f90:60:33:

   60 |    MatrixVectorMultiplication => sgemv
      |                                 1
Error: Explicit interface required for 'sgemv' at (1): value argument

Thank you, everyone, for your assistance in resolving this issue.

2

I think the correct syntax is use SinglePrecision, only: MatrixVectorMultiplication => sgemv.

1 Like
3

I agree with @certik that the only: clause is probably what you want to use. However, perhaps some additional explanation is indicated. In your original code you defined a procedure pointer within the module, and you then assign the procedure pointer at run time. With the only clause, there is no procedure pointer involved, the statement is simply a local alias, and everything is done at compile time. That means that you get the normal compile time argument checking when you reference the procedure.

2 Likes
4

This solution works fine, thanks a alot!

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
6

@Kheirdast awesome, I am glad you got it working!