Segmentation fault using F2003 on simple code

Question

I have written the following code in F2003, trying to implement an adaptive trapezoid method and adaptive simpson's rule method for simple integrals and I am having seg fault problems when I try to run the code using intrinsic functions (exp, log, sin) but it works perfectly for polynomials (x ** 2, for instance).

I have compiled like this:

$ gfortran -Wall -pedantic -fbounds-check integral_module.f03 integrate.f03 -o integral

And I get no warnings.

Any ideas??

integral_module.f03

module integral_module
 implicit none
 public :: integral
 !To test which is more efficient
 integer, public :: counter_int = 0, counter_simpson = 0
contains

recursive function integral(f, a, b, tolerance) &
           result(integral_result)

 intrinsic :: abs

 interface
     function f(x) result(f_result)
         real, intent(in) :: x
         real :: f_result
     end function f
 end interface

 real, intent(in) :: a, b, tolerance
 real :: integral_result
 real :: h, mid
 real :: trapezoid_1, trapezoid_2
 real :: left_area, right_area

 counter_int = counter_int + 1

 h = b - a
 mid = (a + b) / 2

 trapezoid_1 = h * (f(a) + f(b)) / 2.0
 trapezoid_2 = h / 2 * (f(a) + f(mid)) / 2.0 + &
               h / 2 * (f(mid) + f(b)) / 2.0

 if(abs(trapezoid_1 - trapezoid_2) < 3.0 * tolerance) then
     integral_result = trapezoid_2
 else
     left_area = integral(f, a, mid, tolerance / 2)
     right_area = integral(f, mid, b, tolerance / 2)
     integral_result = left_area + right_area
 end if

end function integral

recursive function simpson(f, a, h, tolerance) &
          result(simpson_result)

 intrinsic :: abs

 interface
     function f(x) result(f_result)
         real, intent(in) :: x
         real :: f_result
     end function f
 end interface
real, intent(in) :: a, h, tolerance
 real :: simpson_result
 real :: h_half, a_lower, a_upper
 real :: parabola_1, parabola_2
 real :: left_area, right_area

 counter_simpson = counter_simpson + 1

 h_half = h / 2.0
 a_lower = a - h_half
 a_upper = a + h_half

 parabola_1 = h / 3.0 * (f(a - h) + 4.0 * f(a) + f(a + h))

 parabola_2 = h_half / 3.0 * (f(a_lower - h_half) + 4.0 * f(a_lower) + f(a_lower + h_half)) + &
              h_half / 3.0 * (f(a_upper - h_half) + 4.0 * f(a_upper) + f(a_upper + h_half))

 if(abs(parabola_1 - parabola_2) < 15.0 * tolerance) then
     simpson_result = parabola_2
 else
     left_area = simpson(f, a_lower, h_half, tolerance / 2.0)
     right_area = simpson(f, a_upper, h_half, tolerance / 2.0)
     simpson_result = left_area + right_area
 end if

end function simpson

end module integral_module

And, integrate.f03

module function_module
 implicit none
 public :: non_para_errfun, parabola

contains
 function non_para_errfun(x) result(errfun_result)
     real, intent(in) :: x
     real :: errfun_result
     intrinsic :: exp
     errfun_result = exp(x ** 2.0)
 end function non_para_errfun

 function parabola(x) result(parabola_result)
     real, intent(in) :: x
     real :: parabola_result
     parabola_result = x ** 2.0
 end function parabola

 function brainerd(x) result(brainerd_result)
     real, intent(in) :: x
     real :: brainerd_result
     intrinsic :: exp, sin
     brainerd_result = exp(x ** 2.0) + sin(2.0 * x)
 end function brainerd
end module function_module

module math_module
 implicit none
 intrinsic :: atan
 real, parameter, public :: pi = &
       4.0 * atan(1.0)
end module math_module
program integrate

 use function_module,  f => brainerd
 use integral_module
 use math_module, only : pi
 implicit none
 intrinsic :: sqrt, abs

 real :: x_min, x_max, a, h, ans, ans_simpson
 real, parameter :: tolerance = 0.01

 x_min = 0
 x_max = 2.0 * pi
 a = abs(x_max) - abs(x_min)
 h = (abs(x_max) + abs(x_min)) / 2.0
 !ans = integral(f, x_min, x_max, tolerance)
 ans_simpson = simpson(f, a, h, tolerance)

 !print "(a, f11.6)", &
 !      "The integral is approx. ", ans
 print "(a, f11.6)", &
       "The Simpson Integral is approx: ", ans_simpson
 !print "(a, f11.6)", &
 !      "The exact answer is     ", & 
!      1.0/3.0 * x_max ** 3 - 1.0/3.0 * x_min ** 3

 print *

 !print "(a, i5, a)", "The trapezoid rule was executed ", counter_int, " times."
 print "(a, i5, a)", "The Simpson's rule was executed ", counter_simpson, " times."
end program integrate

Answer 1

You are getting a floating point exception when executing exp(x ** 2.0), since round about x > 9.35 the exponent gets too large for single precision floats.

Switching to double precision gets rid of the overflow.

You'd also want to change the x ** 2.0 to x**2, s.t. the compiler can optimize the integer power instead of doing the much more costly floating point operation.