Incompatibilité de type de données dans fortran

Question

J'ai écrit un algorithme rudimentaire dans Fortran 95 pour calculer le gradient d'une fonction (dont un exemple est prescrit dans le code) en utilisant des différences centrales augmentées d'une procédure connue sous le nom d'extrapolation de Richardson.

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 

En simple précision, il fonctionne bien et me donne des résultats décents. Mais lorsque je tente de changer de double précision comme indiqué dans le code, je reçois une erreur lors de la tentative de compiler affirmant que la déclaration d'affectation

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

produit une incompatibilité de type real(4)/real(8). J'ai essayé plusieurs déclarations différentes de double précision, ajouté chaque constante double précision appropriée dans le code avec d0, et j'obtiens la même erreur à chaque fois. Je suis un peu perplexe quant à la façon dont la fonction f est susceptible de produire un nombre de précision unique.

Answer 1

Le problème est que f n'est pas explicitement défini dans votre programme principal, par conséquent il est implicitement supposé être de simple précision, qui est le type real (4) pour gfortran.

Je suis entièrement d'accord avec le commentaire de High Performance Mark, que vous devriez vraiment utiliser implicit none dans tout votre code fortran, pour vous assurer que tous les objets sont déclarés explicitement. De cette façon, vous auriez obtenu un message d'erreur plus approprié à propos de f n'étant pas explicitement défini.

En outre, vous pourriez envisager deux choses:

• Définir votre fonction dans un module et d'importation que le module dans le programme principal. Il est recommandé de définir tous les sous-programmes/fonctions dans les modules uniquement, afin que le compilateur puisse effectuer des vérifications supplémentaires sur le nombre et le type des arguments lorsque vous appelez la fonction.

• Vous pourriez (encore dans le module) introduire une constante pour la précision et l'utiliser partout, où le type d'un réel doit être spécifié. Prenant l'exemple ci-dessous, en changeant seulement la ligne

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

  dans

  integer, parameter :: dp = kind(1.0) 
  

  vous changeriez toutes vos variables réelles de double en simple précision. Notez également le suffixe _dp pour les constantes littérales au lieu du suffixe d0, qui ajusterait également automatiquement leur précision.

  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