Data type mismatch in fortran

I've written a rudimentary algorithm in Fortran 95 to calculate the gradient of a function (an example of which is prescribed in the code) using central differences augmented with a procedure known as Richardson extrapolation.

function f(n,x)
! The scalar multivariable function to be differentiated

integer :: n
real(kind = kind(1d0)) :: x(n), f

f = x(1)**5.d0 + cos(x(2)) + log(x(3)) - sqrt(x(4))

end function f
!=====!
!=====!
!=====!

program gradient
!==============================================================================!
! Calculates the gradient of the scalar function f at x=0using a finite        !
! difference approximation, with a low order Richardson extrapolation.         !
!==============================================================================!

parameter (n = 4, M = 25)
real(kind = kind(1d0)) :: x(n), xhup(n), xhdown(n), d(M), r(M), dfdxi, h0, h, gradf(n)

h0 = 1.d0
x  = 3.d0

! Loop through each component of the vector x and calculate the appropriate
! derivative
do i = 1,n
! Reset step size
h = h0

    ! Carry out M successive central difference approximations of the derivative
do j = 1,M
        xhup = x
        xhdown = x
        xhup(i) = xhup(i) + h
        xhdown(i) = xhdown(i) - h
        d(j) = ( f(n,xhup) - f(n,xhdown) ) / (2.d0*h)
    h = h / 2.d0
    end do

    r = 0.d0
    do k = 3,M      r(k) = ( 64.d0*d(k) - 20.d0*d(k-1) + d(k-2) ) / 45.d0
        if ( abs(r(k) - r(k-1)) < 0.0001d0 ) then
        dfdxi = r(k)
            exit
        end if
    end do

    gradf(i) = dfdxi
end do

! Print out the gradient
write(*,*) " "
write(*,*) " Grad(f(x)) = "
write(*,*) " "
do i = 1,n
    write(*,*) gradf(i)
end do

end program gradient

In single precision it runs fine and gives me decent results. But when I try to change to double precision as shown in the code, I get an error when trying to compile claiming that the assignment statement

d(j) = ( f(n,xhup) - f(n,xhdown) ) / (2.d0*h)

is producing a type mismatch real(4)/real(8). I have tried several different declarations of double precision, appended every appropriate double precision constant in the code with d0, and I get the same error every time. I'm a little stumped as to how the function f is possibly producing a single precision number.

The problem is that f is not explicitely defined in your main program, therefore it is implicitly assumed to be of single precision, which is the type real(4) for gfortran.

I completely agree to the comment of High Performance Mark, that you really should use implicit none in all your fortran code, to make sure all object are explicitely declared. This way, you would have obtained a more appropriate error message about f not being explicitely defined.

Also, you could consider two more things:

Define your function within a module and import that module in the main program. It is a good practice to define all subroutines/functions within modules only, so that the compiler can make extra checks on number and type of the arguments, when you invoke the function.

You could (again in module) introduce a constant for the precicision and use it everywhere, where the kind of a real must be specified. Taking the example below, by changing only the line

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

into

integer, parameter :: dp = kind(1.0)

you would change all your real variables from double to single precision. Also note the _dp suffix for the literal constants instead of the d0 suffix, which would automatically adjust their precision as well.

module accuracy
  implicit none

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

end module accuracy


module myfunc
  use accuracy
  implicit none

contains

  function f(n,x)
    integer :: n
    real(dp) :: x(n), f
    f = 0.5_dp * x(1)**5 + cos(x(2)) + log(x(3)) - sqrt(x(4))
  end function f

end module myfunc


program gradient
  use myfunc
  implicit none

  real(dp) :: x(n), xhup(n), xhdown(n), d(M), r(M), dfdxi, h0, h, gradf(n)
  :

end program gradient