Fortran double precision program with a simple MKL BLAS routine

Question

In trying to mix precision in a simple program - using both real and double - and use the ddot routine from BLAS, I'm coming up with incorrect output for the double precision piece. Here's the code:

program test

!! adding this statement narrowed the issue down to ddot being considered real(4)
implicit none

integer, parameter :: dp = kind(1.0d0)

!! The following 2 lines were added for the calls to the BLAS routines.
!! This fixed the issue.
real(dp), external :: ddot
real, external :: sdot

real, dimension(3) :: a,b
real(dp), dimension(3) :: d,e

integer :: i

do i = 1,3
    a(i) = 1.0*i
    b(i) = 3.5*i
    d(i) = 1.0d0*i
    e(i) = 3.5d0*i
end do

write (*,200) "sdot real(4) = ", sdot(3,a,1,b,1)  ! should work and return 49.0
write (*,200) "ddot real(4) = ", ddot(3,a,1,b,1)  ! should not work

write (*,200) "sdot real(8) = ", sdot(3,d,1,e,1)  ! should not work
write (*,200) "ddot real(8) = ", ddot(3,d,1,e,1)  ! should work and return 49.0

200 format(a,f5.2)

end program test

I've tried compiling with both gfortran and ifort using the MKL BLAS libraries as follows:

ifort -lmkl_intel_lp64 -lmkl_sequential -lmkl_core

gfortran -lmkl_intel_lp64 -lmkl_sequential -lmkl_core main.f90

The output is:

sdot real(4) = 49.00
ddot real(4) =  0.00
sdot real(8) =  4.10
ddot real(8) =  0.00

How can I get the ddot routine to correctly process the double precision values?

Additionally, adding the -autodouble flag (ifort) or -fdefault-real-8 (gfortran) flag makes both of the ddot routines work, but the sdot routines fail.

Edit: I added the implicit none statement, and the two type statements for the ddot and sdot functions. Without the type specified for the function calls the ddot was being typed implicitly as single precision real.

Answer 1

I haven't used MKL, but perhaps you need a "use" statement so that the compiler knows the interface to the functions? Or to otherwise declare the functions. They are not declared so the compiler is probably assuming that return of ddot is single precision and mis-interpreting the bits.

Turning on the warning option causes the compiler to tell you about the problem. With gfortran, try: -fimplicit-none -Wall -Wline-truncation -Wcharacter-truncation -Wsurprising -Waliasing -Wimplicit-interface -Wunused-parameter -fwhole-file -fcheck=all -std=f2008 -pedantic -fbacktrace

Answer 2

Passing incorrect kind variables is a case of interface mismatch (which is illegal, so in principle the compiler might do anything including starting WW III), so maybe this is messing up the stack and hence following calls also return incorrect results. Try to comment out those incorrect calls (your lines marked with "should not work") and see if that helps.

Also, enable all kinds of debug options you can find, as e.g. the answer by M.S.B. shows for gfortran.