Segmentation fault when passing a function as argument in a subroutine

Question

I try to illustrate how to pass a function to a Newton Raphson procedure. I succeed with a very simple function (called unefonction see below) but it does not work with a function which has got parameters. This second fonction is called gaussienne and it takes one argument, x, and two optional arguments mu and sig. In my newton raphson procedure I called the fonction in this way : f(x). What is strange for me is that during the execution, the program does as if the optional parameters sig and mu were present but they don't... Thus I do not understand ...

Here is the module which contains the functions

module fonction

  implicit none

  ! parametre pour la gaussienne
  double precision :: f_sigma = 1.d0, f_mu = 0.d0

  ! pi accessible uniquement en interne
  double precision, parameter :: pi = 3.14159265359d0

  contains

    double precision function unefonction(x)
        ! fonction : unefonction
        ! renvoie
        !    $\frac{e^x - 10}{x + 2}$

        implicit none

        ! arguments 
        double precision, intent(in) :: x

        unefonction = (exp(x) - 10.) / (x + 2.)

    end function unefonction

! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 

    double precision function gaussienne(x, mu, sig)
        ! fonction gaussienne
        ! utilise les parametres definis dans le module si
        ! mu et sig ne sont pas passes en argument

        implicit none

        ! arguments
        double precision, intent(in)           :: x
        double precision, intent(in), optional :: mu, sig

        ! variables locales
        double precision :: norme, moy, sigma

        ! sigma
        if (present(sig)) then
            write(*,*)"sig present"
            sigma = sig
        else
            sigma = f_sigma
        end if

        ! mu
        if (present(mu)) then
            write(*,*)"mu present"
            moy = mu
        else
            moy = f_mu
        end if

        ! calcul de la gaussienne
        norme = 1.d0 / (sigma * sqrt(2.d0 * pi))
        gaussienne = norme * exp(-(x - moy)**2 / (2.d0 * sigma**2))

    end function gaussienne

end module fonction

Here is the module which contains the newton raphson procedure

module rechercheRacine

    implicit none

    contains

        subroutine newtonRaphson(racine, f, eps, cible)

            ! recherche l'antecedant de cible

            implicit none

            ! arguments
            double precision, intent(inout)        :: racine
            double precision, intent(in), optional :: cible, eps

            ! fonction dont on cherche la racine
            double precision, external :: f

            ! variables locales
            integer          :: compteur
            double precision :: xold, xnew, delta, valcible
            double precision :: threshold, fprim, fdex

            ! precision
            if (present(eps)) then
                threshold = eps
            else
                threshold = 1.d-10
            end if

            ! valeur cible
            if (present(cible)) then
                valcible = cible
            else
                valcible = 0.d0
            end if

            write(*,*) "--------------------------------------------------------"
            write(*,*) "                      NEWTON RAPHSON"
            write(*,*) "--------------------------------------------------------"
            write(*,"('x0    = ',e16.6)") racine
            write(*,"('seuil = ',e16.6)") threshold
            write(*,"('cible = ',e16.6)") valcible
            write(*,*) "--------------------------------------------------------"
            write(*,*) "                        ITERATIONS"
            write(*,*) "--------------------------------------------------------"

            ! initialisation
            compteur = 0
            delta    = 1.d0
            xold     = racine
            write(*, '(i4,4e16.6)') compteur, f(xold), xold, 0., threshold

            ! iterations
            do while (delta > threshold .and. compteur <= 100)

                ! calcul de la fonction en xold
                fdex = f(xold) - valcible

                ! calcul de la derivee numerique
                fprim  = (f(xold + threshold) - f(xold - threshold)) / (2.d0 * threshold) 

                ! application de l'iteration de Newton Raphson
                xnew = xold - fdex / fprim
                delta = abs(xnew - xold)
                compteur = compteur + 1

                ! affichage de la convergence
                write(*, '(i4,4e16.6)') compteur, fdex, xnew, delta, threshold

                ! mise a jour de xstart
                xold = xnew
            end do

            if (delta < threshold) then
                racine = xnew
                write(*, *) '--------------------------------------------------------'
                write(*, *) '                        CONVERGE'
                write(*, *) '--------------------------------------------------------'
                write(*, *) 'A la convergence demandee, une solution est:'
                write(*, "('x = ',e20.10,'    f(x) = ', e20.10)") racine, f(racine)
                write(*, *)
            else
                write(*, *) '--------------------------------------------------------'
                write(*, *) '                      NON CONVERGE'
                write(*, *) '--------------------------------------------------------'
            end if

        end subroutine newtonRaphson

end module rechercheRacine

Here is the main program :

program main

    ! contient la subroutine newtonRaphson
    use rechercheRacine

    ! contient la fonction
    use fonction

    implicit none

    double precision :: racine, eps, cible

    ! appel de la subroutine newtonRaphson
    ! sans la valeur cible : cible (defaut = 0)
    ! sans la precision : eps (defaut 1d-10)
    racine = 1.d0
    call newtonRaphson(racine, unefonction)

    ! --------------------------------------------------------

    ! appel de la subroutine newtonRaphson
    ! avec pour cible 10
    racine = 1.d0
    eps    = 1.d-14
    cible  = 10.d0
    call newtonRaphson(racine, unefonction, eps, cible)

    ! --------------------------------------------------------

    ! parametre de la gaussienne
    f_sigma = 2.d0
    f_mu = 5.d0
    ! appel de la subroutine newtonRaphson
    ! passage des arguments sous la forme clef = valeur
    cible = 0.1d0
    racine = 2.d0
    call newtonRaphson(cible = cible, f = gaussienne, racine = racine)

end program main

The main program works for the function called unefonction but it doesn't work for the gaussienne function.

Here is the error message :

Program received signal SIGSEGV: Segmentation fault - invalid memory reference.

Backtrace for this error:
#0  0x7F1B6F5890F7
#1  0x7F1B6F5896D4
#2  0x7F1B6EEEB49F
#3  0x4009D2 in __fonction_MOD_gaussienne at mod_fonction.f90:54
#4  0x40104D in __rechercheracine_MOD_newtonraphson at mod_racine.f90:59
#5  0x4016BA in MAIN__ at main.f90:40
Erreur de segmentation (core dumped)

I think that the invalid memory reference is due to the fact that the program does as if optional parameters sig and mu were present and thus looks for them while they are not.

Answer 1

Yes, the problem is that the function you pass indeed expects three arguments instead of one only. If you change the external declaration of f in the subroutine newtonRaphson

double precision, external :: f

to an explicit interface (which describes, how you really use it):

interface
  double precision function f(x)
    double precision, intent(in) :: x
  end function f
end interface

your code won't even compile due to the mismatch in the number of parameters.

They are different ways to pass "parameters" to the function f which is called from the routine newtonRaphson:

• You could expect f to have two intent(in) arguments instead of one: Additional to the value x, it could take also a real array, which may be of arbitrary size and may contain the necessary parameters. That would require following interface:

  interface
    double precision function f(x, params)
      double precision, intent(in) :: x
      double precision, intent(in) :: params(:)
    end function f
  end interface
  

  Those functions, which do not need parameters (like unefonction) would just not use the content of the second parameter, while others (like gaussienne) would take their parameters from it.

• You could make newtonRaphson to expect a given extensible type (class) with a type bound procedure returning the value for a given x-value. You can then create abritrary extensions of this type, which may calculate the value for the given x-value based on some parameters stored as fields in the extended type. Then the program could look like below (I stripped several parts), but it would require Fortran 2003 compiler:

  module rechercheRacine
    implicit none
  
    type, abstract :: calculator
    contains
      procedure(getvalue_iface), deferred :: getvalue
    end type calculator
  
    interface
      double precision function getvalue_iface(self, x)
        import :: calculator
        class(calculator), intent(in) :: self
        double precision, intent(in) :: x
      end function getvalue_iface
    end interface
  
  contains
  
    subroutine newtonRaphson(racine, f, eps, cible)
      double precision, intent(inout)        :: racine
      class(calculator), intent(in)          :: f
      double precision, intent(in), optional :: cible, eps
  
      do while (delta > threshold .and. compteur <= 100)
        fdex = f%getvalue(xold) - valcible
        :
      end do
  
    end subroutine newtonRaphson
  
  end module rechercheRacine
  
  
  module fonction
    use rechercheRacine
    implicit none
  
    type, extends(calculator) :: unefonction
    contains
      procedure :: getvalue => unefonction_getvalue
    end type unefonction
  
    type, extends(calculator) :: gaussienne
      double precision :: mu = 0.0d0, sigma = 1.0d0
    contains
      procedure :: getvalue => gaussienne_getvalue
    end type gaussienne
  
  contains
  
    double precision function unefonction_getvalue(self, x)
      class(unefonction), intent(in) :: self
      double precision, intent(in) :: x
      unefonction_getvalue = (exp(x) - 10.) / (x + 2.)
    end function unefonction_getvalue
  
    double precision function gaussienne_getvalue(self, x)
      class(gaussienne), intent(in) :: self
      double precision, intent(in) :: x
  
      :
      gaussienne_getvalue = norme * exp(-(x - moy)**2 / (2.d0 * self%sigma**2))
  
    end function gaussienne_getvalue
  
  end module fonction
  
  
  program main
    use rechercheRacine
    use fonction
    implicit none
  
    type(unefonction) :: fone
    type(gaussienne) :: fgauss
  
    :
    call newtonRaphson(racine, fone)
    call newtonRaphson(cible = cible, f = fgauss, racine = racine)
  
  end program main