How to properly call the SGEMV in Fortran?

Question

I want to perform a Matrix-Vector product in fortran using the SGEMV subroutine from BLAS. I have a code that is similar to this:

program test
integer, parameter :: DP = selected_real_kind(15)
real(kind=DP), dimension (3,3) :: A
real(kind=DP), dimension (3) :: XP,YP
call sgemv(A,XP,YP)

A is a 3x3 Matrix, XP and YP are Vectors. In the included module one can see the following code:

PURE SUBROUTINE SGEMV_F95(A,X,Y,ALPHA,BETA,TRANS)
    ! Fortran77 call:
    ! SGEMV(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
    USE F95_PRECISION, ONLY: WP => SP
    REAL(WP), INTENT(IN), OPTIONAL :: ALPHA
    REAL(WP), INTENT(IN), OPTIONAL :: BETA
    CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: TRANS
    REAL(WP), INTENT(IN) :: A(:,:)
    REAL(WP), INTENT(IN) :: X(:)
    REAL(WP), INTENT(INOUT) :: Y(:)
END SUBROUTINE SGEMV_F95

I understand that the some of the parameters are optional, so where am i wrong in the method call?

Answer 1

When you look at BLAS or LAPACK routines then you should always have a look at the first letter:

• S: single precision
• D: double precision
• C: single precision complex
• Z: double precision complex

You defined your matrix A as well as the vectors XP and YP as a double precision number using the statement:

integer, parameter :: DP = selected_real_kind(15)

So for this, you need to use dgemv or define your precision as single precision.

There is also a difference between calling dgemv and dgemv_f95. dgemv_f95 is part of Intel MKL and not really a common naming. For portability reasons, I would not use that notation but stick to the classic dgemv which is also part of Intel MKL.

DGEMV performs one of the matrix-vector operations

y := alpha*A*x + beta*y,   or   y := alpha*A**T*x + beta*y,

where alpha and beta are scalars, x and y are vectors and A is an m by n matrix.

If you want to know how to call the function, I suggest to have a look here, but it should, in the end, look something like this:

call DGEMV('N',3,3,ALPHA,A,3,XP,1,BETA,YP,1)