Kahan summation algorithm の実装と検証

README.md

参考文献

• 数値計算のつぼ | 共立出版
• カハンの加算アルゴリズム - Wikipedia
• Kahan summation algorithm - Wikipedia

出力

               268435456 < length_data
 2.6843545600000001E+005 < sample_value * length_data
 2.6843545467426203E+005 < sum       ( sample_data(:) ) <- intrinsic
 2.6843545600000001E+005 < sum_kahan ( sample_data(:) )

mod_sum.f08

module mod_sum

    ! <module>s to import
    use, intrinsic :: iso_fortran_env

    ! require all variables to be explicitly declared
    implicit none



    ! accessibility of the <subroutine>s and <function>s in this <module>

    ! kind: subroutine
    private :: sum_kahan_1D_R64

    ! kind: interface
    public :: sum_kahan



    interface sum_kahan
        module procedure :: sum_kahan_1D_R64
    end interface sum_kahan



    ! contained <subroutine>s and <function>s are below
    contains



    pure function sum_kahan_1D_R64 (array) result(sum_computed)

        ! arguments for this <function>
        real(REAL64), intent(in) :: array(:)

        ! return value of this <function>
        real(REAL64) :: sum_computed

        ! variables for this <function>
        integer(INT32) :: size_array
        integer(INT32) :: itr_elm
        real(REAL64)   :: buffer
        real(REAL64)   :: compensator
        real(REAL64)   :: delta

        ! STEP.01
        ! initialize the variable for this <function>
        compensator  = 0.0_REAL64
        size_array   = size( array(:), kind=INT32 )
        sum_computed = array(1)

        ! STEP.02
        ! compute the sum using  Kahan summation algorithm
        do itr_elm = 2_INT32, size_array, 1_INT32
            delta        = array(itr_elm) - compensator
            buffer       = sum_computed + delta
            compensator  = (buffer - sum_computed) - delta
            sum_computed = buffer
        end do

        ! STEP.END
        return

    end function sum_kahan_1D_R64

end module mod_sum

test_sum_kahan.f08

program test

    ! <module>s to import
    use,     intrinsic :: iso_fortran_env
    use, non_intrinsic :: mod_sum




    ! require all variables to be explicitly declared
    implicit none



    ! constants for this <program>
    integer(INT32), parameter :: length_data  = 2_INT32 ** 28_INT32
    real(REAL64),   parameter :: sample_value = 1.0e-3_REAL64



    ! variables for this <program>
    real(REAL64), allocatable :: sample_data(:)



    ! STEP.01
    ! 1. allocate the array to store the sample data
    ! 2. input the sample data
    allocate( sample_data(length_data), source= sample_value )

    ! STEP.02
    ! show the computed result
    print '(I24,1X,A)'       , length_data                  , '< length_data'
    print '(ES24.16E3,1X,A)' , sample_value * length_data   , '< sample_value * length_data'
    print '(ES24.16E3,1X,A)' , sum       ( sample_data(:) ) , '< sum       ( sample_data(:) ) <- intrinsic'
    print '(ES24.16E3,1X,A)' , sum_kahan ( sample_data(:) ) , '< sum_kahan ( sample_data(:) )'

    ! STEP.03
    ! deallocate used array
    deallocate( sample_data )

end program test

test_sum_kahan.sh

gfortran ./mod_sum.f08        -c -ffree-line-length-none -std=f2008 -Wall -Wextra -O0
gfortran ./test_sum_kahan.f08 -c -ffree-line-length-none -std=f2008 -Wall -Wextra -O0
gfortran mod_sum.o test_sum_kahan.o