summation.f90 Source File

!> Provides utility routines for compenstated summation. The purpose
!> of these approaches is to try to minimise the accumulation of
!> floating point errors.
module summation
  implicit none
  private

  public :: kahan_sum, neumaier_sum, klein_sum, pairwise_sum

  interface kahan_sum
     module procedure kahan_sum_real_1d
  end interface kahan_sum

  interface neumaier_sum
     module procedure neumaier_sum_real_1d
  end interface neumaier_sum

  interface klein_sum
     module procedure klein_sum_real_1d
  end interface klein_sum

  interface pairwise_sum
     module procedure pairwise_sum_real_1d
  end interface pairwise_sum

contains

  !> Simple Kahan summation designed to limit the error arising from
  !> summing a large collection of numbers by taking care to keep
  !> track of the required compensation.
  pure real function kahan_sum_real_1d(array) result(the_result)
    implicit none
    real, dimension(:), intent(in) :: array
    integer :: i
    real :: compensation, compensated_element, current_total
    the_result = 0.0
    compensation = 0.0
    do i = 1, size(array)
       compensated_element = array(i) - compensation
       current_total = the_result + compensated_element
       compensation = (current_total - the_result) - compensated_element
       the_result = current_total
    end do
  end function kahan_sum_real_1d

  !> Very similar to Kahan summation but can perform better when the
  !> array contains a large range of inputs (i.e. where there may be
  !> elements that are larger the the cummulattive sum at that point)
  pure real function neumaier_sum_real_1d(array) result(the_result)
    implicit none
    real, dimension(:), intent(in) :: array
    integer :: i
    real :: compensation, current_total
    the_result = 0.0
    compensation = 0.0
    do i = 1, size(array)
       current_total = the_result + array(i)
       ! Here we choose how to calculate the compensation based on
       ! which is large of the running total and the current value
       if ( abs(the_result) >= abs(array(i)) ) then
          compensation = compensation + (the_result - current_total) + array(i)
       else
          compensation = compensation + (array(i) - current_total) + the_result
       end if
       the_result = current_total
    end do
    the_result = the_result + compensation
  end function neumaier_sum_real_1d

  !> Similar approach to Neumaier but introduces two levels of compensation
  pure real function klein_sum_real_1d(array) result(the_result)
    implicit none
    real, dimension(:), intent(in) :: array
    integer :: i
    real :: local_compensation_1, local_compensation_2, compensation_1, compensation_2
    real :: current_total
    the_result = 0.0
    compensation_1 = 0.0
    compensation_2 = 0.0

    do i = 1, size(array)
       ! This part is very similar to neumaier
       current_total = the_result + array(i)
       if ( abs(the_result) >= abs(array(i)) ) then
          local_compensation_1 = (the_result - current_total) + array(i)
       else
          local_compensation_1 = (array(i) - current_total) + the_result
       end if
       the_result = current_total

       ! Here we repeat the above approach but for accumulating the compensation
       current_total = compensation_1 + local_compensation_1
       if ( abs(compensation_1) >= abs(local_compensation_1) ) then
          local_compensation_2 = (compensation_1 - current_total) + local_compensation_1
       else
          local_compensation_2 = (local_compensation_1 - current_total) + compensation_1
       end if

       compensation_1 = current_total
       compensation_2 = compensation_2 + local_compensation_2
    end do
    the_result = the_result + compensation_1 + compensation_2
  end function klein_sum_real_1d

  !> Unlike Kahan summation and similar methods this approach simply
  !> divides the set of numbers in half and sums each half
  !> separately. This is applied recursively until the local problem size
  !> drops below a certain size, at which point regular summation is applied
  !>
  !> This should be more efficient than Kahan style approaches. Here
  !> there should be roughly the same number of additions as a simple
  !> sum, whilst Kahan style routines will have somewhat more (Kahan
  !> is roughly four times as much). We do have to worry a little
  !> about the cost of recursion here. Provided small_case_size is not
  !> too small then hopefully the recursion depth is not that large
  !> and hence we can neglect the cost of recursion.
  recursive real function pairwise_sum_real_1d(array) result(the_result)
    implicit none
    real, dimension(:), intent(in) :: array
    !> The small_case_size parameter controls the problem size for which the
    !> approach switches between recursive calls and a direct summation
    integer, parameter :: small_case_size = 5
    integer :: task_size, mid_point

    task_size = size(array)
    if (task_size <= small_case_size) then
       ! We could use one of the other methods provided here rather
       ! than the language provided `sum` but probably not required
       ! as small_case_size should be small.
       the_result = sum(array)
    else
       mid_point = task_size / 2
       the_result = pairwise_sum(array(:mid_point)) + pairwise_sum(array(mid_point+1:))
    end if
  end function pairwise_sum_real_1d

end module summation