This module contains procedures and generic interfaces for computing sum(x) accurately when x is a long array of wildly varying real or complex values.
More...

Data Types
type	fablocked_type
	This is a concrete derived type whose instances are exclusively used to request the Fast Accurate Blocked summation method of Blanchard, Higham, and Mary, 2020 within the generic interfaces of the ParaMonte library. More...

interface	getDot
	Generate and return the expression `dot_product(x, y)` accurately (almost without roundoff error). More...

interface	getSum
	Generate and return the expression `sum(x)` accurately (almost without roundoff error). More...

type	kahanbabu_type
	This is a concrete derived type whose instances are exclusively used to request the Kahan-Babuska summation method within the generic interfaces of the ParaMonte library. More...

type	nablocked_type
	This is a concrete derived type whose instances are exclusively used to request the Naive Blocked summation method of Blanchard, Higham, and Mary, 2020 within the generic interfaces of the ParaMonte library. More...

type	sum_type

Variables
character(*, SK), parameter	MODULE_NAME = "@pm_mathSum"

type(fablocked_type), parameter	fablocked = fablocked_type()
	This is a scalar `parameter` object of type fablocked_type. More...

type(nablocked_type), parameter	nablocked = nablocked_type()
	This is a scalar `parameter` object of type nablocked_type. More...

type(kahanbabu_type), parameter	kahanbabu = kahanbabu_type()
	This is a scalar `parameter` object of type kahanbabu_type. More...

Detailed Description

This module contains procedures and generic interfaces for computing sum(x) accurately when x is a long array of wildly varying real or complex values.

Some of the algorithms of this module are inspired by the work of Blanchard, Higham, and Mary, 2020, A Class of Fast and Accurate Summation Algorithms.

See also: pm_mathCumSum

Benchmarks:

Benchmark :: The effects of method on runtime efficiency ⛓

: The following program compares the runtime performance of getSum algorithms with the default Fortran sum() intrinsic function.

! Test the performance of `getSum()` with and without the selection `control` argument.
program benchmark
 
    use pm_bench, only: bench_type
    use pm_distUnif, only: setUnifRand
    use pm_mathCumSum, only: setCumSum
    use pm_arrayResize, only: setResized
    use pm_kind, only: SK, IK, LK, RKH, RK, RKG => RKD
    use iso_fortran_env, only: error_unit
 
    implicit none
 
    integer(IK)                         :: ibench
    integer(IK)                         :: iarr
    integer(IK)                         :: arrlen
    integer(IK)                         :: fileUnit
    real(RKG)                           :: dumsum = 0._RKG
    real(RKG)                           :: array(10**8)
    real(RKG)                           :: sumres
    type(bench_type)    , allocatable   :: bench(:)
 
    bench = [ bench_type(name = SK_"sum()", exec = getSumFortran, overhead = setOverhead) &
            , bench_type(name = SK_"fablocked", exec = getSumFAB, overhead = setOverhead) &
            , bench_type(name = SK_"nablocked", exec = getSumNAB, overhead = setOverhead) &
            , bench_type(name = SK_"kahanbabu", exec = getSumKAB, overhead = setOverhead) &
            , bench_type(name = SK_"iteration", exec = getSumIte, overhead = setOverhead) &
            , bench_type(name = SK_"recursion", exec = getSumRec, overhead = setOverhead) &
            ]
 
    write(*,"(*(g0,:,' '))")
    write(*,"(*(g0,:,' '))") "sum() vs. getSum()"
    write(*,"(*(g0,:,' '))")
 
    open(newunit = fileUnit, file = "main.out", status = "replace")
 
        call setUnifRand(array)
        write(fileUnit, "(*(g0,:,','))") "Array Size", (bench(ibench)%name, ibench = 1, size(bench))
        loopOverArraySize: do iarr = 2, 26, 2
 
            arrlen = 2**iarr
            write(*,"(*(g0,:,' '))") "Benchmarking with array size", arrlen
            do ibench = 1, size(bench)
                bench(ibench)%timing = bench(ibench)%getTiming(minsec = 0.07_RK)
            end do
            write(fileUnit,"(*(g0,:,','))") arrlen, (bench(ibench)%timing%mean, ibench = 1, size(bench))
 
        end do loopOverArraySize
        write(*,"(*(g0,:,' '))") dumsum
        write(*,"(*(g0,:,' '))")
 
    close(fileUnit)
 
contains
 
    !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    ! procedure wrappers.
    !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
    subroutine setOverhead()
        call getDummy()
    end subroutine
 
    subroutine getDummy()
        dumsum = dumsum + sumres
    end subroutine
 
    subroutine getSumFortran()
        sumres = sum(array(1 : arrlen))
        call getDummy()
    end subroutine
 
    subroutine getSumFAB()
        use pm_mathSum, only: getSum!, fablocked
        sumres = getSum(array(1 : arrlen))!, fablocked)
        call getDummy()
    end subroutine
 
    subroutine getSumNAB()
        use pm_mathSum, only: getSum, nablocked
        sumres = getSum(array(1 : arrlen), nablocked)
        call getDummy()
    end subroutine
 
    subroutine getSumKAB()
        use pm_mathSum, only: getSum, kahanbabu
        sumres = getSum(array(1 : arrlen), kahanbabu)
        call getDummy()
    end subroutine
 
    subroutine getSumIte()
        use pm_mathSum, only: getSum, iteration
        sumres = getSum(array(1 : arrlen), iteration)
        call getDummy()
    end subroutine
 
    subroutine getSumRec()
        use pm_mathSum, only: getSum, recursion
        sumres = getSum(array(1 : arrlen), recursion)
        call getDummy()
    end subroutine
 
end program benchmark

Example Unix compile command via Intel ifort compiler ⛓

#!/usr/bin/env sh
rm main.exe
ifort -fpp -standard-semantics -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
./main.exe

Example Windows Batch compile command via Intel ifort compiler ⛓

del main.exe
set PATH=..\..\..\lib;%PATH%
ifort /fpp /standard-semantics /O3 /I:..\..\..\include main.F90 ..\..\..\lib\libparamonte*.lib /exe:main.exe
main.exe

Example Unix / MinGW compile command via GNU gfortran compiler ⛓

#!/usr/bin/env sh
rm main.exe
gfortran -cpp -ffree-line-length-none -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
./main.exe

Postprocessing of the benchmark output ⛓

#!/usr/bin/env python
 
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
 
import os
dirname = os.path.basename(os.getcwd()) 
 
fontsize = 14
 
df = pd.read_csv("main.out", delimiter = ",")
colnames = list(df.columns.values)
 
 
 
ax = plt.figure(figsize = 1.25 * np.array([6.4,4.6]), dpi = 200)
ax = plt.subplot()
 
for colname in colnames[1:]:
    plt.plot( df[colnames[0]].values
            , df[colname].values
            , linewidth = 2
            )
 
plt.xticks(fontsize = fontsize)
plt.yticks(fontsize = fontsize)
ax.set_xlabel(colnames[0], fontsize = fontsize)
ax.set_ylabel("Runtime [ seconds ]", fontsize = fontsize)
ax.set_title(" vs. ".join(colnames[1:])+"\nLower is better.", fontsize = fontsize)
ax.set_xscale("log")
ax.set_yscale("log")
plt.minorticks_on()
plt.grid(visible = True, which = "both", axis = "both", color = "0.85", linestyle = "-")
ax.tick_params(axis = "y", which = "minor")
ax.tick_params(axis = "x", which = "minor")
ax.legend   ( colnames[1:]
           #, loc='center left'
           #, bbox_to_anchor=(1, 0.5)
            , fontsize = fontsize
            )
 
plt.tight_layout()
plt.savefig("benchmark." + dirname + ".runtime.png")
 
 
 
ax = plt.figure(figsize = 1.25 * np.array([6.4,4.6]), dpi = 200)
ax = plt.subplot()
 
plt.plot( df[colnames[0]].values
        , np.ones(len(df[colnames[0]].values))
        , linestyle = "--"
       #, color = "black"
        , linewidth = 2
        )
for colname in colnames[2:]:
    plt.plot( df[colnames[0]].values
            , df[colname].values / df[colnames[1]].values
            , linewidth = 2
            )
 
plt.xticks(fontsize = fontsize)
plt.yticks(fontsize = fontsize)
ax.set_xlabel(colnames[0], fontsize = fontsize)
ax.set_ylabel("Runtime compared to {}".format(colnames[1]), fontsize = fontsize)
ax.set_title("Runtime Ratio Comparison. Lower means faster.\nLower than 1 means faster than {}.".format(colnames[1]), fontsize = fontsize)
ax.set_xscale("log")
ax.set_yscale("log")
plt.minorticks_on()
plt.grid(visible = True, which = "both", axis = "both", color = "0.85", linestyle = "-")
ax.tick_params(axis = "y", which = "minor")
ax.tick_params(axis = "x", which = "minor")
ax.legend   ( colnames[1:]
           #, bbox_to_anchor = (1, 0.5)
           #, loc = "center left"
            , fontsize = fontsize
            )
 
plt.tight_layout()
plt.savefig("benchmark." + dirname + ".runtime.ratio.png")

Visualization of the benchmark output ⛓

Benchmark moral ⛓

Among all summation algorithms, fablocked_type appears to offer the most accurate result while also being even faster than the default Fortran sum() and all other implemented summation algorithms for array sizes > 100.

Benchmark :: The effects of method on runtime efficiency ⛓

: The following program compares the runtime performance of getDot algorithms with the default Fortran dot_product() intrinsic function.

! Test the performance of `getDot()` with and without the selection `control` argument.
program benchmark
 
    use pm_bench, only: bench_type
    use pm_distUnif, only: setUnifRand
    use pm_mathCumSum, only: setCumSum
    use pm_arrayResize, only: setResized
    use pm_kind, only: SK, IK, LK, RKH, RK, RKG => RKD
    use iso_fortran_env, only: error_unit
 
    implicit none
 
    integer(IK)                         :: ibench
    integer(IK)                         :: iarr
    integer(IK)                         :: arrlen
    integer(IK)                         :: fileUnit
    real(RKG)                           :: dumsum = 0._RKG
    real(RKG)                           :: array(10**8)
    real(RKG)                           :: dotres
    type(bench_type)    , allocatable   :: bench(:)
    logical(LK)                         :: underflowEnabled
 
    bench = [ bench_type(name = SK_"dot_product()", exec = getDotFortran, overhead = setOverhead) &
            , bench_type(name = SK_"fablocked", exec = getDotFAB, overhead = setOverhead) &
            , bench_type(name = SK_"nablocked", exec = getDotNAB, overhead = setOverhead) &
            , bench_type(name = SK_"kahanbabu", exec = getDotKAB, overhead = setOverhead) &
            , bench_type(name = SK_"iteration", exec = getDotIte, overhead = setOverhead) &
            , bench_type(name = SK_"recursion", exec = getDotRec, overhead = setOverhead) &
            ]
 
    write(*,"(*(g0,:,' '))")
    write(*,"(*(g0,:,' '))") "dot_product() vs. getDot()"
    write(*,"(*(g0,:,' '))")
 
    open(newunit = fileUnit, file = "main.out", status = "replace")
 
        call setUnifRand(array)
        !truth = array; call setCumSum(truth)
        write(fileUnit, "(*(g0,:,','))") "Array Size", (bench(ibench)%name, ibench = 1, size(bench))
        loopOverArraySize: do iarr = 2, 26, 2
 
            arrlen = 2**iarr
            write(*,"(*(g0,:,' '))") "Benchmarking with array size", arrlen
            do ibench = 1, size(bench)
                bench(ibench)%timing = bench(ibench)%getTiming(minsec = 0.07_RK)
            end do
            write(fileUnit,"(*(g0,:,','))") arrlen, (bench(ibench)%timing%mean, ibench = 1, size(bench))
 
        end do loopOverArraySize
        write(*,"(*(g0,:,' '))") dumsum
        write(*,"(*(g0,:,' '))")
 
    close(fileUnit)
 
contains
 
    !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    ! procedure wrappers.
    !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
    subroutine setOverhead()
        call getDummy()
    end subroutine
 
    subroutine getDummy()
        dumsum = dumsum + dotres
    end subroutine
 
    subroutine getDotFortran()
        dotres = dot_product(array(1 : arrlen), array(1 : arrlen))
        call getDummy()
    end subroutine
 
    subroutine getDotFAB()
        use pm_mathSum, only: getDot!, fablocked
        dotres = getDot(array(1 : arrlen), array(1 : arrlen))!, fablocked)
        call getDummy()
    end subroutine
 
    subroutine getDotNAB()
        use pm_mathSum, only: getDot, nablocked
        dotres = getDot(array(1 : arrlen), array(1 : arrlen), nablocked)
        call getDummy()
    end subroutine
 
    subroutine getDotKAB()
        use pm_mathSum, only: getDot, kahanbabu
        dotres = getDot(array(1 : arrlen), array(1 : arrlen), kahanbabu)
        call getDummy()
    end subroutine
 
    subroutine getDotIte()
        use pm_mathSum, only: getDot, iteration
        dotres = getDot(array(1 : arrlen), array(1 : arrlen), iteration)
        call getDummy()
    end subroutine
 
    subroutine getDotRec()
        use pm_mathSum, only: getDot, recursion
        dotres = getDot(array(1 : arrlen), array(1 : arrlen), recursion)
        call getDummy()
    end subroutine
 
end program benchmark

Example Unix compile command via Intel ifort compiler ⛓

#!/usr/bin/env sh
rm main.exe
ifort -fpp -standard-semantics -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
./main.exe

Example Windows Batch compile command via Intel ifort compiler ⛓

del main.exe
set PATH=..\..\..\lib;%PATH%
ifort /fpp /standard-semantics /O3 /I:..\..\..\include main.F90 ..\..\..\lib\libparamonte*.lib /exe:main.exe
main.exe

Example Unix / MinGW compile command via GNU gfortran compiler ⛓

#!/usr/bin/env sh
rm main.exe
gfortran -cpp -ffree-line-length-none -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
./main.exe

Postprocessing of the benchmark output ⛓

#!/usr/bin/env python
 
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
 
import os
dirname = os.path.basename(os.getcwd()) 
 
fontsize = 14
 
df = pd.read_csv("main.out", delimiter = ",")
colnames = list(df.columns.values)
 
 
 
ax = plt.figure(figsize = 1.25 * np.array([6.4,4.6]), dpi = 200)
ax = plt.subplot()
 
for colname in colnames[1:]:
    plt.plot( df[colnames[0]].values
            , df[colname].values
            , linewidth = 2
            )
 
plt.xticks(fontsize = fontsize)
plt.yticks(fontsize = fontsize)
ax.set_xlabel(colnames[0], fontsize = fontsize)
ax.set_ylabel("Runtime [ seconds ]", fontsize = fontsize)
ax.set_title(" vs. ".join(colnames[1:])+"\nLower is better.", fontsize = fontsize)
ax.set_xscale("log")
ax.set_yscale("log")
plt.minorticks_on()
plt.grid(visible = True, which = "both", axis = "both", color = "0.85", linestyle = "-")
ax.tick_params(axis = "y", which = "minor")
ax.tick_params(axis = "x", which = "minor")
ax.legend   ( colnames[1:]
           #, loc='center left'
           #, bbox_to_anchor=(1, 0.5)
            , fontsize = fontsize
            )
 
plt.tight_layout()
plt.savefig("benchmark." + dirname + ".runtime.png")
 
 
 
ax = plt.figure(figsize = 1.25 * np.array([6.4,4.6]), dpi = 200)
ax = plt.subplot()
 
plt.plot( df[colnames[0]].values
        , np.ones(len(df[colnames[0]].values))
        , linestyle = "--"
       #, color = "black"
        , linewidth = 2
        )
for colname in colnames[2:]:
    plt.plot( df[colnames[0]].values
            , df[colname].values / df[colnames[1]].values
            , linewidth = 2
            )
 
plt.xticks(fontsize = fontsize)
plt.yticks(fontsize = fontsize)
ax.set_xlabel(colnames[0], fontsize = fontsize)
ax.set_ylabel("Runtime compared to {}".format(colnames[1]), fontsize = fontsize)
ax.set_title("Runtime Ratio Comparison. Lower means faster.\nLower than 1 means faster than {}.".format(colnames[1]), fontsize = fontsize)
ax.set_xscale("log")
ax.set_yscale("log")
plt.minorticks_on()
plt.grid(visible = True, which = "both", axis = "both", color = "0.85", linestyle = "-")
ax.tick_params(axis = "y", which = "minor")
ax.tick_params(axis = "x", which = "minor")
ax.legend   ( colnames[1:]
           #, bbox_to_anchor = (1, 0.5)
           #, loc = "center left"
            , fontsize = fontsize
            )
 
plt.tight_layout()
plt.savefig("benchmark." + dirname + ".runtime.ratio.png")

Visualization of the benchmark output ⛓

Benchmark moral ⛓

Among all summation algorithms, fablocked_type appears to offer the most accurate result while also being even faster than the default Fortran dot_product() and all other implemented summation algorithms for array sizes > 100.

Final Remarks ⛓

If you believe this algorithm or its documentation can be improved, we appreciate your contribution and help to edit this page's documentation and source file on GitHub.
For details on the naming abbreviations, see this page.
For details on the naming conventions, see this page.
This software is distributed under the MIT license with additional terms outlined below.

If you use any parts or concepts from this library to any extent, please acknowledge the usage by citing the relevant publications of the ParaMonte library.
If you regenerate any parts/ideas from this library in a programming environment other than those currently supported by this ParaMonte library (i.e., other than C, C++, Fortran, MATLAB, Python, R), please also ask the end users to cite this original ParaMonte library.

This software is available to the public under a highly permissive license.
Help us justify its continued development and maintenance by acknowledging its benefit to society, distributing it, and contributing to it.

Copyright: Computational Data Science Lab

Test:: test_pm_mathSum

Author:: Amir Shahmoradi, August 8, 2024, 8:45 PM, NASA Goddard Space Flight Center, Washington, D.C.

Variable Documentation

◆ fablocked

type(fablocked_type), parameter pm_mathSum::fablocked = fablocked_type()

This is a scalar parameter object of type fablocked_type.

For example usage, see the documentation of the target procedure requiring this object.

See also: fablocked
kahanbabu
iteration
recursion
fablocked_type
kahanbabu_type
iteration_type
recursion_type

Final Remarks ⛓

If you believe this algorithm or its documentation can be improved, we appreciate your contribution and help to edit this page's documentation and source file on GitHub.
For details on the naming abbreviations, see this page.
For details on the naming conventions, see this page.
This software is distributed under the MIT license with additional terms outlined below.

If you use any parts or concepts from this library to any extent, please acknowledge the usage by citing the relevant publications of the ParaMonte library.
If you regenerate any parts/ideas from this library in a programming environment other than those currently supported by this ParaMonte library (i.e., other than C, C++, Fortran, MATLAB, Python, R), please also ask the end users to cite this original ParaMonte library.

This software is available to the public under a highly permissive license.
Help us justify its continued development and maintenance by acknowledging its benefit to society, distributing it, and contributing to it.

Copyright: Computational Data Science Lab

Author:: Amir Shahmoradi, August 8, 2024, 8:45 PM, NASA Goddard Space Flight Center, Washington, D.C.

Definition at line 143 of file pm_mathSum.F90.

◆ kahanbabu

type(kahanbabu_type), parameter pm_mathSum::kahanbabu = kahanbabu_type()

This is a scalar parameter object of type kahanbabu_type.

For example usage, see the documentation of the target procedure requiring this object.

See also: fablocked
kahanbabu
iteration
recursion
fablocked_type
kahanbabu_type
iteration_type
recursion_type

Final Remarks ⛓

If you believe this algorithm or its documentation can be improved, we appreciate your contribution and help to edit this page's documentation and source file on GitHub.
For details on the naming abbreviations, see this page.
For details on the naming conventions, see this page.
This software is distributed under the MIT license with additional terms outlined below.

If you use any parts or concepts from this library to any extent, please acknowledge the usage by citing the relevant publications of the ParaMonte library.
If you regenerate any parts/ideas from this library in a programming environment other than those currently supported by this ParaMonte library (i.e., other than C, C++, Fortran, MATLAB, Python, R), please also ask the end users to cite this original ParaMonte library.

This software is available to the public under a highly permissive license.
Help us justify its continued development and maintenance by acknowledging its benefit to society, distributing it, and contributing to it.

Copyright: Computational Data Science Lab

Author:: Amir Shahmoradi, August 8, 2024, 8:45 PM, NASA Goddard Space Flight Center, Washington, D.C.

Definition at line 266 of file pm_mathSum.F90.

◆ MODULE_NAME

character(*, SK), parameter pm_mathSum::MODULE_NAME = "@pm_mathSum"

Definition at line 80 of file pm_mathSum.F90.

◆ nablocked

type(nablocked_type), parameter pm_mathSum::nablocked = nablocked_type()

This is a scalar parameter object of type nablocked_type.

For example usage, see the documentation of the target procedure requiring this object.

See also: fablocked
kahanbabu
iteration
recursion
fablocked_type
kahanbabu_type
iteration_type
recursion_type

Final Remarks ⛓

If you believe this algorithm or its documentation can be improved, we appreciate your contribution and help to edit this page's documentation and source file on GitHub.
For details on the naming abbreviations, see this page.
For details on the naming conventions, see this page.
This software is distributed under the MIT license with additional terms outlined below.

If you use any parts or concepts from this library to any extent, please acknowledge the usage by citing the relevant publications of the ParaMonte library.
If you regenerate any parts/ideas from this library in a programming environment other than those currently supported by this ParaMonte library (i.e., other than C, C++, Fortran, MATLAB, Python, R), please also ask the end users to cite this original ParaMonte library.

This software is available to the public under a highly permissive license.
Help us justify its continued development and maintenance by acknowledging its benefit to society, distributing it, and contributing to it.

Copyright: Computational Data Science Lab

Author:: Amir Shahmoradi, August 8, 2024, 8:45 PM, NASA Goddard Space Flight Center, Washington, D.C.

Definition at line 206 of file pm_mathSum.F90.

Data Types

Variables

Detailed Description

Variable Documentation

◆ fablocked

◆ kahanbabu

◆ MODULE_NAME

◆ nablocked