A simple algorithm for adaptive numerical integration of function in a finite interval can easily be programmed with a recursive function in Fortran90 or C. The original program is taken from the collection of examples coming with the F compiler.

The main routine sets the interval boundaries, calls the integral function with the quadrature algorithm and prints the numerical and the exact results, which is known for the given simple function.

main.f95main.c
program integrate

   use function_module
   use integral_module

   real, parameter :: pi = &
3.1415926535897932384626433832795028841972
   real :: x_min, x_max
   real :: answer

   x_min = -4.0
   x_max = 4.0
   answer = integral (func, x_min, x_max, 0.001)
   print "(a, f20.12)",  &
         "The integral is approximately ",  &
          answer
   print "(a, f20.12)",  &
         "The exact answer is           ",  &
          sqrt (pi)

end program integrate
#include <stdio.h>

double  integral(
    double (*f)(double), /* function to integrate */
    double a,            /* left interval boundary  */
    double b,            /* right interval boundary */
    double tolerance);   /* error tolerance */

double func(double x);

int main(int argc,char** argv)
{
  double x_min, x_max;
  double answer = 0.0;
  double tm;
  int i;

  x_min = -4.0;
  x_max = 4.0;

  answer = integral(func, x_min, x_max, 0.001);

  printf("The integral is approximately %.12f\n", answer);
  printf("The exact answer is           %.12f\n", sqrt(PI));

  return 0;
}

The function module evaluates the function at a given point.

function.f95function.c
module function_module

public :: func

contains

function func (x) result (f_result)

   real, intent (in) :: x
   real :: f_result

   f_result = exp (-x ** 2)

end function func

end module function_module
#include <math.h>

/* function to integrate */
double func(double x)
{
  return  exp(-x*x);
}

Finally the integral module contains the quadrature algorithm as a recursive function.

integral_module.f95integral_module.c
module integral_module
public :: integral

contains

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

   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 :: one_trapezoid_area, two_trapezoid_area
   real :: left_area, right_area

   h = b - a
   mid = (a + b) /2
   one_trapezoid_area = h * (f(a) + f(b)) / 2.0
   two_trapezoid_area = h/2 * (f(a) + f(mid)) / 2.0 + &
                        h/2 * (f(mid) + f(b)) / 2.0
   if (abs(one_trapezoid_area - two_trapezoid_area)  &
         < 3.0 * tolerance) then
      integral_result = two_trapezoid_area
   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

end module integral_module
double  integral(
                 double (*f)(),    /* function to integrate */
                 double a,         /* left interval boundary  */
                 double b,         /* right interval boundary */
                 double tolerance) /* error tolerance */
{
  double integral_result;       /* value to return */
  double h;                     /* interval size */
  double mid;                   /* center of interval */
  double one_trapezoid_area;    /* area of one trapezoid */
  double two_trapezoid_area;    /* area of two trapezoids */
  double left_area;             /* integral of left half interval */
  double right_area;            /*     "       right    "         */


  h = b - a;
  mid = (a+b)/2;
  one_trapezoid_area = h * (f(a) + f(b)) / 2.0;
  two_trapezoid_area = h/2 * (f(a) + f(mid)) / 2.0 +
    h/2 * (f(mid) + f(b)) / 2.0;


  if (fabs(one_trapezoid_area - two_trapezoid_area) 
            < 3.0 * tolerance){
    /* error acceptable   */
    integral_result = two_trapezoid_area;

  }else{
    /* error not acceptable */
    /* recursiv function calls for left and right areas */

      {
          left_area = integral(f, a, mid, tolerance / 2);
          right_area = integral(f, mid, b, tolerance / 2);
      }

      integral_result = left_area + right_area;
    }

  return integral_result;
}

The integral over the function f is approximated by two quadrature formulas: The trapezoidal rule (one_trapezoid_area) and the midpoint rule (two_trapezoid_area). The difference between these values is used to estimate the error. If the error is less than the given tolerance, the result is accepted and returned as the function result. Otherwise the interval is cut into halves as well as the demanded tolerance and the integral function is recursively called for both new intervals. The  OpenMP 3.0 standard offers this  feature, which is well suited to parallelize a recursive algorithm very elegantly.

main_task.c
#include <stdio.h>
#include <math.h>

#ifdef _OPENMP
#include <omp.h>
#endif

double  integral_par(
          double (*f)(double), /* function to integrate */
          double a,            /* left interval boundary  */
          double b,            /* right interval boundary */
          double tolerance);   /* error tolerance */

double  integral(
          double (*f)(double), /* function to integrate */
          double a,            /* left interval boundary  */
          double b,            /* right interval boundary */
          double tolerance)    /* error tolerance */
{
  double answer = 0.0;
  
#pragma omp parallel
  {
    answer = integral_par(f, a, b, tolerance);
  } /* omp parallel */
  return answer;
}
	
integral_par.c
#include <stdio.h>
#include <math.h>

#ifdef _OPENMP
#include <omp.h>
#endif

double integral_par( double (*f)(double), /* function to integrate */ double a, /* left interval boundary */ double b, /* right interval boundary */ double tolerance) /* error tolerance */ { double integral_result; /* value to return */ double h; /* interval size */ double mid; /* center of interval */ double one_trapezoid_area; /* area of one trapezoid */ double two_trapezoid_area; /* area of two trapezoids */ double left_area; /* integral of left half interval */ double right_area; /* " right " */ h = b - a; mid = (a+b)/2; one_trapezoid_area = h * (f(a) + f(b)) / 2.0; two_trapezoid_area = h/2 * (f(a) + f(mid)) / 2.0 + h/2 * (f(mid) + f(b)) / 2.0; if (fabs(one_trapezoid_area - two_trapezoid_area) < 3.0 * tolerance){ /* error acceptable */ integral_result = two_trapezoid_area; }else{ /* error not acceptable */ /* put recursiv function calls for left and right areas */ /* into task queue */ #pragma omp task shared(left_area) { left_area = integral_par(f, a, mid, tolerance / 2); } #pragma omp task shared(right_area) { right_area = integral_par(f, mid, b, tolerance / 2); }
#pragma omp taskwait integral_result = left_area + right_area; } return integral_result; }

 

Here we use the OpenMP Tasking model. Can this function be parallelized elegantly and efficiently anyway?
The recursive function call can be easily replaced by a stack mechanism, used for storing intervals and corresponding tolerances, and a loop.

stack.f95stack.h
module data_m
        type data_t
                real :: left, right, tol
                type(data_t), pointer :: pred
        end type data_t
end module data_m

module stack_m
        use data_m
        type stack_t
                integer :: count
                type (data_t), pointer :: stackpointer
        end type stack_t
contains
        subroutine new_stack ( stack )
                implicit none
                type (stack_t) :: stack
                nullify ( stack%stackpointer )
                stack%count = 0
        end subroutine new_stack

        function empty_stack ( stack ) result ( res )
                implicit none
                type (stack_t), intent(in) :: stack
                logical :: res
                res = stack%count == 0
        end function empty_stack

        subroutine push ( stack, left, right, tol )
                implicit none
                type (stack_t), intent(inout) :: stack
                real, intent(in) :: left, right, tol
                type (data_t), pointer :: new_data

                allocate (new_data)
                new_data%left = left
                new_data%right = right
                new_data%tol = tol

                if ( stack%count == 0 ) then ! empty stack so far
                        nullify (new_data%pred )
                        stack%stackpointer => new_data
                else ! append
                        new_data%pred => stack%stackpointer
                        stack%stackpointer => new_data
                end if
                stack%count = stack%count + 1
        end subroutine push

        subroutine pop ( stack, left, right, tol )
                implicit none
                type (stack_t), intent(inout) :: stack
                real, intent(out) :: left, right, tol
                type (data_t), pointer :: tmp

                if ( stack%count > 1 ) then
                        left = stack%stackpointer%left
                        right = stack%stackpointer%right 
                        tol = stack%stackpointer%tol
                        tmp => stack%stackpointer%pred
                        deallocate ( stack%stackpointer )
                        stack%stackpointer => tmp
                        stack%count = stack%count - 1
                else if ( stack%count == 1 ) then
                        left = stack%stackpointer%left
                        right = stack%stackpointer%right 
                        tol = stack%stackpointer%tol
                        deallocate ( stack%stackpointer )
                        nullify ( stack%stackpointer )
                        stack%count = 0
                endif
        end subroutine pop

end module stack_m 
stack.h:

/* the stack structure */
struct stack_s{               
  int el_count;            /* count of elements on stack */    
  int el_size;             /* size of an element */
  int mem_reserve;         /* allocated memory for stack */
  void* elements;          /* pointer to begin of stack */
};

typedef struct stack_s* stack_t;


void create_stack(stack_t* stack, int element_size);
int empty_stack(stack_t stack);
void push_stack(stack_t stack, void* element);
void pop_stack(stack_t stack, void* element);

stack.c
stack.c:

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "stack.h"

#define INITIAL_SIZE 128   /* initial size of new stacks */


/******************************************
 * create new stack
 ******************************************/
void 
create_stack(
             stack_t* stack,     /* stack to create */
             int element_size)   /* size of a stack element */
{
  int initial_size = INITIAL_SIZE;

  /* allocate memory for new stack struct */
  (*stack) = (stack_t) malloc(sizeof(struct stack_s));
  if (!(*stack)){
    fprintf(stderr, "error: could not allocate memory for stack.. Abort.\n"); 
    exit(1);
  }    

  /* allocate memory for stack elements */
  (*stack)->elements = (void*) malloc(element_size * initial_size);
  (*stack)->mem_reserve = initial_size; 
  if (!(*stack)->elements){
    fprintf(stderr, "error: could not allocate memory for stack.. Abort.\n");
    exit(1);
  }

  (*stack)->el_size = element_size;
  (*stack)->el_count = 0;

}

/*****************************************
 * check if the stack is empty 
 *****************************************/
int 
empty_stack
(stack_t stack)
{
  return stack->el_count <= 0;
}


/*****************************************
 * push a element on stack
 *****************************************/
void 
push_stack(
           stack_t stack,    /* target stack */
           void* element)    /* element to push */
{
  int i, new_reserve;
  int log2_count;

  /* check if we need more memory for stack */    
  if (stack->el_count >= stack->mem_reserve){

      /* calculate new size for the stack
         it should be a power of two */
      for (i = stack->el_count, log2_count = 0; 
           i > 0; 
           i>>1, log2_count++);
      new_reserve = 1 << log2_count;
 
      /* reallocate memory for phase thread tables 
         and nullify new values */
      stack->elements = (void *) realloc(stack->elements, 
                      stack->el_size * new_reserve);
      if (!stack->elements){
        fprintf(stderr, "error: can't reallocate stack.. Aborting\n");
        exit(1);
      }

      stack->mem_reserve = new_reserve;
  }
  
  /* now push the element on top of the stack */
  memcpy((char*)stack->elements + stack->el_count*stack->el_size, 
            element, stack->el_size);
  stack->el_count++;

}


/*****************************************
 * pop an element from stack
 *****************************************/
void 
pop_stack(
          stack_t stack,    /* target stack */
          void* element)    /* where poped el. should be stored */
{
  if (stack->el_count <= 0){
    fprintf(stderr, "error: trying to pop from empty stack.\n");
    exit(2);
  }

  stack->el_count--;
  memcpy(element, 
         (char*)stack->elements + stack->el_count*stack->el_size, 
         stack->el_size);
}




The initial interval is put on the stack. While there still are intervals on the stack, there is something to do. If the integration over the interval lying on top of the stack is successful, the stack is shrinking, if not, the two new intervals are put on the stack.
Once the stack module is available the new integral function is not much more complicated than the original version.

integral_module.f95integral_module.c
module integral_module
use stack_m
public :: integral

contains

function integral (f, ah, bh, tolh)  &
      result (integral_result)

   interface
   function f (x) result (f_result)
      real, intent (in) :: x
      real :: f_result
   end function f
   end interface
   real, intent (in) :: ah, bh, tolh
   real :: a, b, tolerance
   real :: integral_result
   real :: h, mid
   real :: one_trapezoid_area, two_trapezoid_area
   real :: left_area, right_area
   type (stack_t) :: stack

   call new_stack ( stack )
   call push ( stack, ah, bh, tolh )
   
   integral_result = 0.0
   do
      if ( empty_stack ( stack ) ) exit
      call pop ( stack, a, b, tolerance )
   
      h = b - a
      mid = (a + b) /2
      one_trapezoid_area = h * (f(a) + f(b)) / 2.0
      two_trapezoid_area = h/2 * (f(a) + f(mid)) / 2.0 + &
                           h/2 * (f(mid) + f(b)) / 2.0
      if (abs(one_trapezoid_area - two_trapezoid_area)  &
            < 3.0 * tolerance) then
         integral_result = integral_result + two_trapezoid_area
      else
         call push ( stack, a, mid, tolerance / 2 )
         call push ( stack, mid, b, tolerance / 2 )
      end if

   end do

end function integral

end module integral_module
 
#include "stack.h"

/* a part (element) of shared work */
typedef struct _work_t{
  double a;
  double b;
  double tol;
} work_t;

/**************************************/
 **************************************/
double  
integral(
         double (*f)(double), /* function to integrate */
         double a,            /* left interval boundary  */
         double b,            /* right interval boundary */
         double tol)          /* error tolerance */
{
  double integral_result;       /* value to return */
  double h;                     /* interval size */
  double mid;                   /* center of interval */
  double one_trapezoid_area;    /* area of one trapezoid */
  double two_trapezoid_area;    /* area of two trapezoids */
  double left_area;             /* integral of left half interval */
  double right_area;            /*     "       right    "         */
  double tmp;                   

  stack_t stack;
  work_t work;

                                /* prepare stack */
  work.a = a;
  work.b = b;
  work.tol = tol;

  create_stack(&stack, sizeof(work_t));
  push_stack(stack, &work);

  integral_result = 0.0;

  while(!empty_stack(stack)){

    pop_stack(stack, &work);
    b = work.b;
    a = work.a;
    tolerance = work.tol;

    h = b - a;
    mid = (a+b)/2;
    one_trapezoid_area = h * (f(a) + f(b)) / 2.0;
    two_trapezoid_area = h/2 * (f(a) + f(mid)) / 2.0 +
      h/2 * (f(mid) + f(b)) / 2.0;


    if (fabs(one_trapezoid_area - two_trapezoid_area) < 3.0 * tolerance){
            /* error acceptable   */
      integral_result += two_trapezoid_area;

    }else{
      /* error not acceptable */

      /* push (a,mid, tol/2) onto stack */
      tmp = work.b;
      work.b = mid;
      work.tol = work.tol/2;

      push_stack(stack, &work);

      /* push (mid, b, tol/2) onto stack */
      work.a = mid;
      work.b = tmp;

      push_stack(stack, &work);
    }

  } /* while */

  return integral_result;
}


The first parallelization of this algorithm is very easy: All stack accesses have to be put into critical regions and the summation of the integral results as well.

integral_module.f95
module integral_module
use stack_m
public :: integral

contains

function integral (f, ah, bh, tolh)  &
      result (integral_result)

   interface
   function f (x) result (f_result)
      real, intent (in) :: x
      real :: f_result
   end function f
   end interface
   real, intent (in) :: ah, bh, tolh
   real :: a, b, tolerance
   real :: integral_result
   real :: h, mid
   real :: one_trapezoid_area, two_trapezoid_area
   real :: left_area, right_area
   type (stack_t) :: stack
   logical :: ready

   call new_stack ( stack )
   call push ( stack, ah, bh, tolh )
   integral_result = 0.0

!$omp parallel default(none) &
!$omp shared(stack,integral_result,f) &
!$omp private(a,b,tolerance,h,mid,one_trapezoid_area,two_trapezoid_area,ready)
   ready = .false.
   
   do
!$omp critical (stack)  
      if ( empty_stack ( stack ) ) then
         ready=.true.
      else
         call pop ( stack, a, b, tolerance )
      end if
!$omp end critical (stack)    
      if ( ready ) exit
      
      h = b - a
      mid = (a + b) /2
      one_trapezoid_area = h * (f(a) + f(b)) / 2.0
      two_trapezoid_area = h/2 * (f(a) + f(mid)) / 2.0 + &
                           h/2 * (f(mid) + f(b)) / 2.0
      if (abs(one_trapezoid_area - two_trapezoid_area)  &
            < 3.0 * tolerance) then
!$omp critical (result)  
         integral_result = integral_result + two_trapezoid_area
!$omp end critical (result)  
      else
!$omp critical (stack)  
         call push ( stack, a, mid, tolerance / 2 )
         call push ( stack, mid, b, tolerance / 2 )
!$omp end critical (stack)    
      end if

   end do
!$omp end parallel

end function integral

end module integral_module

But there still is one open problem: If somewhere in the middle of the adaptive integration process, the stack has less entries than there are threads, the supernumerous threads will exit the loop and wait at the end of the parallel region. This might already happen in the very beginning.
Therefore a counter is introduced, counting the number of threads currently catively working on an interval. Only if this counter has been set to zero, a thread, which encounters an empty stack will exit the loop.

integral_module.f95integral_module.c
module integral_module
use stack_m
public :: integral

contains

function integral (f, ah, bh, tolh)  &
      result (integral_result)

   interface
   function f (x) result (f_result)
      real, intent (in) :: x
      real :: f_result
   end function f
   end interface
   real, intent (in) :: ah, bh, tolh
   real :: a, b, tolerance
   real :: integral_result
   real :: h, mid
   real :: one_trapezoid_area, two_trapezoid_area
   real :: left_area, right_area
   type (stack_t) :: stack
   logical :: idle, ready
   integer :: busy

   call new_stack ( stack )
   call push ( stack, ah, bh, tolh )
   integral_result = 0.0
   busy = 0
   ready = .false.

!$omp parallel default(none) &
!$omp shared(stack,integral_result,f,busy) &
!$omp private(a,b,tolerance,h,mid,one_trapezoid_area,&
!$omp         two_trapezoid_area,idle,ready)
   idle = .true.
   
   do
!$omp critical (stack)  
      if ( empty_stack ( stack ) ) then
         if ( .not. idle ) then
            idle=.true.
            busy = busy - 1
         end if
         if ( busy .eq. 0 ) ready = .true.
      else
         call pop ( stack, a, b, tolerance )
         if ( idle ) then
            idle = .false.
            busy = busy + 1 
         end if
      end if
!$omp end critical (stack)    
      if ( idle ) then
         if ( ready ) exit
         !call idle_loop ( 0.001 )
         cycle
      end if
      
      h = b - a
      mid = (a + b) /2
      one_trapezoid_area = h * (f(a) + f(b)) / 2.0
      two_trapezoid_area = h/2 * (f(a) + f(mid)) / 2.0 + &
                           h/2 * (f(mid) + f(b)) / 2.0
      if (abs(one_trapezoid_area - two_trapezoid_area)  &
            < 3.0 * tolerance) then
!$omp critical (result)  
         integral_result = integral_result + two_trapezoid_area
!$omp end critical (result)  
      else
!$omp critical (stack)  
         call push ( stack, a, mid, tolerance / 2 )
         call push ( stack, mid, b, tolerance / 2 )
!$omp end critical (stack)    
      end if

   end do
!$omp end parallel

end function integral

end module integral_module

#include "stack.h"

typedef struct _work_t{
  double a;
  double b;
  double tol;
} work_t;


double  
integral(
     double (*f)(double), /* function to integrate */
     double ah,           /* left interval boundary  */
     double bh,           /* right interval boundary */
     double tolh)         /* error tolerance */

{
  double a, b, tolerance;       /* vars poped from stack */
  double integral_result;       /* value to return */
  double h;                     /* interval size */
  double mid;                   /* center of interval */
  double one_trapezoid_area;    /* area of one trapezoid */
  double two_trapezoid_area;    /* area of two trapezoids */
  double left_area;             /* integral of left half interval */
  double right_area;            /*     "       right    "         */

  stack_t stack;
  work_t work;

  int ready, idle, busy;

  /* prepare stack */
  work.a = ah;
  work.b = bh;
  work.tol = tolh;

  create_stack(&stack, sizeof(work_t));
  push_stack(stack, &work);

  integral_result = 0.0;

  busy = 0;

#pragma omp parallel default(none) \
    shared(stack, integral_result,f,busy) \
    private(a,b,tolerance, work, h, mid, one_trapezoid_area, \
            two_trapezoid_area, idle, ready)
  {

    ready = 0;
    idle = 1;

    while(!ready){

#pragma omp critical (stack)
      {
        if (!empty_stack(stack)){
          /* we have new work */ 
          pop_stack(stack, &work);

          if (idle){
            /* say others i'm busy */
            busy += 1;
            idle = 0;
          }

        }else{
          /* no new work on stack */
          if (!idle){
            busy -= 1;
            idle = 1;
          }

          /* nobody has anything to do; let us leave the loop */
          if (busy == 0)
            ready = 1;

        }
      } /* end critical(stack) */

      if (idle)
        continue;

      b = work.b;
      a = work.a;
      tolerance = work.tol;

      h = b - a;
      mid = (a+b)/2;
      one_trapezoid_area = h * (f(a) + f(b)) / 2.0;
      two_trapezoid_area = h/2 * (f(a) + f(mid)) / 2.0 +
        h/2 * (f(mid) + f(b)) / 2.0;


      if (fabs(one_trapezoid_area - two_trapezoid_area) < 3.0 * tolerance){
        /* error acceptable   */

#pragma omp critical (result)
        integral_result += two_trapezoid_area;

      }else{
        /* error not acceptable */
        /* push new work into stack */

        work.a = a;
        work.b = mid;
        work.tol = tolerance/2;
#pragma omp critical (stack)
        {
          push_stack(stack, &work);
      
          work.a = mid;
          work.b = b;
          push_stack(stack, &work);

        }
      }

    } /* while */

  } /* end omp parallel */

  return integral_result;
}


The summation of the integral results can be further optimized, by using private variables for the partial sums, which each thread can accumulate without a critical region.
Only at the very end, these partial sums are then summed up in a critical region, which is only entered once by each thread.

...
!$omp ... private(...,private_result)
...
      if (abs(one_trapezoid_area - two_trapezoid_area)  &
            < 3.0 * tolerance) then
         private_result = private_result + two_trapezoid_area
      else
!$omp critical (stack)  
         call push ( stack, a, fa, mid, fmid, tolerance / 2 )
         call push ( stack, mid, fmid, b, fb, tolerance / 2 )
!$omp end critical (stack)    
      end if

   end do
!$omp critical (result)  
         integral_result = integral_result + private_result
!$omp end critical (result)       
   
!$omp end parallel

end function integral

end module integral_module

Finally it should be mentioned that the number of calls to the integrated function can be drastically reduced. For real life problems these function evaluations are usually by far the most time consuming parts. By storing function values at the interval limits in the stack as well, a lot of redundant computations can be avoided. Whereas so far three function evaluations have been executed during each call of the recursive function respectively during each loop iteration, it is only necessary to calculate the new function value at the midpoint of the current interval, if the function values at the interval boundaries a and b are already known from previous integration steps. Just in the beginning the first function evaluations at the initial interval boundaries have to be put on the stack.
The modifications to the serial version are only minor and can be easily adapted to the parallel versions as well.

integral_module.f95
module integral_module
use stack_m
public :: integral

contains

function integral (f, ah, bh, tolh)  &
      result (integral_result)

   implicit none
   interface
   function f (x) result (f_result)
      real, intent (in) :: x
      real :: f_result
   end function f
   end interface
   real, intent (in) :: ah, bh, tolh
   real :: a, b, tolerance, fa, fb, fmid
   real :: integral_result
   real :: h, mid
   real :: one_trapezoid_area, two_trapezoid_area
   real :: left_area, right_area
   type (stack_t) :: stack

   call new_stack ( stack )
   call push ( stack, ah, f(ah), bh, f(bh), tolh )
   
   integral_result = 0.0
   do
      if ( empty_stack ( stack ) ) exit
      call pop ( stack, a, fa, b, fb, tolerance )
   
      h = b - a
      mid = (a + b) /2
      fmid = f(mid)
      one_trapezoid_area = h * (fa + fb) / 2.0
      two_trapezoid_area = h/2 * (fa + fmid) / 2.0 + &
                           h/2 * (fmid + fb) / 2.0
      if (abs(one_trapezoid_area - two_trapezoid_area)  &
            < 3.0 * tolerance) then
         integral_result = integral_result + two_trapezoid_area
      else
         call push ( stack, a, fa, mid, fmid, tolerance / 2 )
         call push ( stack, mid, fmid, b, fb, tolerance / 2 )
      end if

   end do

end function integral

end module integral_module