In Fortran 95 habe ich einen rudimentären Algorithmus geschrieben, um den Gradienten einer Funktion (für die im Code ein Beispiel vorgeschrieben ist) unter Verwendung zentraler Differenzen zu berechnen, die durch eine als Richardson-Extrapolation bekannte Prozedur ergänzt werden.

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 einfacher Präzision läuft es gut und gibt mir gute Ergebnisse. Aber wenn ich versuche, auf doppelte Genauigkeit zu wechseln, wie im Code gezeigt, erhalte ich eine Fehlermeldung, wenn ich versuche, die Zuweisungsanweisung zu kompilieren

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

erzeugt eine Typinkongruenzreal(4)/real(8). Ich habe verschiedene Deklarationen mit doppelter Genauigkeit ausprobiert und jede geeignete Konstante mit doppelter Genauigkeit im Code mit angehängtd0, und ich erhalte jedes Mal den gleichen Fehler. Ich bin ein wenig ratlos, wie das funktioniertf produziert möglicherweise eine einzelne Genauigkeitszahl.