accsum

accsum Commit Details

Date:2011-03-21 11:02:44 (7 years 1 month ago)
Author:Michael Baudin
Commit:15
Parents: 14
Message: * Fixed fast algorithms, added unit tests, help pages.
Changes:
R/help/en_US/accsum_priestx.xml → /help/en_US/support/accsum_priestx.xml
R/help/en_US/accsum_orderdynamic.xml → /help/en_US/support/accsum_orderdynamic.xml
A/src/cleaner_src.sce
A/tests/unit_tests/fdcs.tst
A/tests/unit_tests/fscs.tst
A/tests/unit_tests/fcompsum.tst
M/help/en_US/accsum_scs.xml
M/macros/accsum_scs.sci
M/sci_gateway/c/sci_accsum_fcompsum.c
M/help/en_US/accsum_fcompsum.xml
M/macros/accsum_shuffle.sci
M/help/en_US/accsum_twosum.xml
M/macros/accsum_priestx.sci
M/help/en_US/pseudomacros/accsum_fdcs.sci
M/help/en_US/accsum_dblcompsum.xml
M/macros/accsum_sumcond.sci
M/demos/ademo.sce
M/src/c/accsum.c
M/macros/accsum_order.sci
M/help/en_US/accsum_fasttwosum.xml
M/macros/accsum_wilkinson.sci
M/help/en_US/accsum_fdcs.xml
M/src/c/accsum.h
M/sci_gateway/c/sci_accsum_fdcs.c
M/help/en_US/pseudomacros/accsum_fscs.sci
M/help/en_US/accsum_compsum.xml
M/macros/accsum_compsum.sci
M/macros/accsum_straight.sci
M/sci_gateway/c/sci_accsum_fscs.c
M/help/en_US/accsum_fscs.xml
M/macros/accsum_twosum.sci
M/macros/accsum_dblcompsum.sci
M/help/en_US/accsum_dcs.xml
M/changelog.txt
M/macros/accsum_dcs.sci
M/help/en_US/pseudomacros/accsum_fcompsum.sci
M/macros/accsum_fasttwosum.sci
M/macros/accsum_orderdynamic.sci

File differences

tests/unit_tests/fcompsum.tst
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
// Copyright (C) 2011 - Michael Baudin
// Copyright (C) 2011 - DIGITEO - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
//
// <-- JVM NOT MANDATORY -->
// <-- ENGLISH IMPOSED -->
[s,e] = accsum_fcompsum ( [2 1] );
assert_checkequal(s,3);
assert_checkequal(e,0);
//
[s,e] = accsum_fcompsum ( [1 2] );
assert_checkequal(s,3);
assert_checkequal(e,0);
//
[s,e] = accsum_fcompsum ( [2 1] );
assert_checkequal(s,3);
assert_checkequal(e,0);
//
[s,e] = accsum_fcompsum ( [1 2] );
assert_checkequal(s,3);
assert_checkequal(e,0);
//
x = accsum_wilkinson(10);
[s,e] = accsum_fcompsum ( x );
assert_checkalmostequal(s,1024,1.e-12);
assert_checkalmostequal(e,-3.786e-14,[],1.e-10);
tests/unit_tests/fscs.tst
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
// Copyright (C) 2011 - Michael Baudin
// Copyright (C) 2011 - DIGITEO - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
//
// <-- JVM NOT MANDATORY -->
// <-- ENGLISH IMPOSED -->
[s,e] = accsum_fscs ( [2 1] );
assert_checkequal(s,3);
assert_checkequal(e,0);
//
[s,e] = accsum_fscs ( [1 2] );
assert_checkequal(s,3);
assert_checkequal(e,0);
//
[s,e] = accsum_fscs ( [2 1] );
assert_checkequal(s,3);
assert_checkequal(e,0);
//
[s,e] = accsum_fscs ( [1 2] );
assert_checkequal(s,3);
assert_checkequal(e,0);
//
x = accsum_wilkinson(10);
[s,e] = accsum_fscs ( x );
assert_checkalmostequal(s,1024,1.e-12);
assert_checkalmostequal(e,-3.786e-14,[],1.e-10);
tests/unit_tests/fdcs.tst
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
// Copyright (C) 2011 - Michael Baudin
// Copyright (C) 2011 - DIGITEO - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
//
// <-- JVM NOT MANDATORY -->
// <-- ENGLISH IMPOSED -->
[s,e] = accsum_fdcs ( [2 1] );
assert_checkequal(s,3);
assert_checkequal(e,0);
//
[s,e] = accsum_fdcs ( [1 2] );
assert_checkequal(s,3);
assert_checkequal(e,0);
//
[s,e] = accsum_fdcs ( [2 1] );
assert_checkequal(s,3);
assert_checkequal(e,0);
//
[s,e] = accsum_fdcs ( [1 2] );
assert_checkequal(s,3);
assert_checkequal(e,0);
//
x = accsum_wilkinson(10);
[s,e] = accsum_fdcs ( x );
assert_checkalmostequal(s,1024,1.e-12);
assert_checkalmostequal(e,-3.786e-14,[],1.e-10);
macros/accsum_dblcompsum.sci
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
// Copyright (C) 2010 - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
function [s,e] = accsum_dblcompsum ( x )
// The doubly compensated sum of a matrix.
//
// Calling Sequence
// [s,e] = accsum_dblcompsum ( x )
//
// Parameters
// x : a n-by-m matrix of doubles
// s : a 1-by-1 matrix of doubles, the highest significant digits of the sum
// e : a 1-by-1, matrix of doubles, the lowest significant digits of the sum
//
// Description
// Returns the sum of x.
// Algorithm 4.2 in SNAA (Doubly Compensated Summation)
// Due to Priest (and Kahan)
// Assumes that Scilab uses round-to-nearest.
//
// Examples
// [s,e] = accsum_dblcompsum ( [2 1] ) // 3
// [s,e] = accsum_dblcompsum ( [1 2] ) // 3
// // is 1+(%eps/2) but is 1 without algo
// [s,e] = accsum_dblcompsum ( [1 %eps/2 ] )
// // is 1+(%eps/2) but is 1 without algo
// [s,e] = accsum_dblcompsum ( [%eps/2 1] )
//
// Authors
// Michael Baudin, 2010
//
// Bibliography
// "Stability and numerical accuracy of algorithms", Nicolas Higham
// "Handbook of Floating Point Computations", Muller et al
// Sort x in decreasing magnitude order.
[B,k] = gsort(abs(x))
x = x(k)
// Perform the sum
s = x(1)
e = 0
n = size(x,"*")
for i = 2 : n
y = e + x(i)
u = x(i) - (y - e)
t = y + s
v = y - (t - s)
z = u + v
s = t + z
e = z - (s - t)
end
endfunction
// Copyright (C) 2010 - 2011 - Michael Baudin
// Copyright (C) 2011 - DIGITEO - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
function [s,e] = accsum_dblcompsum ( x )
// The doubly compensated sum of a matrix.
//
// Calling Sequence
// s = accsum_dblcompsum ( x )
// [s,e] = accsum_dblcompsum ( x )
//
// Parameters
// x : a n-by-m matrix of doubles
// s : a 1-by-1 matrix of doubles, the highest significant digits of the sum
// e : a 1-by-1, matrix of doubles, the lowest significant digits of the sum
//
// Description
// Returns the sum of x.
// Algorithm 4.2 in SNAA (Doubly Compensated Summation)
// Due to Priest (and Kahan)
// Assumes that Scilab uses round-to-nearest.
//
// Examples
// [s,e] = accsum_dblcompsum ( [2 1] ) // 3
// [s,e] = accsum_dblcompsum ( [1 2] ) // 3
// // is 1+(%eps/2) but is 1 without algo
// [s,e] = accsum_dblcompsum ( [1 %eps/2 ] )
// // is 1+(%eps/2) but is 1 without algo
// [s,e] = accsum_dblcompsum ( [%eps/2 1] )
//
// Authors
// Michael Baudin, 2010
//
// Bibliography
// "Stability and numerical accuracy of algorithms", Nicolas Higham
// "Handbook of Floating Point Computations", Muller et al
// Sort x in decreasing magnitude order.
[B,k] = gsort(abs(x))
x = x(k)
// Perform the sum
s = x(1)
e = 0
n = size(x,"*")
for i = 2 : n
y = e + x(i)
u = x(i) - (y - e)
t = y + s
v = y - (t - s)
z = u + v
s = t + z
e = z - (s - t)
end
endfunction
macros/accsum_sumcond.sci
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
// Copyright (C) 2010 - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
function [s,c] = accsum_sumcond ( x)
// Condition number of the sum function.
//
// Calling Sequence
// x = accsum_sumcond ( )
//
// Parameters
// x : a m-by-n, matrix of doubles
// s : a 1-by-1, matrix of doubles, the sum
// c : a 1-by-1, matrix of doubles, the condition number
//
// Description
// Condition number of the sum function.
//
// Examples
// // A straightforward example of an ill-conditionned data for the sum.
// xl = 10^(1:15);
// x = [-xl xl+0.1];
// sum(x)
// [s,c] = accsum_sumcond(x)
// for o = 1 : 8
// xo = accsum_order ( x , o );
// [s,e] = accsum_DCS ( xo );
// disp([o s e])
// end
//
// Authors
// Michael Baudin, 2010
//
// Bibliography
// "Stability and numerical accuracy of algorithms", Nicolas Higham
// "Handbook of Floating Point Computations", Muller et al
s = sum(x)
varf = sum(abs(x))
c = varf / abs(s)
endfunction
// Copyright (C) 2010 - 2011 - Michael Baudin
// Copyright (C) 2011 - DIGITEO - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
function [s,c] = accsum_sumcond ( x )
// Condition number of the sum function.
//
// Calling Sequence
// s = accsum_sumcond ( x )
// [s,c] = accsum_sumcond ( x )
//
// Parameters
// x : a m-by-n, matrix of doubles
// s : a 1-by-1, matrix of doubles, the sum
// c : a 1-by-1, matrix of doubles, the condition number
//
// Description
// Condition number of the sum function.
//
// Examples
// // A straightforward example of an ill-conditionned data for the sum.
// xl = 10^(1:15);
// x = [-xl xl+0.1];
// sum(x)
// [s,c] = accsum_sumcond(x)
// for o = 1 : 8
// xo = accsum_order ( x , o );
// [s,e] = accsum_DCS ( xo );
// disp([o s e])
// end
//
// Authors
// Michael Baudin, 2010
//
// Bibliography
// "Stability and numerical accuracy of algorithms", Nicolas Higham
// "Handbook of Floating Point Computations", Muller et al
s = sum(x)
varf = sum(abs(x))
c = varf / abs(s)
endfunction
macros/accsum_dcs.sci
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
// Copyright (C) 2010 - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
function [s,e] = accsum_dcs ( x )
// A Doubly Self Compensated Sum algorithm
//
// Calling Sequence
// [s,e] = accsum_dcs ( x )
//
// Parameters
// x : a m-by-n, matrix of doubles
// s : a 1-by-1 matrix of doubles, the highest significant digits of the sum
// e : a 1-by-1, matrix of doubles, the lowest significant digits of the sum
//
// Description
// A Doubly Self Compensated Sum algorithm.
// Uses accsum_fasttwosum.
//
// Examples
// [s,e] = accsum_dcs ( [2 1] ) // 3
// [s,e] = accsum_dcs ( [1 2] ) // 3
// [s,e] = accsum_dcs ( [2 1] ) // 3
// [s,e] = accsum_dcs ( [1 2] ) // 3
// x = accsum_wilkinson(10); size(x,"*")
// s = sum(x)
// [s,e] = accsum_dcs ( x )
// [s,e] = accsum_dcs ( x )
//
// Authors
// Michael Baudin, 2010
//
// Bibliography
// "Stability and numerical accuracy of algorithms", Nicolas Higham
// "Handbook of Floating Point Computations", Muller et al
// https://hpcrd.lbl.gov/SCG/ocean/NRS/ECMWF/img14.htm
// https://hpcrd.lbl.gov/SCG/ocean/NRS/SCSsum.F
s = 0
e = 0
n = size(x,"*")
for i = 1 : n
[s1,e1] = accsum_fasttwosum (e, x(i))
[s2,e2] = accsum_fasttwosum (s, s1)
[s,e] = accsum_fasttwosum (s2, e1 + e2)
end
endfunction
// Copyright (C) 2010 - 2011 - Michael Baudin
// Copyright (C) 2011 - DIGITEO - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
function [s,e] = accsum_dcs ( x )
// A Doubly Self Compensated Sum algorithm
//
// Calling Sequence
// s = accsum_dcs ( x )
// [s,e] = accsum_dcs ( x )
//
// Parameters
// x : a m-by-n, matrix of doubles
// s : a 1-by-1 matrix of doubles, the highest significant digits of the sum
// e : a 1-by-1, matrix of doubles, the lowest significant digits of the sum
//
// Description
// A Doubly Self Compensated Sum algorithm.
// Uses accsum_fasttwosum.
//
// Examples
// [s,e] = accsum_dcs ( [2 1] ) // 3
// [s,e] = accsum_dcs ( [1 2] ) // 3
// [s,e] = accsum_dcs ( [2 1] ) // 3
// [s,e] = accsum_dcs ( [1 2] ) // 3
// x = accsum_wilkinson(10); size(x,"*")
// s = sum(x)
// [s,e] = accsum_dcs ( x )
// [s,e] = accsum_dcs ( x )
//
// Authors
// Michael Baudin, 2010
//
// Bibliography
// "Stability and numerical accuracy of algorithms", Nicolas Higham
// "Handbook of Floating Point Computations", Muller et al
// https://hpcrd.lbl.gov/SCG/ocean/NRS/ECMWF/img14.htm
// https://hpcrd.lbl.gov/SCG/ocean/NRS/SCSsum.F
s = 0
e = 0
n = size(x,"*")
for i = 1 : n
[s1,e1] = accsum_fasttwosum (e, x(i))
[s2,e2] = accsum_fasttwosum (s, s1)
[s,e] = accsum_fasttwosum (s2, e1 + e2)
end
endfunction
macros/accsum_order.sci
1
1
2
23
34
45
// Copyright (C) 2010 - Michael Baudin
// Copyright (C) 2010 - 2011 - Michael Baudin
// Copyright (C) 2011 - DIGITEO - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
macros/accsum_fasttwosum.sci
1
1
2
23
34
45
......
1011
1112
1213
14
1315
1416
1517
// Copyright (C) 2010 - Michael Baudin
// Copyright (C) 2010 - 2011 - Michael Baudin
// Copyright (C) 2011 - DIGITEO - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// The fast2sum sum of a and b.
//
// Calling Sequence
// s = accsum_fasttwosum ( a , b )
// [s,t] = accsum_fasttwosum ( a , b )
//
// Parameters
macros/accsum_wilkinson.sci
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
// Copyright (C) 2010 - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
function x = accsum_wilkinson ( r )
// A test vector by Wilkinson.
//
// Calling Sequence
// x = accsum_wilkinson ( r )
//
// Parameters
// r : a 1-by-1 matrix of doubles, a positive integer
// x : a 1-by-r, matrix of doubles
//
// Description
// Returns a vector x of size 2^r where r is a flint.
// Exercise 4.2 in SNAA
//
// Examples
// // Example with r = 3
// x = accsum_wilkinson (3)
// // Expected result
// e = [1 , 1-%eps , 1-2*%eps , 1-2*%eps , 1-4*%eps , 1-4*%eps , 1-4*%eps , 1-4*%eps]'
//
// // Use with larger r
// x = accsum_wilkinson(10);
// s=sum(x)
// [s,e] = accsum_compsum ( x )
// [s,e] = accsum_dblcompsum ( x )
//
// Authors
// Michael Baudin, 2010
//
// Bibliography
// "Stability and numerical accuracy of algorithms", Nicolas Higham
// "Handbook of Floating Point Computations", Muller et al
n = 2^r
t = 53
u = 2^(-t) // would give the same as u = %eps
x(1) = 1
x(2) = 1 - u
for k = 2 : r
x(2^(k-1) + 1 : 2^k) = 1 - 2^(k-1) * u
end
endfunction
// Copyright (C) 2010 - 2011 - Michael Baudin
// Copyright (C) 2011 - DIGITEO - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
function x = accsum_wilkinson ( r )
// A test vector by Wilkinson.
//
// Calling Sequence
// x = accsum_wilkinson ( r )
//
// Parameters
// r : a 1-by-1 matrix of doubles, a positive integer
// x : a 1-by-r, matrix of doubles
//
// Description
// Returns a vector x of size 2^r where r is a flint.
// Exercise 4.2 in SNAA
//
// Examples
// // Example with r = 3
// x = accsum_wilkinson (3)
// // Expected result
// e = [1 , 1-%eps , 1-2*%eps , 1-2*%eps , 1-4*%eps , 1-4*%eps , 1-4*%eps , 1-4*%eps]'
//
// // Use with larger r
// x = accsum_wilkinson(10);
// s=sum(x)
// [s,e] = accsum_compsum ( x )
// [s,e] = accsum_dblcompsum ( x )
//
// Authors
// Michael Baudin, 2010
//
// Bibliography
// "Stability and numerical accuracy of algorithms", Nicolas Higham
// "Handbook of Floating Point Computations", Muller et al
n = 2^r
t = 53
u = 2^(-t) // would give the same as u = %eps
x(1) = 1
x(2) = 1 - u
for k = 2 : r
x(2^(k-1) + 1 : 2^k) = 1 - 2^(k-1) * u
end
endfunction
macros/accsum_orderdynamic.sci
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
// Copyright (C) 2010 - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
function s = accsum_orderdynamic ( x , order )
// Try to find an order which makes the sum perform bad/good.
//
// Calling Sequence
// s = accsum_orderdynamic ( x , order )
//
// Parameters
// x : a m-by-n, matrix of doubles
// order : a 1-by-1 matrix of doubles, a positive integer, the order. order = 1 : search for the entry which makes the smallest |s+y(i)|, good for accuracy, order = 2 : search for the entry which makes the largest |s+y(i)|, bad for accuracy
// s : the sum
//
// Description
// Try to find an order which makes the sum perform bad/good.
//
// Examples
// // Test on simple data
// x = [-9.9 5.5 2.2 -3.3 -6.6 0 1.1]
// sum(x)
// s = accsum_orderdynamic ( x , 1 )
// s = accsum_orderdynamic ( x , 2 )
//
// // Test on difficult data
// path = accsum_getpath ( );
// filename=fullfile(path,"demos","etaana.dat");
// x=fscanfMat(filename);
// // 7.390949249267578125000
// sum(x)
// // 0.3437500000000000000000
// s = accsum_orderdynamic ( x , 1 )
// // - 73.015625
// s = accsum_orderdynamic ( x , 2 )
//
// Authors
// Michael Baudin, 2010
//
// Bibliography
// "Stability and numerical accuracy of algorithms", Nicolas Higham
// "Handbook of Floating Point Computations", Muller et al
// https://hpcrd.lbl.gov/SCG/ocean/NRS/ECMWF/img16.htm
s = 0
z = x
// We remove an entry of z at each iteration of the loop:
// in the end, the array z is empty.
for k = 1 : size(x,"*")
[B,i] = gsort(abs(s+z),"g","i")
if ( order == 1 ) then
i = i(1)
else
i = i($)
end
s = s + z(i)
z(i) = []
end
endfunction
// Copyright (C) 2010 - 2011 - Michael Baudin
// Copyright (C) 2011 - DIGITEO - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
function s = accsum_orderdynamic ( x , order )
// Try to find an order which makes the sum perform bad/good.
//
// Calling Sequence
// s = accsum_orderdynamic ( x , order )
//
// Parameters
// x : a m-by-n, matrix of doubles
// order : a 1-by-1 matrix of doubles, a positive integer, the order. order = 1 : search for the entry which makes the smallest |s+y(i)|, good for accuracy, order = 2 : search for the entry which makes the largest |s+y(i)|, bad for accuracy
// s : a 1-by-1 matrix of doubles, the sum
//
// Description
// Try to find an order which makes the sum perform bad/good.
//
// Examples
// // Test on simple data
// x = [-9.9 5.5 2.2 -3.3 -6.6 0 1.1]
// sum(x)
// s = accsum_orderdynamic ( x , 1 )
// s = accsum_orderdynamic ( x , 2 )
//
// // Test on difficult data
// path = accsum_getpath ( );
// filename=fullfile(path,"demos","etaana.dat");
// x=fscanfMat(filename);
// // 7.390949249267578125000
// sum(x)
// // 0.3437500000000000000000
// s = accsum_orderdynamic ( x , 1 )
// // - 73.015625
// s = accsum_orderdynamic ( x , 2 )
//
// Authors
// Michael Baudin, 2010
//
// Bibliography
// "Stability and numerical accuracy of algorithms", Nicolas Higham
// "Handbook of Floating Point Computations", Muller et al
// https://hpcrd.lbl.gov/SCG/ocean/NRS/ECMWF/img16.htm
s = 0
z = x
// We remove an entry of z at each iteration of the loop:
// in the end, the array z is empty.
for k = 1 : size(x,"*")
[B,i] = gsort(abs(s+z),"g","i")
if ( order == 1 ) then
i = i(1)
else
i = i($)
end
s = s + z(i)
z(i) = []
end
endfunction
macros/accsum_compsum.sci
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
// Copyright (C) 2010 - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
function [s,e] = accsum_compsum ( x )
// The compensated sum of a matrix.
//
// Calling Sequence
// [s,e] = accsum_compsum ( x )
//
// Parameters
// x : a n-by-m matrix of doubles
// s : a 1-by-1 matrix of doubles, the highest significant digits of the sum
// e : a 1-by-1, matrix of doubles, the lowest significant digits of the sum
//
// Description
// Returns the sum of x.
// Algorithm 4.4 in HFCC
// Due to Knuth
// No assumption on a , b
// Assumes that Scilab uses round-to-nearest.
//
// Examples
// [s,e] = accsum_compsum ( [2 1] ) // 3
// [s,e] = accsum_compsum ( [1 2] ) // 3
// // is 1+(%eps/2) but is 1 without algo
// [s,e] = accsum_compsum ( [1 %eps/2 ] )
// // is 1+(%eps/2) but is 1 without algo
// [s,e] = accsum_compsum ( [%eps/2 1] )
//
// Authors
// Michael Baudin, 2010
//
// Bibliography
// "Stability and numerical accuracy of algorithms", Nicolas Higham
// "Handbook of Floating Point Computations", Muller et al
s = 0
e = 0
n = size(x,"*")
for i = 1 : n
t = s
y = x(i) + e
s = t + y
e = (t - s) + y
end
endfunction
// Copyright (C) 2010 - 2011 - Michael Baudin
// Copyright (C) 2011 - DIGITEO - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
function [s,e] = accsum_compsum ( x )
// The compensated sum of a matrix.
//
// Calling Sequence
// s = accsum_compsum ( x )
// [s,e] = accsum_compsum ( x )
//
// Parameters
// x : a n-by-m matrix of doubles
// s : a 1-by-1 matrix of doubles, the highest significant digits of the sum
// e : a 1-by-1, matrix of doubles, the lowest significant digits of the sum
//
// Description
// Returns s and e such that s+e is the sum of x.
// Algorithm 4.4 in HFCC
// Due to Knuth
// No assumption on a , b
// Assumes that Scilab uses round-to-nearest.
//
// Examples
// [s,e] = accsum_compsum ( [2 1] ) // 3
// [s,e] = accsum_compsum ( [1 2] ) // 3
// // is 1+(%eps/2) but is 1 without algo
// [s,e] = accsum_compsum ( [1 %eps/2 ] )
// // is 1+(%eps/2) but is 1 without algo
// [s,e] = accsum_compsum ( [%eps/2 1] )
//
// Authors
// Michael Baudin, 2010
//
// Bibliography
// "Stability and numerical accuracy of algorithms", Nicolas Higham
// "Handbook of Floating Point Computations", Muller et al
s = 0
e = 0
n = size(x,"*")
for i = 1 : n
t = s
y = x(i) + e
s = t + y
e = (t - s) + y
end
endfunction
macros/accsum_straight.sci
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
// Copyright (C) 2010 - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
function s = accsum_straight ( x )
// The straightforward sum of a matrix.
//
// Calling Sequence
// s = accsum_straight ( x )
//
// Parameters
// x : a n-by-m matrix of doubles
// s : a 1-by-1, matrix of doubles, the sum
//
// Description
// Computes the sum.
//
// Uses an unvectorized loop.
// This is a slow algorithm, used for comparison purposes
// only.
//
// Caution!
// This function may not return the same
// value as the "sum" function of Scilab!
// This is because the Intel MKL or ATLAS may
// reorder the data, processing it block by block
//
// Examples
// x = 1:10;
// accsum_straight(x)
//
// Authors
// Michael Baudin, 2010
//
// Bibliography
// "Stability and numerical accuracy of algorithms", Nicolas Higham
// "Handbook of Floating Point Computations", Muller et al
n = size(x,"*")
s = 0
for i = 1 : n
s = s + x(i)
end
endfunction
// Copyright (C) 2010 - 2011 - Michael Baudin
// Copyright (C) 2011 - DIGITEO - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
function s = accsum_straight ( x )
// The straightforward sum of a matrix.
//
// Calling Sequence
// s = accsum_straight ( x )
//
// Parameters
// x : a n-by-m matrix of doubles
// s : a 1-by-1, matrix of doubles, the sum
//
// Description
// Computes the sum.
//
// Uses an unvectorized loop.
// This is a slow algorithm, used for comparison purposes
// only.
//
// Caution!
// This function may not return the same
// value as the "sum" function of Scilab!
// This is because the Intel MKL or ATLAS may
// reorder the data, processing it block by block
//
// Examples
// x = 1:10;
// accsum_straight(x)
//
// Authors
// Michael Baudin, 2010
//
// Bibliography
// "Stability and numerical accuracy of algorithms", Nicolas Higham
// "Handbook of Floating Point Computations", Muller et al
n = size(x,"*")
s = 0
for i = 1 : n
s = s + x(i)
end
endfunction
macros/accsum_scs.sci
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
// Copyright (C) 2010 - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
function [s,e] = accsum_scs ( x )
// A Self Compensated Sum algorithm
//
// Calling Sequence
// [s,e] = accsum_scs ( x )
//
// Parameters
// x : a m-by-n, matrix of doubles
// s : a 1-by-1 matrix of doubles, the highest significant digits of the sum
// e : a 1-by-1, matrix of doubles, the lowest significant digits of the sum
//
// Description
// A Self Compensated Sum algorithm.
// Uses accsum_fasttwosum.
//
// Examples
// [s,e] = accsum_scs ( [2 1] ) // 3
// [s,e] = accsum_scs ( [1 2] ) // 3
// [s,e] = accsum_DCS ( [2 1] ) // 3
// [s,e] = accsum_DCS ( [1 2] ) // 3
// x = accsum_wilkinson(10); size(x,"*")
// s = sum(x)
// [s,e] = accsum_scs ( x )
// [s,e] = accsum_DCS ( x )
//
// Authors
// Michael Baudin, 2010
//
// Bibliography
// "Stability and numerical accuracy of algorithms", Nicolas Higham
// "Handbook of Floating Point Computations", Muller et al
// https://hpcrd.lbl.gov/SCG/ocean/NRS/ECMWF/img14.htm
// https://hpcrd.lbl.gov/SCG/ocean/NRS/SCSsum.F
s = 0
e = 0
n = size(x,"*")
for i = 1 : n
[s,e] = accsum_fasttwosum ( s , x(i) + e )
end
endfunction
// Copyright (C) 2010 - 2011 - Michael Baudin
// Copyright (C) 2011 - DIGITEO - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
function [s,e] = accsum_scs ( x )
// A Self Compensated Sum algorithm
//
// Calling Sequence
// s = accsum_scs ( x )
// [s,e] = accsum_scs ( x )
//
// Parameters
// x : a m-by-n, matrix of doubles
// s : a 1-by-1 matrix of doubles, the highest significant digits of the sum
// e : a 1-by-1, matrix of doubles, the lowest significant digits of the sum
//
// Description
// A Self Compensated Sum algorithm.
// Uses accsum_fasttwosum.
//
// Examples
// [s,e] = accsum_scs ( [2 1] ) // 3
// [s,e] = accsum_scs ( [1 2] ) // 3
// [s,e] = accsum_DCS ( [2 1] ) // 3
// [s,e] = accsum_DCS ( [1 2] ) // 3
// x = accsum_wilkinson(10); size(x,"*")
// s = sum(x)
// [s,e] = accsum_scs ( x )
// [s,e] = accsum_DCS ( x )
//
// Authors
// Michael Baudin, 2010
//
// Bibliography
// "Stability and numerical accuracy of algorithms", Nicolas Higham
// "Handbook of Floating Point Computations", Muller et al
// https://hpcrd.lbl.gov/SCG/ocean/NRS/ECMWF/img14.htm
// https://hpcrd.lbl.gov/SCG/ocean/NRS/SCSsum.F
s = 0
e = 0
n = size(x,"*")
for i = 1 : n
[s,e] = accsum_fasttwosum ( s , x(i) + e )
end
endfunction
macros/accsum_shuffle.sci
1
1
2
23
34
45
// Copyright (C) 2010 - Michael Baudin
// Copyright (C) 2010 - 2011 - Michael Baudin
// Copyright (C) 2011 - DIGITEO - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
macros/accsum_priestx.sci
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
// Copyright (C) 2010 - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
function x = accsum_priestx ( )
// A difficult example for SCS by Priest.
//
// Calling Sequence
// x = accsum_priestx ( )
//
// Parameters
// x : a 6-by-1, matrix of doubles
//
// Description
// "On properties of floating point arithmetics:
// numerical stability and the cost of accurate computations"
// Douglas Priest, 1992
// p62
// "What is the smallest n for which simply compensated
// summation is not guaranteed to yield a result
// with a small forward relative error in an arithmetic which
// conforms to the IEEE 754 standard. The smallest known
// to the author is ... Then the exact sum is 2, but the sum
// computed by the SCS is 3."
//
// Examples
// x = accsum_priestx ( )
// [s,e] = accsum_SCS ( x ) // = 3 + 0
// [s,e] = accsum_DCS ( x ) // = 2 + 0
//
// Authors
// Michael Baudin, 2010
//
// Bibliography
// "Stability and numerical accuracy of algorithms", Nicolas Higham
// "Handbook of Floating Point Computations", Muller et al
// "On properties of floating point arithmetics: numerical stability and the cost of accurate computations", Douglas Priest, 1992
t = 53
x(1) = 2^(t+1)
x(2) = 2^(t+1) - 2
x(3:6) = -(2^t - 1)
endfunction
// Copyright (C) 2010 - 2011 - Michael Baudin
// Copyright (C) 2011 - DIGITEO - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
function x = accsum_priestx ( )
// A difficult example for SCS by Priest.
//
// Calling Sequence
// x = accsum_priestx ( )
//
// Parameters
// x : a 6-by-1, matrix of doubles
//
// Description
// "On properties of floating point arithmetics:
// numerical stability and the cost of accurate computations"
// Douglas Priest, 1992
// p62
// "What is the smallest n for which simply compensated
// summation is not guaranteed to yield a result
// with a small forward relative error in an arithmetic which
// conforms to the IEEE 754 standard. The smallest known
// to the author is ... Then the exact sum is 2, but the sum
// computed by the SCS is 3."
//
// Examples
// x = accsum_priestx ( )
// [s,e] = accsum_SCS ( x ) // = 3 + 0
// [s,e] = accsum_DCS ( x ) // = 2 + 0
//
// Authors
// Michael Baudin, 2010
//
// Bibliography
// "Stability and numerical accuracy of algorithms", Nicolas Higham
// "Handbook of Floating Point Computations", Muller et al
// "On properties of floating point arithmetics: numerical stability and the cost of accurate computations", Douglas Priest, 1992
t = 53
x(1) = 2^(t+1)
x(2) = 2^(t+1) - 2
x(3:6) = -(2^t - 1)
endfunction
macros/accsum_twosum.sci
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
// Copyright (C) 2010 - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
function [s,t] = accsum_twosum ( a , b )
// The twosum sum of a and b.
//
// Calling Sequence
// [s,t] = accsum_twosum ( a , b )
//
// Parameters
// a : a 1-by-1 matrix of doubles
// b : a 1-by-1 matrix of doubles
// s : a 1-by-1 matrix of doubles, the highest significant digits of the sum
// t : a 1-by-1, matrix of doubles, the lowest significant digits of the sum
//
// Description
// Returns the sum of a and b
// so that a + b = s + t exactly.
// where s contains the highest significant digits
// and t the lowest significant digits.
// Algorithm 4.4 in HFCC
// Due to Knuth
// No assumption on a , b
// Assumes that Scilab uses round-to-nearest.
//
// Examples
// [s,t] = accsum_twosum ( 2 , 1 ) // 3
// [s,t] = accsum_twosum ( 1 , 2 ) // 3
// // is 1+(%eps/2) but is 1 without algo
// [s,t] = accsum_twosum ( 1 , %eps/2 )
// // is 1+(%eps/2) but is 1 without algo
// [s,t] = accsum_twosum ( %eps/2 , 1 )
//
// Authors
// Michael Baudin, 2010
//
// Bibliography
// "Stability and numerical accuracy of algorithms", Nicolas Higham
// "Handbook of Floating Point Computations", Muller et al
s = a + b
ap = s - b
bp = s - ap
da = a - ap
db = b - bp
t = da + db
endfunction
// Copyright (C) 2010 - 2011 - Michael Baudin
// Copyright (C) 2011 - DIGITEO - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
function [s,t] = accsum_twosum ( a , b )
// The twosum sum of a and b.
//
// Calling Sequence
// s = accsum_twosum ( a , b )
// [s,t] = accsum_twosum ( a , b )
//
// Parameters
// a : a 1-by-1 matrix of doubles
// b : a 1-by-1 matrix of doubles
// s : a 1-by-1 matrix of doubles, the highest significant digits of the sum
// t : a 1-by-1, matrix of doubles, the lowest significant digits of the sum
//
// Description
// Returns the sum of a and b
// so that a + b = s + t exactly.
// where s contains the highest significant digits
// and t the lowest significant digits.
// Algorithm 4.4 in HFCC
// Due to Knuth
// No assumption on a , b
// Assumes that Scilab uses round-to-nearest.
//
// Examples
// [s,t] = accsum_twosum ( 2 , 1 ) // 3
// [s,t] = accsum_twosum ( 1 , 2 ) // 3
// // is 1+(%eps/2) but is 1 without algo
// [s,t] = accsum_twosum ( 1 , %eps/2 )
// // is 1+(%eps/2) but is 1 without algo
// [s,t] = accsum_twosum ( %eps/2 , 1 )
//
// Authors
// Michael Baudin, 2010
//
// Bibliography
// "Stability and numerical accuracy of algorithms", Nicolas Higham
// "Handbook of Floating Point Computations", Muller et al
s = a + b
ap = s - b
bp = s - ap
da = a - ap
db = b - bp
t = da + db
endfunction
changelog.txt
77
88
99
10
1011
1112
* Added gateway for compiled algorithms.
* Created fdcs, fscs, fcomsum.
* Downcased DCS and SCS.
* Fixed fast algorithms, added unit tests, help pages.
-- Michael Baudin <michael.baudin@scilab.org> March 2011
demos/ademo.sce
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
// Copyright (C) 2010 - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
// A comment from the paper:
// "Among the 6 different orders, the one with decreasing
// magnitude order produces the correct answer. Intuitively in this
// order, the numbers with similar magnitude bug opposite signes are added
// together first, and then the result is added to the accumulation
// without much loss of precision. This order is the recommended
// order in summation of large arrays."
// Test this on sea surface height data :
// https://hpcrd.lbl.gov/SCG/ocean/NRS/
// "Using Accurate Arithmetics to Improve Numerical Reproducibility and Stability in
// Parallel Applications". Yun He and Chris H.Q. Ding. Journal of Supercomputing,
// Vol.18, Issue 3, 259-277, March 2001. Also Proceedings of International Conference
// on Supercomputing (ICS'00), May 2000, 225-234.
//
// Thanks to A. Edelman for pointing this paper.
path = accsum_getpath ( );
filename=fullfile(path,"demos","etaana.dat");
x=fscanfMat(filename);
size(x)
// Display the spread of the data
// - 1.074D+14 0.0009623 1.434D+13 1.080D+14
// This is a massive spread !
// The values are large, of opposite signs and are so that the
// sum is small. This is a worst case for the sum, as it
// maximizes the possibilities of cancellation.
[min(x) mean(x) st_deviation(x) max(x)]
s=sum(x)
[s,e] = accsum_compsum ( x )
[s,e] = accsum_dblcompsum ( x )
[s,e] = accsum_SCS ( x )
[s,e] = accsum_DCS ( x )
// sum(x) = 7.390625
// [s,e] = accsum_compsum ( x ) = 0.3823695 + 0.
// [s,e] = accsum_dblcompsum ( x ) = 0.3579858 + 0.
// [s,e] = accsum_SCS ( x ) = 0.3823695 + 0.
// [s,e] = accsum_DCS ( x ) = 0.3882289 + 0.
// Reproduce experiment at : https://hpcrd.lbl.gov/SCG/ocean/NRS/ECMWF/img17.htm
// Sensitivity regarding to summation order
mprintf("%-5s %-20s %-20s %-20s %-20s\n","Order","Straight","sum","SCS","DCS");
path = accsum_getpath ( );
filename=fullfile(path,"demos","etaana.dat");
for o = 1 : 8
x=fscanfMat(filename);
x = accsum_order ( x , o );
s0=accsum_straight(x);
s1=sum(x);
[s2,e] = accsum_SCS ( x );
[s3,e] = accsum_DCS ( x );
mprintf("#%-5d %+20.17f %+20.17f %+20.17f %+20.17f\n",o,s0,s1,s2,s3)
end
// Exact result is : +0.35798583924770355
// Order Straight sum SCS DCS
//#1 +34.41476821899414100 +7.39094924926757810 +0.38236951828002930 +0.38822889328002930
//#2 -70.64062500000000000 -0.64062500000000000 +0.35937500000000000 +0.35937500000000000
//#3 -73.01562500000000000 -0.39062500000000000 +0.35937500000000000 +0.35937500000000000
//#4 +13.85937500000000000 +1.85937500000000000 +0.35937500000000000 +0.35937500000000000
//#5 +14.31865964829921700 +1.31250000000000000 +0.35798583924770355 +0.35798583924770355
//#6 -36.25424389541149100 -2.24288563430309300 +0.35812444984912872 +0.35812444984912872
//#7 -66.64062500000000000 +22.10937500000000000 +0.35937500000000000 +0.35937500000000000
//#8 -0.09179687500000000 +0.30664062500000000 +0.39257812500000000 +0.36132812500000000
// The order #5, combined with simply (or double) compensated summation algo gives the
// correct answer.
// Conclusions:
// * Catastrophic cancellation may appear, for example when the spread of the data is extreme.
// * The compensated sum is efficient for stabilizing results.
// * When the values are all positive, ordering the values by increasing
// magnitude would give the best result.
// * Ordering the data by decreasing magnitude helps to further reduce the error,
// when the values have mixed signs.
// * The sum function is sensitive to the use of the compiler optimizations,
// the Intel MKL or ATLAS libraries, which changes the orders of the computation.
// In Scilab, we use a straightforward sum function :
// http://gitweb.scilab.org/?p=scilab.git;a=blob;f=scilab/modules/elementary_functions/src/fortran/dsum.f;h=3e216367a176911f46a16a771ead8cfb8620f74c;hb=HEAD
// The reason of the difference between the Scilab macro "Straight" and the sum
// function is the use of compiler optimizations.
// Furthermore, the processor may use intermediate extended precision.
//////////////////////////////////////////////////////
// Make a Monte-Carlo simulation for the sum
lines(0);
kmax = 100;
path = accsum_getpath ( );
filename=fullfile(path,"demos","etaana.dat");
x=fscanfMat(filename);
for k = 1 : kmax
if ( modulo(k,10) == 0 ) then
disp([k min(sumk) mean(sumk) max(sumk)]);
end
x = accsum_order ( x , 8 );
s=accsum_straight(x);
sumk(k) = s;
end
disp([min(sumk) mean(sumk) max(sumk)])
//========= E N D === O F === D E M O =========//
//
// Load this script into the editor
//
filename = "ademo.sce";
dname = get_absolute_file_path(filename);
editor ( fullfile(dname,filename) );
// Copyright (C) 2010 - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
// A comment from the paper:
// "Among the 6 different orders, the one with decreasing
// magnitude order produces the correct answer. Intuitively in this
// order, the numbers with similar magnitude bug opposite signes are added
// together first, and then the result is added to the accumulation
// without much loss of precision. This order is the recommended
// order in summation of large arrays."
// Test this on sea surface height data :
// https://hpcrd.lbl.gov/SCG/ocean/NRS/
// "Using Accurate Arithmetics to Improve Numerical Reproducibility and Stability in
// Parallel Applications". Yun He and Chris H.Q. Ding. Journal of Supercomputing,
// Vol.18, Issue 3, 259-277, March 2001. Also Proceedings of International Conference
// on Supercomputing (ICS'00), May 2000, 225-234.
//
// Thanks to A. Edelman for pointing this paper.
path = accsum_getpath ( );
filename=fullfile(path,"demos","etaana.dat");
x=fscanfMat(filename);
size(x)
// Display the spread of the data
// - 1.074D+14 0.0009623 1.434D+13 1.080D+14
// This is a massive spread !
// The values are large, of opposite signs and are so that the
// sum is small. This is a worst case for the sum, as it
// maximizes the possibilities of cancellation.
[min(x) mean(x) st_deviation(x) max(x)]
s=sum(x)
[s,e] = accsum_compsum ( x )
[s,e] = accsum_dblcompsum ( x )
[s,e] = accsum_scs ( x )
[s,e] = accsum_dcs ( x )
// sum(x) = 7.390625
// [s,e] = accsum_compsum ( x ) = 0.3823695 + 0.
// [s,e] = accsum_dblcompsum ( x ) = 0.3579858 + 0.
// [s,e] = accsum_scs ( x ) = 0.3823695 + 0.
// [s,e] = accsum_dcs ( x ) = 0.3882289 + 0.
// Reproduce experiment at : https://hpcrd.lbl.gov/SCG/ocean/NRS/ECMWF/img17.htm
// Sensitivity regarding to summation order
mprintf("%-5s %-20s %-20s %-20s %-20s\n","Order","Straight","sum","SCS","DCS");
path = accsum_getpath ( );
filename=fullfile(path,"demos","etaana.dat");
for o = 1 : 8
x=fscanfMat(filename);
x = accsum_order ( x , o );
s0=accsum_straight(x);
s1=sum(x);
[s2,e] = accsum_scs ( x );
[s3,e] = accsum_dcs ( x );
mprintf("#%-5d %+20.17f %+20.17f %+20.17f %+20.17f\n",o,s0,s1,s2,s3)
end
// Exact result is : +0.35798583924770355
// Order Straight sum SCS DCS
//#1 +34.41476821899414100 +7.39094924926757810 +0.38236951828002930 +0.38822889328002930
//#2 -70.64062500000000000 -0.64062500000000000 +0.35937500000000000 +0.35937500000000000
//#3 -73.01562500000000000 -0.39062500000000000 +0.35937500000000000 +0.35937500000000000
//#4 +13.85937500000000000 +1.85937500000000000 +0.35937500000000000 +0.35937500000000000
//#5 +14.31865964829921700 +1.31250000000000000 +0.35798583924770355 +0.35798583924770355
//#6 -36.25424389541149100 -2.24288563430309300 +0.35812444984912872 +0.35812444984912872
//#7 -66.64062500000000000 +22.10937500000000000 +0.35937500000000000 +0.35937500000000000
//#8 -0.09179687500000000 +0.30664062500000000 +0.39257812500000000 +0.36132812500000000
// The order #5, combined with simply (or double) compensated summation algo gives the
// correct answer.
// Conclusions:
// * Catastrophic cancellation may appear, for example when the spread of the data is extreme.
// * The compensated sum is efficient for stabilizing results.
// * When the values are all positive, ordering the values by increasing
// magnitude would give the best result.
// * Ordering the data by decreasing magnitude helps to further reduce the error,
// when the values have mixed signs.
// * The sum function is sensitive to the use of the compiler optimizations,
// the Intel MKL or ATLAS libraries, which changes the orders of the computation.
// In Scilab, we use a straightforward sum function :
// http://gitweb.scilab.org/?p=scilab.git;a=blob;f=scilab/modules/elementary_functions/src/fortran/dsum.f;h=3e216367a176911f46a16a771ead8cfb8620f74c;hb=HEAD
// The reason of the difference between the Scilab macro "Straight" and the sum
// function is the use of compiler optimizations.
// Furthermore, the processor may use intermediate extended precision.
//////////////////////////////////////////////////////
// Make a Monte-Carlo simulation for the sum
lines(0);
kmax = 100;
path = accsum_getpath ( );
filename=fullfile(path,"demos","etaana.dat");
x=fscanfMat(filename);
for k = 1 : kmax
if ( modulo(k,10) == 0 ) then
disp([k min(sumk) mean(sumk) max(sumk)]);
end
x = accsum_order ( x , 8 );
s=accsum_straight(x);
sumk(k) = s;
end
disp([min(sumk) mean(sumk) max(sumk)])
//========= E N D === O F === D E M O =========//
//
// Load this script into the editor
//
filename = "ademo.sce";
dname = get_absolute_file_path(filename);
editor ( fullfile(dname,filename) );
src/cleaner_src.sce
1
2
3
4
5
6
7
8
9
10
11
12
13
mode(-1)
// ====================================================================
// Copyright (C) INRIA - Serge Steer
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
// ====================================================================
src_path=get_absolute_file_path('cleaner_src.sce');
exec(src_path+'c/cleaner.sce',-1)
clear src_path
src/c/accsum.c
1
1
2
23
34
45
......
2324
2425
2526
26
27
28
29
2730
2831
2932
......
3134
3235
3336
34
37
3538
36
3739
3840
3941
......
4143
4244
4345
44
46
4547
4648
47
49
4850
49
51
5052
5153
5254
53
55
56
57
5458
5559
5660
......
5862
5963
6064
61
65
6266
63
6467
68
6569
6670
6771
68
72
6973
7074
71
75
7276
77
7378
7479
7580
76
77
78
81
7982
8083
8184
......
130133
131134
132135
133
136
134137
135138
136139
137
140
141
142
138143
139144
140145
......
143148
144149
145150
151
152
153
154
146155
147
156
148157
149158
150159
151
160
152161
153
162
154163
155164
156165
// Copyright (C) 2010 - Michael Baudin
// Copyright (C) 2010-2011 - Michael Baudin
// Copyright (C) 2011 - DIGITEO - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// on Supercomputing (ICS'00), May 2000, 225-234.
// Sets s as the sum of x entries, with a double compensated summation.
// Returns s and e such that s+e is the sum of x.
// s : the highest significant digits of the sum
// e : the lowest significant digits of the sum
// A Doubly Self Compensated Sum algorithm.
// Uses accsum_fasttwosum.
// Bibliography
// "Handbook of Floating Point Computations", Muller et al
// https://hpcrd.lbl.gov/SCG/ocean/NRS/ECMWF/img14.htm
// https://hpcrd.lbl.gov/SCG/ocean/NRS/SCSsum.F
void accsum_fdcs(double *x, int n, double *s)
void accsum_fdcs(double *x, int n, double *s, double *e)
{
double e;
double s1;
double e1;
double s2;
int i;
*s = 0;
e = 0;
*e = 0;
for (i = 0; i< n; i++)
{
accsum_fasttwosum (e, x[i], &s1, &e1 );
accsum_fasttwosum (*e, x[i], &s1, &e1 );
accsum_fasttwosum (*s, s1, &s2, &e2);
accsum_fasttwosum (s2, e1 + e2, s, &e);
accsum_fasttwosum (s2, e1 + e2, s, e);
}
}
// Sets s as the sum of x entries, with a simply compensated summation.
// Returns s and e such that s+e is the sum of x.
// s : the highest significant digits of the sum
// e : the lowest significant digits of the sum
// A Self Compensated Sum algorithm.
// Uses accsum_fasttwosum.
// Bibliography
// "Handbook of Floating Point Computations", Muller et al
// https://hpcrd.lbl.gov/SCG/ocean/NRS/ECMWF/img14.htm
// https://hpcrd.lbl.gov/SCG/ocean/NRS/SCSsum.F
void accsum_fscs(double *x, int n, double *s)
void accsum_fscs(double *x, int n, double *s, double *e)
{
double e;
double s1;
double e1;
int i;
*s = 0;
e = 0;
*e = 0;
for (i = 0; i< n; i++)
{
accsum_fasttwosum ( *s , x[i] + e , &s1, &e);
accsum_fasttwosum ( *s , x[i] + *e , &s1, &e1);
*s = s1;
*e = e1;
}
}
// Sets s+e as the sum of a and b.
// Returns the sum of a and b
// so that a + b = s + t exactly.
// Returns s and t so that a + b = s + t exactly.
// where s contains the highest significant digits
// and t the lowest significant digits.
//
*t = da + db;
}
void accsum_fcompsum ( double * x , int n , double * s )
void accsum_fcompsum ( double * x , int n , double * s , double * e )
{
// The compensated sum of a matrix.
// Description
// Returns the sum of x.
// Returns s and e such that s+e is the sum of x.
// s : the highest significant digits of the sum
// e : the lowest significant digits of the sum
// Algorithm 4.4 in HFCC
// Due to Knuth
// No assumption on a , b
// "Stability and numerical accuracy of algorithms", Nicolas Higham
// "Handbook of Floating Point Computations", Muller et al
int i;
double t;
double y;
*s = 0;
e = 0;
*e = 0;
for (i = 0; i< n; i++)
{
t = *s;
y = x[i] + e;
y = x[i] + *e;
*s = t + y;
e = (t - *s) + y;
*e = (t - *s) + y;
}
}
src/c/accsum.h
3232
3333
3434
35
35
3636
3737
3838
39
39
4040
41
41
4242
4343
4444
__BEGIN_DECLS
// Sets s as the sum of x entries, with a double compensated summation.
void accsum_fdcs(double *x, int n, double *s);
void accsum_fdcs(double *x, int n, double *s, double *e);
// Sets s+e as the sum of a and b
void accsum_fasttwosum(double a, double b, double *s, double *e);
// Sets s as the sum of x entries, with a simply compensated summation.
void accsum_fscs(double *x, int n, double *s);
void accsum_fscs(double *x, int n, double *s, double *e);
// The compensated sum of a matrix.
void accsum_fcompsum ( double * x , int n , double * s );
void accsum_fcompsum ( double * x , int n , double * s , double *e);
__END_DECLS
help/en_US/accsum_orderdynamic.xml
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
<?xml version="1.0" encoding="UTF-8"?>
<!--
*
* This help file was generated from accsum_orderdynamic.sci using help_from_sci().
*
-->
<refentry version="5.0-subset Scilab" xml:id="accsum_orderdynamic" xml:lang="en"
xmlns="http://docbook.org/ns/docbook"
xmlns:xlink="http://www.w3.org/1999/xlink"
xmlns:svg="http://www.w3.org/2000/svg"
xmlns:ns3="http://www.w3.org/1999/xhtml"
xmlns:mml="http://www.w3.org/1998/Math/MathML"
xmlns:db="http://docbook.org/ns/docbook">
<refnamediv>
<refname>accsum_orderdynamic</refname><refpurpose>Try to find an order which makes the sum perform bad/good.</refpurpose>
</refnamediv>
<refsynopsisdiv>
<title>Calling Sequence</title>
<synopsis>
s = accsum_orderdynamic ( x , order )
</synopsis>
</refsynopsisdiv>
<refsection>
<title>Parameters</title>
<variablelist>
<varlistentry><term>x :</term>
<listitem><para> a m-by-n, matrix of doubles</para></listitem></varlistentry>
<varlistentry><term>order :</term>
<listitem><para> a 1-by-1 matrix of doubles, a positive integer, the order. order = 1 : search for the entry which makes the smallest |s+y(i)|, good for accuracy, order = 2 : search for the entry which makes the largest |s+y(i)|, bad for accuracy</para></listitem></varlistentry>
<varlistentry><term>s :</term>
<listitem><para> a 1-by-1 matrix of doubles, the sum</para></listitem></varlistentry>
</variablelist>
</refsection>
<refsection>
<title>Description</title>
<para>
Try to find an order which makes the sum perform bad/good.
</para>
<para>
</para>
</refsection>
<refsection>
<title>Examples</title>
<programlisting role="example"><![CDATA[
// Test on simple data
x = [-9.9 5.5 2.2 -3.3 -6.6 0 1.1]
sum(x)
s = accsum_orderdynamic ( x , 1 )
s = accsum_orderdynamic ( x , 2 )
// Test on difficult data
path = accsum_getpath ( );
filename=fullfile(path,"demos","etaana.dat");
x=fscanfMat(filename);
// 7.390949249267578125000
sum(x)
// 0.3437500000000000000000
s = accsum_orderdynamic ( x , 1 )
// - 73.015625
s = accsum_orderdynamic ( x , 2 )
]]></programlisting>
</refsection>
<refsection>
<title>Authors</title>
<simplelist type="vert">
<member>Michael Baudin, 2010</member>
</simplelist>
</refsection>
<refsection>
<title>Bibliography</title>
<para>"Stability and numerical accuracy of algorithms", Nicolas Higham</para>
<para>"Handbook of Floating Point Computations", Muller et al</para>
<para>https://hpcrd.lbl.gov/SCG/ocean/NRS/ECMWF/img16.htm</para>
</refsection>
</refentry>
help/en_US/accsum_priestx.xml
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
<?xml version="1.0" encoding="UTF-8"?>
<!--
*
* This help file was generated from accsum_priestx.sci using help_from_sci().
*
-->
<refentry version="5.0-subset Scilab" xml:id="accsum_priestx" xml:lang="en"
xmlns="http://docbook.org/ns/docbook"
xmlns:xlink="http://www.w3.org/1999/xlink"
xmlns:svg="http://www.w3.org/2000/svg"
xmlns:ns3="http://www.w3.org/1999/xhtml"
xmlns:mml="http://www.w3.org/1998/Math/MathML"
xmlns:db="http://docbook.org/ns/docbook">
<refnamediv>
<refname>accsum_priestx</refname><refpurpose>A difficult example for SCS by Priest.</refpurpose>
</refnamediv>
<refsynopsisdiv>
<title>Calling Sequence</title>
<synopsis>
x = accsum_priestx ( )
</synopsis>
</refsynopsisdiv>
<refsection>
<title>Parameters</title>
<variablelist>
<varlistentry><term>x :</term>
<listitem><para> a 6-by-1, matrix of doubles</para></listitem></varlistentry>
</variablelist>
</refsection>
<refsection>
<title>Description</title>
<para>
"On properties of floating point arithmetics:
numerical stability and the cost of accurate computations"
Douglas Priest, 1992
p62
"What is the smallest n for which simply compensated
summation is not guaranteed to yield a result
with a small forward relative error in an arithmetic which
conforms to the IEEE 754 standard. The smallest known
to the author is ... Then the exact sum is 2, but the sum
computed by the SCS is 3."
</para>
<para>
</para>
</refsection>
<refsection>
<title>Examples</title>
<programlisting role="example"><![CDATA[
x = accsum_priestx ( )
[s,e] = accsum_SCS ( x ) // = 3 + 0
[s,e] = accsum_DCS ( x ) // = 2 + 0
]]></programlisting>
</refsection>
</refentry>
help/en_US/accsum_scs.xml
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
<?xml version="1.0" encoding="UTF-8"?>
<!--
*
* This help file was generated from accsum_scs.sci using help_from_sci().
*
-->
<refentry version="5.0-subset Scilab" xml:id="accsum_scs" xml:lang="en"
xmlns="http://docbook.org/ns/docbook"
xmlns:xlink="http://www.w3.org/1999/xlink"
xmlns:svg="http://www.w3.org/2000/svg"
xmlns:ns3="http://www.w3.org/1999/xhtml"
xmlns:mml="http://www.w3.org/1998/Math/MathML"
xmlns:db="http://docbook.org/ns/docbook">
<refnamediv>
<refname>accsum_scs</refname><refpurpose>A Self Compensated Sum algorithm</refpurpose>
</refnamediv>
<refsynopsisdiv>
<title>Calling Sequence</title>
<synopsis>
[s,e] = accsum_scs ( x )
</synopsis>
</refsynopsisdiv>
<refsection>
<title>Parameters</title>
<variablelist>
<varlistentry><term>x :</term>
<listitem><para> a m-by-n, matrix of doubles</para></listitem></varlistentry>
<varlistentry><term>s :</term>
<listitem><para> a 1-by-1 matrix of doubles, the highest significant digits of the sum</para></listitem></varlistentry>
<varlistentry><term>e :</term>
<listitem><para> a 1-by-1, matrix of doubles, the lowest significant digits of the sum</para></listitem></varlistentry>
</variablelist>
</refsection>
<refsection>
<title>Description</title>
<para>
A Self Compensated Sum algorithm.
Uses accsum_fasttwosum.
</para>
<para>
</para>
</refsection>
<refsection>
<title>Examples</title>
<programlisting role="example"><![CDATA[
[s,e] = accsum_scs ( [2 1] ) // 3
[s,e] = accsum_scs ( [1 2] ) // 3
[s,e] = accsum_DCS ( [2 1] ) // 3
[s,e] = accsum_DCS ( [1 2] ) // 3
x = accsum_wilkinson(10); size(x,"*")
s = sum(x)
[s,e] = accsum_scs ( x )
[s,e] = accsum_DCS ( x )
]]></programlisting>
</refsection>
<refsection>
<title>Authors</title>
<simplelist type="vert">
<member>Michael Baudin, 2010</member>
</simplelist>
</refsection>
<refsection>
<title>Bibliography</title>
<para>"Stability and numerical accuracy of algorithms", Nicolas Higham</para>
<para>"Handbook of Floating Point Computations", Muller et al</para>
<para>https://hpcrd.lbl.gov/SCG/ocean/NRS/ECMWF/img14.htm</para>
<para>https://hpcrd.lbl.gov/SCG/ocean/NRS/SCSsum.F</para>
</refsection>
</refentry>
<?xml version="1.0" encoding="UTF-8"?>
<!--
*
* This help file was generated from accsum_scs.sci using help_from_sci().
*
-->
<refentry version="5.0-subset Scilab" xml:id="accsum_scs" xml:lang="en"
xmlns="http://docbook.org/ns/docbook"
xmlns:xlink="http://www.w3.org/1999/xlink"
xmlns:svg="http://www.w3.org/2000/svg"
xmlns:ns3="http://www.w3.org/1999/xhtml"
xmlns:mml="http://www.w3.org/1998/Math/MathML"
xmlns:db="http://docbook.org/ns/docbook">
<refnamediv>
<refname>accsum_scs</refname><refpurpose>A Self Compensated Sum algorithm</refpurpose>
</refnamediv>
<refsynopsisdiv>
<title>Calling Sequence</title>
<synopsis>
s = accsum_scs ( x )
[s,e] = accsum_scs ( x )
</synopsis>
</refsynopsisdiv>
<refsection>
<title>Parameters</title>
<variablelist>
<varlistentry><term>x :</term>
<listitem><para> a m-by-n, matrix of doubles</para></listitem></varlistentry>
<varlistentry><term>s :</term>
<listitem><para> a 1-by-1 matrix of doubles, the highest significant digits of the sum</para></listitem></varlistentry>
<varlistentry><term>e :</term>
<listitem><para> a 1-by-1, matrix of doubles, the lowest significant digits of the sum</para></listitem></varlistentry>
</variablelist>
</refsection>
<refsection>
<title>Description</title>
<para>
A Self Compensated Sum algorithm.
Uses accsum_fasttwosum.
</para>
<para>
</para>
</refsection>
<refsection>
<title>Examples</title>
<programlisting role="example"><![CDATA[
[s,e] = accsum_scs ( [2 1] ) // 3
[s,e] = accsum_scs ( [1 2] ) // 3
[s,e] = accsum_DCS ( [2 1] ) // 3
[s,e] = accsum_DCS ( [1 2] ) // 3
x = accsum_wilkinson(10); size(x,"*")
s = sum(x)
[s,e] = accsum_scs ( x )
[s,e] = accsum_DCS ( x )
]]></programlisting>
</refsection>
<refsection>
<title>Authors</title>
<simplelist type="vert">
<member>Michael Baudin, 2010</member>
</simplelist>
</refsection>
<refsection>
<title>Bibliography</title>
<para>"Stability and numerical accuracy of algorithms", Nicolas Higham</para>
<para>"Handbook of Floating Point Computations", Muller et al</para>
<para>https://hpcrd.lbl.gov/SCG/ocean/NRS/ECMWF/img14.htm</para>
<para>https://hpcrd.lbl.gov/SCG/ocean/NRS/SCSsum.F</para>
</refsection>
</refentry>
help/en_US/accsum_dcs.xml
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
<?xml version="1.0" encoding="UTF-8"?>
<!--
*
* This help file was generated from accsum_dcs.sci using help_from_sci().
*
-->
<refentry version="5.0-subset Scilab" xml:id="accsum_dcs" xml:lang="en"
xmlns="http://docbook.org/ns/docbook"
xmlns:xlink="http://www.w3.org/1999/xlink"
xmlns:svg="http://www.w3.org/2000/svg"
xmlns:ns3="http://www.w3.org/1999/xhtml"
xmlns:mml="http://www.w3.org/1998/Math/MathML"
xmlns:db="http://docbook.org/ns/docbook">
<refnamediv>
<refname>accsum_dcs</refname><refpurpose>A Doubly Self Compensated Sum algorithm</refpurpose>
</refnamediv>
<refsynopsisdiv>
<title>Calling Sequence</title>
<synopsis>
[s,e] = accsum_dcs ( x )
</synopsis>
</refsynopsisdiv>
<refsection>
<title>Parameters</title>
<variablelist>
<varlistentry><term>x :</term>
<listitem><para> a m-by-n, matrix of doubles</para></listitem></varlistentry>
<varlistentry><term>s :</term>
<listitem><para> a 1-by-1 matrix of doubles, the highest significant digits of the sum</para></listitem></varlistentry>
<varlistentry><term>e :</term>
<listitem><para> a 1-by-1, matrix of doubles, the lowest significant digits of the sum</para></listitem></varlistentry>
</variablelist>
</refsection>
<refsection>
<title>Description</title>
<para>
A Doubly Self Compensated Sum algorithm.
Uses accsum_fasttwosum.
</para>
<para>
</para>
</refsection>
<refsection>
<title>Examples</title>
<programlisting role="example"><![CDATA[
[s,e] = accsum_dcs ( [2 1] ) // 3
[s,e] = accsum_dcs ( [1 2] ) // 3
[s,e] = accsum_dcs ( [2 1] ) // 3
[s,e] = accsum_dcs ( [1 2] ) // 3
x = accsum_wilkinson(10); size(x,"*")
s = sum(x)
[s,e] = accsum_dcs ( x )
[s,e] = accsum_dcs ( x )
]]></programlisting>
</refsection>
<refsection>
<title>Authors</title>
<simplelist type="vert">
<member>Michael Baudin, 2010</member>
</simplelist>
</refsection>
<refsection>
<title>Bibliography</title>
<para>"Stability and numerical accuracy of algorithms", Nicolas Higham</para>
<para>"Handbook of Floating Point Computations", Muller et al</para>
<para>https://hpcrd.lbl.gov/SCG/ocean/NRS/ECMWF/img14.htm</para>
<para>https://hpcrd.lbl.gov/SCG/ocean/NRS/SCSsum.F</para>
</refsection>
</refentry>
<?xml version="1.0" encoding="UTF-8"?>
<!--
*
* This help file was generated from accsum_dcs.sci using help_from_sci().
*
-->
<refentry version="5.0-subset Scilab" xml:id="accsum_dcs" xml:lang="en"
xmlns="http://docbook.org/ns/docbook"
xmlns:xlink="http://www.w3.org/1999/xlink"
xmlns:svg="http://www.w3.org/2000/svg"
xmlns:ns3="http://www.w3.org/1999/xhtml"
xmlns:mml="http://www.w3.org/1998/Math/MathML"
xmlns:db="http://docbook.org/ns/docbook">
<refnamediv>
<refname>accsum_dcs</refname><refpurpose>A Doubly Self Compensated Sum algorithm</refpurpose>
</refnamediv>
<refsynopsisdiv>
<title>Calling Sequence</title>
<synopsis>
s = accsum_dcs ( x )
[s,e] = accsum_dcs ( x )
</synopsis>
</refsynopsisdiv>
<refsection>
<title>Parameters</title>
<variablelist>
<varlistentry><term>x :</term>
<listitem><para> a m-by-n, matrix of doubles</para></listitem></varlistentry>
<varlistentry><term>s :</term>
<listitem><para> a 1-by-1 matrix of doubles, the highest significant digits of the sum</para></listitem></varlistentry>
<varlistentry><term>e :</term>
<listitem><para> a 1-by-1, matrix of doubles, the lowest significant digits of the sum</para></listitem></varlistentry>
</variablelist>
</refsection>
<refsection>
<title>Description</title>
<para>
A Doubly Self Compensated Sum algorithm.
Uses accsum_fasttwosum.
</para>
<para>
</para>
</refsection>
<refsection>
<title>Examples</title>
<programlisting role="example"><![CDATA[
[s,e] = accsum_dcs ( [2 1] ) // 3
[s,e] = accsum_dcs ( [1 2] ) // 3
[s,e] = accsum_dcs ( [2 1] ) // 3
[s,e] = accsum_dcs ( [1 2] ) // 3
x = accsum_wilkinson(10); size(x,"*")
s = sum(x)
[s,e] = accsum_dcs ( x )
[s,e] = accsum_dcs ( x )
]]></programlisting>
</refsection>
<refsection>
<title>Authors</title>
<simplelist type="vert">
<member>Michael Baudin, 2010</member>
</simplelist>
</refsection>
<refsection>
<title>Bibliography</title>
<para>"Stability and numerical accuracy of algorithms", Nicolas Higham</para>
<para>"Handbook of Floating Point Computations", Muller et al</para>
<para>https://hpcrd.lbl.gov/SCG/ocean/NRS/ECMWF/img14.htm</para>
<para>https://hpcrd.lbl.gov/SCG/ocean/NRS/SCSsum.F</para>
</refsection>
</refentry>
help/en_US/accsum_fcompsum.xml
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
<?xml version="1.0" encoding="UTF-8"?>
<!--
*
* This help file was generated from accsum_fcompsum.sci using help_from_sci().
*
-->
<refentry version="5.0-subset Scilab" xml:id="accsum_fcompsum" xml:lang="en"
xmlns="http://docbook.org/ns/docbook"
xmlns:xlink="http://www.w3.org/1999/xlink"
xmlns:svg="http://www.w3.org/2000/svg"
xmlns:ns3="http://www.w3.org/1999/xhtml"
xmlns:mml="http://www.w3.org/1998/Math/MathML"
xmlns:db="http://docbook.org/ns/docbook">
<refnamediv>
<refname>accsum_fcompsum</refname><refpurpose>The compensated sum of a matrix.</refpurpose>
</refnamediv>
<refsynopsisdiv>
<title>Calling Sequence</title>
<synopsis>
s = accsum_fcompsum ( x )
</synopsis>
</refsynopsisdiv>
<refsection>
<title>Parameters</title>
<variablelist>
<varlistentry><term>x :</term>
<listitem><para> a n-by-m matrix of doubles</para></listitem></varlistentry>
<varlistentry><term>s :</term>
<listitem><para> a 1-by-1 matrix of doubles, the highest significant digits of the sum</para></listitem></varlistentry>
</variablelist>
</refsection>
<refsection>
<title>Description</title>
<para>
Returns the sum of x.
Algorithm 4.4 in HFCC
Due to Knuth
No assumption on a , b
Assumes that we uses round-to-nearest.
This is a fast implementation, based on compiled source code..
</para>
<para>
</para>
</refsection>
<refsection>
<title>Examples</title>
<programlisting role="example"><![CDATA[
s = accsum_fcompsum ( [2 1] ) // 3
s = accsum_fcompsum ( [1 2] ) // 3
// is 1+(%eps/2) but is 1 without algo
s = accsum_fcompsum ( [1 %eps/2 ] )
// is 1+(%eps/2) but is 1 without algo
s = accsum_fcompsum ( [%eps/2 1] )
]]></programlisting>
</refsection>
<refsection>
<title>Authors</title>
<simplelist type="vert">
<member>Michael Baudin, 2010-2011</member>
</simplelist>
</refsection>
<refsection>
<title>Bibliography</title>
<para>"Stability and numerical accuracy of algorithms", Nicolas Higham</para>
<para>"Handbook of Floating Point Computations", Muller et al</para>
</refsection>
</refentry>
<?xml version="1.0" encoding="UTF-8"?>
<!--
*
* This help file was generated from accsum_fcompsum.sci using help_from_sci().
*
-->
<refentry version="5.0-subset Scilab" xml:id="accsum_fcompsum" xml:lang="en"
xmlns="http://docbook.org/ns/docbook"
xmlns:xlink="http://www.w3.org/1999/xlink"
xmlns:svg="http://www.w3.org/2000/svg"
xmlns:ns3="http://www.w3.org/1999/xhtml"
xmlns:mml="http://www.w3.org/1998/Math/MathML"
xmlns:db="http://docbook.org/ns/docbook">
<refnamediv>
<refname>accsum_fcompsum</refname><refpurpose>The compensated sum of a matrix.</refpurpose>
</refnamediv>
<refsynopsisdiv>
<title>Calling Sequence</title>
<synopsis>
s = accsum_fcompsum ( x )
[s,e] = accsum_fcompsum ( x )
</synopsis>
</refsynopsisdiv>
<refsection>
<title>Parameters</title>
<variablelist>
<varlistentry><term>x :</term>
<listitem><para> a n-by-m matrix of doubles</para></listitem></varlistentry>
<varlistentry><term>s :</term>
<listitem><para> a 1-by-1 matrix of doubles, the highest significant digits of the sum</para></listitem></varlistentry>
<varlistentry><term>e :</term>
<listitem><para> a 1-by-1, matrix of doubles, the lowest significant digits of the sum</para></listitem></varlistentry>
</variablelist>
</refsection>
<refsection>
<title>Description</title>
<para>
Returns the sum of x.
Algorithm 4.4 in HFCC
Due to Knuth
No assumption on a , b
Assumes that we uses round-to-nearest.
This is a fast implementation, based on compiled source code..
</para>
<para>
</para>
</refsection>
<refsection>
<title>Examples</title>
<programlisting role="example"><![CDATA[
s = accsum_fcompsum ( [2 1] ) // 3
s = accsum_fcompsum ( [1 2] ) // 3
// is 1+(%eps/2) but is 1 without algo
s = accsum_fcompsum ( [1 %eps/2 ] )
// is 1+(%eps/2) but is 1 without algo
s = accsum_fcompsum ( [%eps/2 1] )
]]></programlisting>
</refsection>
<refsection>
<title>Authors</title>
<simplelist type="vert">
<member>Michael Baudin, 2010-2011</member>
</simplelist>
</refsection>
<refsection>
<title>Bibliography</title>
<para>"Stability and numerical accuracy of algorithms", Nicolas Higham</para>
<para>"Handbook of Floating Point Computations", Muller et al</para>
</refsection>
</refentry>
help/en_US/pseudomacros/accsum_fcompsum.sci
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
// Copyright (C) 2010-2011 - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
function s = accsum_fcompsum ( x )
// The compensated sum of a matrix.
//
// Calling Sequence
// s = accsum_fcompsum ( x )
//
// Parameters
// x : a n-by-m matrix of doubles
// s : a 1-by-1 matrix of doubles, the highest significant digits of the sum
//
// Description
// Returns the sum of x.
// Algorithm 4.4 in HFCC
// Due to Knuth
// No assumption on a , b
// Assumes that we uses round-to-nearest.
// This is a fast implementation, based on compiled source code..
//
// Examples
// s = accsum_fcompsum ( [2 1] ) // 3
// s = accsum_fcompsum ( [1 2] ) // 3
// // is 1+(%eps/2) but is 1 without algo
// s = accsum_fcompsum ( [1 %eps/2 ] )
// // is 1+(%eps/2) but is 1 without algo
// s = accsum_fcompsum ( [%eps/2 1] )
//
// Authors
// Michael Baudin, 2010-2011
//
// Bibliography
// "Stability and numerical accuracy of algorithms", Nicolas Higham
// "Handbook of Floating Point Computations", Muller et al
endfunction
// Copyright (C) 2010-2011 - Michael Baudin
// Copyright (C) 2011 - DIGITEO - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
function [s,e] = accsum_fcompsum ( x )
// The compensated sum of a matrix.
//
// Calling Sequence
// s = accsum_fcompsum ( x )
// [s,e] = accsum_fcompsum ( x )
//
// Parameters
// x : a n-by-m matrix of doubles
// s : a 1-by-1 matrix of doubles, the highest significant digits of the sum
// e : a 1-by-1, matrix of doubles, the lowest significant digits of the sum
//
// Description
// Returns the sum of x.
// Algorithm 4.4 in HFCC
// Due to Knuth
// No assumption on a , b
// Assumes that we uses round-to-nearest.
// This is a fast implementation, based on compiled source code..
//
// Examples
// s = accsum_fcompsum ( [2 1] ) // 3
// s = accsum_fcompsum ( [1 2] ) // 3
// // is 1+(%eps/2) but is 1 without algo
// s = accsum_fcompsum ( [1 %eps/2 ] )
// // is 1+(%eps/2) but is 1 without algo
// s = accsum_fcompsum ( [%eps/2 1] )
//
// Authors
// Michael Baudin, 2010-2011
//
// Bibliography
// "Stability and numerical accuracy of algorithms", Nicolas Higham
// "Handbook of Floating Point Computations", Muller et al
endfunction
help/en_US/pseudomacros/accsum_fscs.sci
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
// Copyright (C) 2010-2011 - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
function s = accsum_fscs ( x )
// A Self Compensated Sum algorithm
//
// Calling Sequence
// s = accsum_fscs ( x )
//
// Parameters
// x : a m-by-n, matrix of doubles
// s : a 1-by-1 matrix of doubles, the highest significant digits of the sum
//
// Description
// A Self Compensated Sum algorithm.
// Uses the fasttwosum algorithm.
// This is a fast implementation, based on compiled source code.
//
// Examples
// s = accsum_fscs ( [2 1] ) // 3
// s = accsum_fscs ( [1 2] ) // 3
// x = accsum_wilkinson(10); size(x,"*")
// s = sum(x)
// s = accsum_fscs ( x )
//
// Authors
// Michael Baudin, 2010-2011
//
// Bibliography
// "Stability and numerical accuracy of algorithms", Nicolas Higham
// "Handbook of Floating Point Computations", Muller et al
// https://hpcrd.lbl.gov/SCG/ocean/NRS/ECMWF/img14.htm
// https://hpcrd.lbl.gov/SCG/ocean/NRS/SCSsum.F
endfunction
// Copyright (C) 2010-2011 - Michael Baudin
// Copyright (C) 2011 - DIGITEO - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
function [s,e] = accsum_fscs ( x )
// A Self Compensated Sum algorithm
//
// Calling Sequence
// s = accsum_fscs ( x )
// [s,e] = accsum_fscs ( x )
//
// Parameters
// x : a m-by-n, matrix of doubles
// s : a 1-by-1 matrix of doubles, the highest significant digits of the sum
// e : a 1-by-1, matrix of doubles, the lowest significant digits of the sum
//
// Description
// A Self Compensated Sum algorithm.
// Uses the fasttwosum algorithm.
// This is a fast implementation, based on compiled source code.
//
// Examples
// s = accsum_fscs ( [2 1] ) // 3
// s = accsum_fscs ( [1 2] ) // 3
// x = accsum_wilkinson(10); size(x,"*")
// s = sum(x)
// s = accsum_fscs ( x )
//
// Authors
// Michael Baudin, 2010-2011
//
// Bibliography
// "Stability and numerical accuracy of algorithms", Nicolas Higham
// "Handbook of Floating Point Computations", Muller et al
// https://hpcrd.lbl.gov/SCG/ocean/NRS/ECMWF/img14.htm
// https://hpcrd.lbl.gov/SCG/ocean/NRS/SCSsum.F
endfunction
help/en_US/pseudomacros/accsum_fdcs.sci
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
// Copyright (C) 2010-2011 - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
function s = accsum_fdcs ( x )
// A Doubly Self Compensated Sum algorithm
//
// Calling Sequence
// s = accsum_fdcs ( x )
//
// Parameters
// x : a m-by-n, matrix of doubles
// s : a 1-by-1 matrix of doubles, the highest significant digits of the sum
//
// Description
// A Doubly Self Compensated Sum algorithm.
// Uses the fasttwosum algorithm.
// This is a fast implementation, based on compiled source code.
//
// Examples
// s = accsum_fdcs ( [2 1] ) // 3
// s = accsum_fdcs ( [1 2] ) // 3
// s = accsum_fdcs ( [2 1] ) // 3
// s = accsum_fdcs ( [1 2] ) // 3
// x = accsum_wilkinson(10); size(x,"*")
// s = sum(x)
// s = accsum_fdcs ( x )
// s = accsum_fdcs ( x )
//
// Authors
// Michael Baudin, 2010-2011
//
// Bibliography
// "Stability and numerical accuracy of algorithms", Nicolas Higham
// "Handbook of Floating Point Computations", Muller et al
// https://hpcrd.lbl.gov/SCG/ocean/NRS/ECMWF/img14.htm
// https://hpcrd.lbl.gov/SCG/ocean/NRS/SCSsum.F
endfunction
// Copyright (C) 2010-2011 - Michael Baudin
// Copyright (C) 2011 - DIGITEO - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
function [s,e] = accsum_fdcs ( x )
// A Doubly Self Compensated Sum algorithm
//
// Calling Sequence
// s = accsum_fdcs ( x )
// [s,e] = accsum_fdcs ( x )
//
// Parameters
// x : a m-by-n, matrix of doubles
// s : a 1-by-1 matrix of doubles, the highest significant digits of the sum
// e : a 1-by-1, matrix of doubles, the lowest significant digits of the sum
//
// Description
// A Doubly Self Compensated Sum algorithm.
// Uses the fasttwosum algorithm.
// This is a fast implementation, based on compiled source code.
//
// Examples
// s = accsum_fdcs ( [2 1] ) // 3
// s = accsum_fdcs ( [1 2] ) // 3
// s = accsum_fdcs ( [2 1] ) // 3
// s = accsum_fdcs ( [1 2] ) // 3
// x = accsum_wilkinson(10); size(x,"*")
// s = sum(x)
// s = accsum_fdcs ( x )
// s = accsum_fdcs ( x )
//
// Authors
// Michael Baudin, 2010-2011
//
// Bibliography
// "Stability and numerical accuracy of algorithms", Nicolas Higham
// "Handbook of Floating Point Computations", Muller et al
// https://hpcrd.lbl.gov/SCG/ocean/NRS/ECMWF/img14.htm
// https://hpcrd.lbl.gov/SCG/ocean/NRS/SCSsum.F
endfunction
help/en_US/accsum_fasttwosum.xml
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
<?xml version="1.0" encoding="UTF-8"?>
<!--
*
* This help file was generated from accsum_fasttwosum.sci using help_from_sci().
*
-->
<refentry version="5.0-subset Scilab" xml:id="accsum_fasttwosum" xml:lang="en"
xmlns="http://docbook.org/ns/docbook"
xmlns:xlink="http://www.w3.org/1999/xlink"
xmlns:svg="http://www.w3.org/2000/svg"
xmlns:ns3="http://www.w3.org/1999/xhtml"
xmlns:mml="http://www.w3.org/1998/Math/MathML"
xmlns:db="http://docbook.org/ns/docbook">
<refnamediv>
<refname>accsum_fasttwosum</refname><refpurpose>The fast2sum sum of a and b.</refpurpose>
</refnamediv>
<refsynopsisdiv>
<title>Calling Sequence</title>
<synopsis>
[s,t] = accsum_fasttwosum ( a , b )
</synopsis>
</refsynopsisdiv>
<refsection>
<title>Parameters</title>
<variablelist>
<varlistentry><term>a :</term>
<listitem><para> a 1-by-1 matrix of doubles</para></listitem></varlistentry>
<varlistentry><term>b :</term>
<listitem><para> a 1-by-1 matrix of doubles</para></listitem></varlistentry>
<varlistentry><term>s :</term>
<listitem><para> a 1-by-1 matrix of doubles, the highest significant digits of the sum</para></listitem></varlistentry>
<varlistentry><term>t :</term>
<listitem><para> a 1-by-1, matrix of doubles, the lowest significant digits of the sum</para></listitem></varlistentry>
</variablelist>
</refsection>
<refsection>
<title>Description</title>
<para>
Returns the sum of a and b
so that a + b = s + t exactly.
where s contains the highest significant digits
and t the lowest significant digits.
</para>
<para>
Assumes that |a|&gt;=|b|. If not, the arguments are switched.
</para>
<para>
Algorithm 4.3 in HFCC
Due to Dekker (and Kahan)
Assumes that Scilab uses round-to-nearest.
</para>
<para>
</para>
</refsection>
<refsection>
<title>Examples</title>
<programlisting role="example"><![CDATA[
[s,t] = accsum_fasttwosum ( 2 , 1 ) // 3
// is 1+(%eps/2) but is 1 without algo
[s,t] = accsum_fasttwosum ( 1 , %eps/2 )
]]></programlisting>
</refsection>
<refsection>
<title>Authors</title>
<simplelist type="vert">
<member>Michael Baudin, 2010</member>
</simplelist>
</refsection>
<refsection>
<title>Bibliography</title>
<para>"Stability and numerical accuracy of algorithms", Nicolas Higham</para>
<para>"Handbook of Floating Point Computations", Muller et al</para>
</refsection>
</refentry>
<?xml version="1.0" encoding="UTF-8"?>
<!--
*
* This help file was generated from accsum_fasttwosum.sci using help_from_sci().
*
-->
<refentry version="5.0-subset Scilab" xml:id="accsum_fasttwosum" xml:lang="en"
xmlns="http://docbook.org/ns/docbook"
xmlns:xlink="http://www.w3.org/1999/xlink"
xmlns:svg="http://www.w3.org/2000/svg"
xmlns:ns3="http://www.w3.org/1999/xhtml"
xmlns:mml="http://www.w3.org/1998/Math/MathML"
xmlns:db="http://docbook.org/ns/docbook">
<refnamediv>
<refname>accsum_fasttwosum</refname><refpurpose>The fast2sum sum of a and b.</refpurpose>
</refnamediv>
<refsynopsisdiv>
<title>Calling Sequence</title>
<synopsis>
s = accsum_fasttwosum ( a , b )
[s,t] = accsum_fasttwosum ( a , b )
</synopsis>
</refsynopsisdiv>
<refsection>
<title>Parameters</title>
<variablelist>
<varlistentry><term>a :</term>
<listitem><para> a 1-by-1 matrix of doubles</para></listitem></varlistentry>
<varlistentry><term>b :</term>
<listitem><para> a 1-by-1 matrix of doubles</para></listitem></varlistentry>
<varlistentry><term>s :</term>
<listitem><para> a 1-by-1 matrix of doubles, the highest significant digits of the sum</para></listitem></varlistentry>
<varlistentry><term>t :</term>
<listitem><para> a 1-by-1, matrix of doubles, the lowest significant digits of the sum</para></listitem></varlistentry>
</variablelist>
</refsection>
<refsection>
<title>Description</title>
<para>
Returns the sum of a and b
so that a + b = s + t exactly.
where s contains the highest significant digits
and t the lowest significant digits.
</para>
<para>
Assumes that |a|&gt;=|b|. If not, the arguments are switched.
</para>
<para>
Algorithm 4.3 in HFCC
Due to Dekker (and Kahan)
Assumes that Scilab uses round-to-nearest.
</para>
<para>
</para>
</refsection>
<refsection>
<title>Examples</title>
<programlisting role="example"><![CDATA[
[s,t] = accsum_fasttwosum ( 2 , 1 ) // 3
// is 1+(%eps/2) but is 1 without algo
[s,t] = accsum_fasttwosum ( 1 , %eps/2 )
]]></programlisting>
</refsection>
<refsection>
<title>Authors</title>
<simplelist type="vert">
<member>Michael Baudin, 2010</member>
</simplelist>
</refsection>
<refsection>
<title>Bibliography</title>
<para>"Stability and numerical accuracy of algorithms", Nicolas Higham</para>
<para>"Handbook of Floating Point Computations", Muller et al</para>
</refsection>
</refentry>
help/en_US/accsum_fscs.xml
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
<?xml version="1.0" encoding="UTF-8"?>
<!--
*
* This help file was generated from accsum_fscs.sci using help_from_sci().
*
-->
<refentry version="5.0-subset Scilab" xml:id="accsum_fscs" xml:lang="en"
xmlns="http://docbook.org/ns/docbook"
xmlns:xlink="http://www.w3.org/1999/xlink"
xmlns:svg="http://www.w3.org/2000/svg"
xmlns:ns3="http://www.w3.org/1999/xhtml"
xmlns:mml="http://www.w3.org/1998/Math/MathML"
xmlns:db="http://docbook.org/ns/docbook">
<refnamediv>
<refname>accsum_fscs</refname><refpurpose>A Self Compensated Sum algorithm</refpurpose>
</refnamediv>
<refsynopsisdiv>
<title>Calling Sequence</title>
<synopsis>
s = accsum_fscs ( x )
</synopsis>
</refsynopsisdiv>
<refsection>
<title>Parameters</title>
<variablelist>
<varlistentry><term>x :</term>
<listitem><para> a m-by-n, matrix of doubles</para></listitem></varlistentry>
<varlistentry><term>s :</term>
<listitem><para> a 1-by-1 matrix of doubles, the highest significant digits of the sum</para></listitem></varlistentry>
</variablelist>
</refsection>
<refsection>
<title>Description</title>
<para>
A Self Compensated Sum algorithm.
Uses the fasttwosum algorithm.
This is a fast implementation, based on compiled source code.
</para>
<para>
</para>
</refsection>
<refsection>
<title>Examples</title>
<programlisting role="example"><![CDATA[
s = accsum_fscs ( [2 1] ) // 3
s = accsum_fscs ( [1 2] ) // 3
x = accsum_wilkinson(10); size(x,"*")
s = sum(x)
s = accsum_fscs ( x )
]]></programlisting>
</refsection>
<refsection>
<title>Authors</title>
<simplelist type="vert">
<member>Michael Baudin, 2010-2011</member>
</simplelist>
</refsection>
<refsection>
<title>Bibliography</title>
<para>"Stability and numerical accuracy of algorithms", Nicolas Higham</para>
<para>"Handbook of Floating Point Computations", Muller et al</para>
<para>https://hpcrd.lbl.gov/SCG/ocean/NRS/ECMWF/img14.htm</para>
<para>https://hpcrd.lbl.gov/SCG/ocean/NRS/SCSsum.F</para>
</refsection>
</refentry>
<?xml version="1.0" encoding="UTF-8"?>
<!--
*
* This help file was generated from accsum_fscs.sci using help_from_sci().
*
-->
<refentry version="5.0-subset Scilab" xml:id="accsum_fscs" xml:lang="en"
xmlns="http://docbook.org/ns/docbook"
xmlns:xlink="http://www.w3.org/1999/xlink"
xmlns:svg="http://www.w3.org/2000/svg"
xmlns:ns3="http://www.w3.org/1999/xhtml"
xmlns:mml="http://www.w3.org/1998/Math/MathML"
xmlns:db="http://docbook.org/ns/docbook">
<refnamediv>
<refname>accsum_fscs</refname><refpurpose>A Self Compensated Sum algorithm</refpurpose>
</refnamediv>
<refsynopsisdiv>
<title>Calling Sequence</title>
<synopsis>
s = accsum_fscs ( x )
[s,e] = accsum_fscs ( x )
</synopsis>
</refsynopsisdiv>
<refsection>
<title>Parameters</title>
<variablelist>
<varlistentry><term>x :</term>
<listitem><para> a m-by-n, matrix of doubles</para></listitem></varlistentry>
<varlistentry><term>s :</term>
<listitem><para> a 1-by-1 matrix of doubles, the highest significant digits of the sum</para></listitem></varlistentry>
<varlistentry><term>e :</term>
<listitem><para> a 1-by-1, matrix of doubles, the lowest significant digits of the sum</para></listitem></varlistentry>
</variablelist>
</refsection>
<refsection>
<title>Description</title>
<para>
A Self Compensated Sum algorithm.
Uses the fasttwosum algorithm.
This is a fast implementation, based on compiled source code.
</para>
<para>
</para>
</refsection>
<refsection>
<title>Examples</title>
<programlisting role="example"><![CDATA[
s = accsum_fscs ( [2 1] ) // 3
s = accsum_fscs ( [1 2] ) // 3
x = accsum_wilkinson(10); size(x,"*")
s = sum(x)
s = accsum_fscs ( x )
]]></programlisting>
</refsection>
<refsection>
<title>Authors</title>
<simplelist type="vert">
<member>Michael Baudin, 2010-2011</member>
</simplelist>
</refsection>
<refsection>
<title>Bibliography</title>
<para>"Stability and numerical accuracy of algorithms", Nicolas Higham</para>
<para>"Handbook of Floating Point Computations", Muller et al</para>
<para>https://hpcrd.lbl.gov/SCG/ocean/NRS/ECMWF/img14.htm</para>
<para>https://hpcrd.lbl.gov/SCG/ocean/NRS/SCSsum.F</para>
</refsection>
</refentry>
help/en_US/accsum_twosum.xml
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
<?xml version="1.0" encoding="UTF-8"?>
<!--
*
* This help file was generated from accsum_twosum.sci using help_from_sci().
*
-->
<refentry version="5.0-subset Scilab" xml:id="accsum_twosum" xml:lang="en"
xmlns="http://docbook.org/ns/docbook"
xmlns:xlink="http://www.w3.org/1999/xlink"
xmlns:svg="http://www.w3.org/2000/svg"
xmlns:ns3="http://www.w3.org/1999/xhtml"
xmlns:mml="http://www.w3.org/1998/Math/MathML"
xmlns:db="http://docbook.org/ns/docbook">
<refnamediv>
<refname>accsum_twosum</refname><refpurpose>The twosum sum of a and b.</refpurpose>
</refnamediv>
<refsynopsisdiv>
<title>Calling Sequence</title>
<synopsis>
[s,t] = accsum_twosum ( a , b )
</synopsis>
</refsynopsisdiv>
<refsection>
<title>Parameters</title>
<variablelist>
<varlistentry><term>a :</term>
<listitem><para> a 1-by-1 matrix of doubles</para></listitem></varlistentry>
<varlistentry><term>b :</term>
<listitem><para> a 1-by-1 matrix of doubles</para></listitem></varlistentry>
<varlistentry><term>s :</term>
<listitem><para> a 1-by-1 matrix of doubles, the highest significant digits of the sum</para></listitem></varlistentry>
<varlistentry><term>t :</term>
<listitem><para> a 1-by-1, matrix of doubles, the lowest significant digits of the sum</para></listitem></varlistentry>
</variablelist>
</refsection>
<refsection>
<title>Description</title>
<para>
Returns the sum of a and b
so that a + b = s + t exactly.
where s contains the highest significant digits
and t the lowest significant digits.
Algorithm 4.4 in HFCC
Due to Knuth
No assumption on a , b
Assumes that Scilab uses round-to-nearest.
</para>
<para>
</para>
</refsection>
<refsection>
<title>Examples</title>
<programlisting role="example"><![CDATA[
[s,t] = accsum_twosum ( 2 , 1 ) // 3
[s,t] = accsum_twosum ( 1 , 2 ) // 3
// is 1+(%eps/2) but is 1 without algo
[s,t] = accsum_twosum ( 1 , %eps/2 )
// is 1+(%eps/2) but is 1 without algo
[s,t] = accsum_twosum ( %eps/2 , 1 )
]]></programlisting>
</refsection>
<refsection>
<title>Authors</title>
<simplelist type="vert">
<member>Michael Baudin, 2010</member>
</simplelist>
</refsection>
<refsection>
<title>Bibliography</title>
<para>"Stability and numerical accuracy of algorithms", Nicolas Higham</para>
<para>"Handbook of Floating Point Computations", Muller et al</para>
</refsection>
</refentry>
<?xml version="1.0" encoding="UTF-8"?>
<!--
*
* This help file was generated from accsum_twosum.sci using help_from_sci().
*
-->
<refentry version="5.0-subset Scilab" xml:id="accsum_twosum" xml:lang="en"
xmlns="http://docbook.org/ns/docbook"
xmlns:xlink="http://www.w3.org/1999/xlink"
xmlns:svg="http://www.w3.org/2000/svg"
xmlns:ns3="http://www.w3.org/1999/xhtml"
xmlns:mml="http://www.w3.org/1998/Math/MathML"
xmlns:db="http://docbook.org/ns/docbook">
<refnamediv>
<refname>accsum_twosum</refname><refpurpose>The twosum sum of a and b.</refpurpose>
</refnamediv>
<refsynopsisdiv>
<title>Calling Sequence</title>
<synopsis>
s = accsum_twosum ( a , b )
[s,t] = accsum_twosum ( a , b )
</synopsis>
</refsynopsisdiv>
<refsection>
<title>Parameters</title>
<variablelist>
<varlistentry><term>a :</term>
<listitem><para> a 1-by-1 matrix of doubles</para></listitem></varlistentry>
<varlistentry><term>b :</term>
<listitem><para> a 1-by-1 matrix of doubles</para></listitem></varlistentry>
<varlistentry><term>s :</term>
<listitem><para> a 1-by-1 matrix of doubles, the highest significant digits of the sum</para></listitem></varlistentry>
<varlistentry><term>t :</term>
<listitem><para> a 1-by-1, matrix of doubles, the lowest significant digits of the sum</para></listitem></varlistentry>
</variablelist>
</refsection>
<refsection>
<title>Description</title>
<para>
Returns the sum of a and b
so that a + b = s + t exactly.
where s contains the highest significant digits
and t the lowest significant digits.
Algorithm 4.4 in HFCC
Due to Knuth
No assumption on a , b
Assumes that Scilab uses round-to-nearest.
</para>
<para>
</para>
</refsection>
<refsection>
<title>Examples</title>
<programlisting role="example"><![CDATA[
[s,t] = accsum_twosum ( 2 , 1 ) // 3
[s,t] = accsum_twosum ( 1 , 2 ) // 3
// is 1+(%eps/2) but is 1 without algo
[s,t] = accsum_twosum ( 1 , %eps/2 )
// is 1+(%eps/2) but is 1 without algo
[s,t] = accsum_twosum ( %eps/2 , 1 )
]]></programlisting>
</refsection>
<refsection>
<title>Authors</title>
<simplelist type="vert">
<member>Michael Baudin, 2010</member>
</simplelist>
</refsection>
<refsection>
<title>Bibliography</title>
<para>"Stability and numerical accuracy of algorithms", Nicolas Higham</para>
<para>"Handbook of Floating Point Computations", Muller et al</para>
</refsection>
</refentry>
help/en_US/accsum_fdcs.xml
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
<?xml version="1.0" encoding="UTF-8"?>
<!--
*
* This help file was generated from accsum_fdcs.sci using help_from_sci().
*
-->
<refentry version="5.0-subset Scilab" xml:id="accsum_fdcs" xml:lang="en"
xmlns="http://docbook.org/ns/docbook"
xmlns:xlink="http://www.w3.org/1999/xlink"
xmlns:svg="http://www.w3.org/2000/svg"
xmlns:ns3="http://www.w3.org/1999/xhtml"
xmlns:mml="http://www.w3.org/1998/Math/MathML"
xmlns:db="http://docbook.org/ns/docbook">
<refnamediv>
<refname>accsum_fdcs</refname><refpurpose>A Doubly Self Compensated Sum algorithm</refpurpose>
</refnamediv>
<refsynopsisdiv>
<title>Calling Sequence</title>
<synopsis>
s = accsum_fdcs ( x )
</synopsis>
</refsynopsisdiv>
<refsection>
<title>Parameters</title>
<variablelist>
<varlistentry><term>x :</term>
<listitem><para> a m-by-n, matrix of doubles</para></listitem></varlistentry>
<varlistentry><term>s :</term>
<listitem><para> a 1-by-1 matrix of doubles, the highest significant digits of the sum</para></listitem></varlistentry>
</variablelist>
</refsection>
<refsection>
<title>Description</title>
<para>
A Doubly Self Compensated Sum algorithm.
Uses the fasttwosum algorithm.
This is a fast implementation, based on compiled source code.
</para>
<para>
</para>
</refsection>
<refsection>
<title>Examples</title>
<programlisting role="example"><![CDATA[
s = accsum_fdcs ( [2 1] ) // 3
s = accsum_fdcs ( [1 2] ) // 3
s = accsum_fdcs ( [2 1] ) // 3
s = accsum_fdcs ( [1 2] ) // 3
x = accsum_wilkinson(10); size(x,"*")
s = sum(x)
s = accsum_fdcs ( x )
s = accsum_fdcs ( x )
]]></programlisting>
</refsection>
<refsection>
<title>Authors</title>
<simplelist type="vert">
<member>Michael Baudin, 2010-2011</member>
</simplelist>
</refsection>
<refsection>
<title>Bibliography</title>
<para>"Stability and numerical accuracy of algorithms", Nicolas Higham</para>
<para>"Handbook of Floating Point Computations", Muller et al</para>
<para>https://hpcrd.lbl.gov/SCG/ocean/NRS/ECMWF/img14.htm</para>
<para>https://hpcrd.lbl.gov/SCG/ocean/NRS/SCSsum.F</para>
</refsection>
</refentry>
<?xml version="1.0" encoding="UTF-8"?>
<!--
*
* This help file was generated from accsum_fdcs.sci using help_from_sci().
*
-->
<refentry version="5.0-subset Scilab" xml:id="accsum_fdcs" xml:lang="en"
xmlns="http://docbook.org/ns/docbook"
xmlns:xlink="http://www.w3.org/1999/xlink"
xmlns:svg="http://www.w3.org/2000/svg"
xmlns:ns3="http://www.w3.org/1999/xhtml"
xmlns:mml="http://www.w3.org/1998/Math/MathML"
xmlns:db="http://docbook.org/ns/docbook">
<refnamediv>
<refname>accsum_fdcs</refname><refpurpose>A Doubly Self Compensated Sum algorithm</refpurpose>
</refnamediv>
<refsynopsisdiv>
<title>Calling Sequence</title>
<synopsis>
s = accsum_fdcs ( x )
[s,e] = accsum_fdcs ( x )
</synopsis>
</refsynopsisdiv>
<refsection>
<title>Parameters</title>
<variablelist>
<varlistentry><term>x :</term>
<listitem><para> a m-by-n, matrix of doubles</para></listitem></varlistentry>
<varlistentry><term>s :</term>
<listitem><para> a 1-by-1 matrix of doubles, the highest significant digits of the sum</para></listitem></varlistentry>
<varlistentry><term>e :</term>
<listitem><para> a 1-by-1, matrix of doubles, the lowest significant digits of the sum</para></listitem></varlistentry>
</variablelist>
</refsection>
<refsection>
<title>Description</title>
<para>
A Doubly Self Compensated Sum algorithm.
Uses the fasttwosum algorithm.
This is a fast implementation, based on compiled source code.
</para>
<para>
</para>
</refsection>
<refsection>
<title>Examples</title>
<programlisting role="example"><![CDATA[
s = accsum_fdcs ( [2 1] ) // 3
s = accsum_fdcs ( [1 2] ) // 3
s = accsum_fdcs ( [2 1] ) // 3
s = accsum_fdcs ( [1 2] ) // 3
x = accsum_wilkinson(10); size(x,"*")
s = sum(x)
s = accsum_fdcs ( x )
s = accsum_fdcs ( x )
]]></programlisting>
</refsection>
<refsection>
<title>Authors</title>
<simplelist type="vert">
<member>Michael Baudin, 2010-2011</member>
</simplelist>
</refsection>
<refsection>
<title>Bibliography</title>
<para>"Stability and numerical accuracy of algorithms", Nicolas Higham</para>
<para>"Handbook of Floating Point Computations", Muller et al</para>
<para>https://hpcrd.lbl.gov/SCG/ocean/NRS/ECMWF/img14.htm</para>
<para>https://hpcrd.lbl.gov/SCG/ocean/NRS/SCSsum.F</para>
</refsection>
</refentry>
help/en_US/support/accsum_priestx.xml
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
<?xml version="1.0" encoding="UTF-8"?>
<!--
*
* This help file was generated from accsum_priestx.sci using help_from_sci().
*
-->
<refentry version="5.0-subset Scilab" xml:id="accsum_priestx" xml:lang="en"
xmlns="http://docbook.org/ns/docbook"
xmlns:xlink="http://www.w3.org/1999/xlink"
xmlns:svg="http://www.w3.org/2000/svg"
xmlns:ns3="http://www.w3.org/1999/xhtml"
xmlns:mml="http://www.w3.org/1998/Math/MathML"
xmlns:db="http://docbook.org/ns/docbook">
<refnamediv>
<refname>accsum_priestx</refname><refpurpose>A difficult example for SCS by Priest.</refpurpose>
</refnamediv>
<refsynopsisdiv>
<title>Calling Sequence</title>
<synopsis>
x = accsum_priestx ( )
</synopsis>
</refsynopsisdiv>
<refsection>
<title>Parameters</title>
<variablelist>
<varlistentry><term>x :</term>
<listitem><para> a 6-by-1, matrix of doubles</para></listitem></varlistentry>
</variablelist>
</refsection>
<refsection>
<title>Description</title>
<para>
"On properties of floating point arithmetics:
numerical stability and the cost of accurate computations"
Douglas Priest, 1992
p62
"What is the smallest n for which simply compensated
summation is not guaranteed to yield a result
with a small forward relative error in an arithmetic which
conforms to the IEEE 754 standard. The smallest known
to the author is ... Then the exact sum is 2, but the sum
computed by the SCS is 3."
</para>
<para>
</para>
</refsection>
<refsection>
<title>Examples</title>
<programlisting role="example"><![CDATA[
x = accsum_priestx ( )
[s,e] = accsum_SCS ( x ) // = 3 + 0
[s,e] = accsum_DCS ( x ) // = 2 + 0
]]></programlisting>
</refsection>
</refentry>
help/en_US/support/accsum_orderdynamic.xml
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
<?xml version="1.0" encoding="UTF-8"?>
<!--
*
* This help file was generated from accsum_orderdynamic.sci using help_from_sci().
*
-->
<refentry version="5.0-subset Scilab" xml:id="accsum_orderdynamic" xml:lang="en"
xmlns="http://docbook.org/ns/docbook"
xmlns:xlink="http://www.w3.org/1999/xlink"
xmlns:svg="http://www.w3.org/2000/svg"
xmlns:ns3="http://www.w3.org/1999/xhtml"
xmlns:mml="http://www.w3.org/1998/Math/MathML"
xmlns:db="http://docbook.org/ns/docbook">
<refnamediv>
<refname>accsum_orderdynamic</refname><refpurpose>Try to find an order which makes the sum perform bad/good.</refpurpose>
</refnamediv>
<refsynopsisdiv>
<title>Calling Sequence</title>
<synopsis>
s = accsum_orderdynamic ( x , order )
</synopsis>
</refsynopsisdiv>
<refsection>
<title>Parameters</title>
<variablelist>
<varlistentry><term>x :</term>
<listitem><para> a m-by-n, matrix of doubles</para></listitem></varlistentry>
<varlistentry><term>order :</term>
<listitem><para> a 1-by-1 matrix of doubles, a positive integer, the order. order = 1 : search for the entry which makes the smallest |s+y(i)|, good for accuracy, order = 2 : search for the entry which makes the largest |s+y(i)|, bad for accuracy</para></listitem></varlistentry>
<varlistentry><term>s :</term>
<listitem><para> a 1-by-1 matrix of doubles, the sum</para></listitem></varlistentry>
</variablelist>
</refsection>
<refsection>
<title>Description</title>
<para>
Try to find an order which makes the sum perform bad/good.
</para>
<para>
</para>
</refsection>
<refsection>
<title>Examples</title>
<programlisting role="example"><![CDATA[
// Test on simple data
x = [-9.9 5.5 2.2 -3.3 -6.6 0 1.1]
sum(x)
s = accsum_orderdynamic ( x , 1 )
s = accsum_orderdynamic ( x , 2 )
// Test on difficult data
path = accsum_getpath ( );
filename=fullfile(path,"demos","etaana.dat");
x=fscanfMat(filename);
// 7.390949249267578125000
sum(x)
// 0.3437500000000000000000
s = accsum_orderdynamic ( x , 1 )
// - 73.015625
s = accsum_orderdynamic ( x , 2 )
]]></programlisting>
</refsection>
<refsection>
<title>Authors</title>
<simplelist type="vert">
<member>Michael Baudin, 2010</member>
</simplelist>
</refsection>
<refsection>
<title>Bibliography</title>
<para>"Stability and numerical accuracy of algorithms", Nicolas Higham</para>
<para>"Handbook of Floating Point Computations", Muller et al</para>
<para>https://hpcrd.lbl.gov/SCG/ocean/NRS/ECMWF/img16.htm</para>
</refsection>
</refentry>
help/en_US/accsum_compsum.xml
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
<?xml version="1.0" encoding="UTF-8"?>
<!--
*
* This help file was generated from accsum_compsum.sci using help_from_sci().
*
-->
<refentry version="5.0-subset Scilab" xml:id="accsum_compsum" xml:lang="en"
xmlns="http://docbook.org/ns/docbook"
xmlns:xlink="http://www.w3.org/1999/xlink"
xmlns:svg="http://www.w3.org/2000/svg"
xmlns:ns3="http://www.w3.org/1999/xhtml"
xmlns:mml="http://www.w3.org/1998/Math/MathML"
xmlns:db="http://docbook.org/ns/docbook">
<refnamediv>
<refname>accsum_compsum</refname><refpurpose>The compensated sum of a matrix.</refpurpose>
</refnamediv>
<refsynopsisdiv>
<title>Calling Sequence</title>
<synopsis>
[s,e] = accsum_compsum ( x )
</synopsis>
</refsynopsisdiv>
<refsection>
<title>Parameters</title>
<variablelist>
<varlistentry><term>x :</term>
<listitem><para> a n-by-m matrix of doubles</para></listitem></varlistentry>
<varlistentry><term>s :</term>
<listitem><para> a 1-by-1 matrix of doubles, the highest significant digits of the sum</para></listitem></varlistentry>
<varlistentry><term>e :</term>
<listitem><para> a 1-by-1, matrix of doubles, the lowest significant digits of the sum</para></listitem></varlistentry>
</variablelist>
</refsection>
<refsection>
<title>Description</title>
<para>
Returns the sum of x.
Algorithm 4.4 in HFCC
Due to Knuth
No assumption on a , b
Assumes that Scilab uses round-to-nearest.
</para>
<para>
</para>
</refsection>
<refsection>
<title>Examples</title>
<programlisting role="example"><![CDATA[
[s,e] = accsum_compsum ( [2 1] ) // 3
[s,e] = accsum_compsum ( [1 2] ) // 3
// is 1+(%eps/2) but is 1 without algo
[s,e] = accsum_compsum ( [1 %eps/2 ] )
// is 1+(%eps/2) but is 1 without algo
[s,e] = accsum_compsum ( [%eps/2 1] )
]]></programlisting>
</refsection>
<refsection>
<title>Authors</title>
<simplelist type="vert">
<member>Michael Baudin, 2010</member>
</simplelist>
</refsection>
<refsection>
<title>Bibliography</title>
<para>"Stability and numerical accuracy of algorithms", Nicolas Higham</para>
<para>"Handbook of Floating Point Computations", Muller et al</para>
</refsection>
</refentry>
<?xml version="1.0" encoding="UTF-8"?>
<!--
*
* This help file was generated from accsum_compsum.sci using help_from_sci().
*
-->
<refentry version="5.0-subset Scilab" xml:id="accsum_compsum" xml:lang="en"
xmlns="http://docbook.org/ns/docbook"
xmlns:xlink="http://www.w3.org/1999/xlink"
xmlns:svg="http://www.w3.org/2000/svg"
xmlns:ns3="http://www.w3.org/1999/xhtml"
xmlns:mml="http://www.w3.org/1998/Math/MathML"
xmlns:db="http://docbook.org/ns/docbook">
<refnamediv>
<refname>accsum_compsum</refname><refpurpose>The compensated sum of a matrix.</refpurpose>
</refnamediv>
<refsynopsisdiv>
<title>Calling Sequence</title>
<synopsis>
s = accsum_compsum ( x )
[s,e] = accsum_compsum ( x )
</synopsis>
</refsynopsisdiv>
<refsection>
<title>Parameters</title>
<variablelist>
<varlistentry><term>x :</term>
<listitem><para> a n-by-m matrix of doubles</para></listitem></varlistentry>
<varlistentry><term>s :</term>
<listitem><para> a 1-by-1 matrix of doubles, the highest significant digits of the sum</para></listitem></varlistentry>
<varlistentry><term>e :</term>
<listitem><para> a 1-by-1, matrix of doubles, the lowest significant digits of the sum</para></listitem></varlistentry>
</variablelist>
</refsection>
<refsection>
<title>Description</title>
<para>
Returns s and e such that s+e is the sum of x.
Algorithm 4.4 in HFCC
Due to Knuth
No assumption on a , b
Assumes that Scilab uses round-to-nearest.
</para>
<para>
</para>
</refsection>
<refsection>
<title>Examples</title>
<programlisting role="example"><![CDATA[
[s,e] = accsum_compsum ( [2 1] ) // 3
[s,e] = accsum_compsum ( [1 2] ) // 3
// is 1+(%eps/2) but is 1 without algo
[s,e] = accsum_compsum ( [1 %eps/2 ] )
// is 1+(%eps/2) but is 1 without algo
[s,e] = accsum_compsum ( [%eps/2 1] )
]]></programlisting>
</refsection>
<refsection>
<title>Authors</title>
<simplelist type="vert">
<member>Michael Baudin, 2010</member>
</simplelist>
</refsection>
<refsection>
<title>Bibliography</title>
<para>"Stability and numerical accuracy of algorithms", Nicolas Higham</para>
<para>"Handbook of Floating Point Computations", Muller et al</para>
</refsection>
</refentry>
help/en_US/accsum_dblcompsum.xml
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
<?xml version="1.0" encoding="UTF-8"?>
<!--
*
* This help file was generated from accsum_dblcompsum.sci using help_from_sci().
*
-->
<refentry version="5.0-subset Scilab" xml:id="accsum_dblcompsum" xml:lang="en"
xmlns="http://docbook.org/ns/docbook"
xmlns:xlink="http://www.w3.org/1999/xlink"
xmlns:svg="http://www.w3.org/2000/svg"
xmlns:ns3="http://www.w3.org/1999/xhtml"
xmlns:mml="http://www.w3.org/1998/Math/MathML"
xmlns:db="http://docbook.org/ns/docbook">
<refnamediv>
<refname>accsum_dblcompsum</refname><refpurpose>The doubly compensated sum of a matrix.</refpurpose>
</refnamediv>
<refsynopsisdiv>
<title>Calling Sequence</title>
<synopsis>
[s,e] = accsum_dblcompsum ( x )
</synopsis>
</refsynopsisdiv>
<refsection>
<title>Parameters</title>
<variablelist>
<varlistentry><term>x :</term>
<listitem><para> a n-by-m matrix of doubles</para></listitem></varlistentry>
<varlistentry><term>s :</term>
<listitem><para> a 1-by-1 matrix of doubles, the highest significant digits of the sum</para></listitem></varlistentry>
<varlistentry><term>e :</term>
<listitem><para> a 1-by-1, matrix of doubles, the lowest significant digits of the sum</para></listitem></varlistentry>
</variablelist>
</refsection>
<refsection>
<title>Description</title>
<para>
Returns the sum of x.
Algorithm 4.2 in SNAA (Doubly Compensated Summation)
Due to Priest (and Kahan)
Assumes that Scilab uses round-to-nearest.
</para>
<para>
</para>
</refsection>
<refsection>
<title>Examples</title>
<programlisting role="example"><![CDATA[
[s,e] = accsum_dblcompsum ( [2 1] ) // 3
[s,e] = accsum_dblcompsum ( [1 2] ) // 3
// is 1+(%eps/2) but is 1 without algo
[s,e] = accsum_dblcompsum ( [1 %eps/2 ] )
// is 1+(%eps/2) but is 1 without algo
[s,e] = accsum_dblcompsum ( [%eps/2 1] )
]]></programlisting>
</refsection>
<refsection>
<title>Authors</title>
<simplelist type="vert">
<member>Michael Baudin, 2010</member>
</simplelist>
</refsection>
<refsection>
<title>Bibliography</title>
<para>"Stability and numerical accuracy of algorithms", Nicolas Higham</para>
<para>"Handbook of Floating Point Computations", Muller et al</para>
</refsection>
</refentry>
<?xml version="1.0" encoding="UTF-8"?>
<!--
*
* This help file was generated from accsum_dblcompsum.sci using help_from_sci().
*
-->
<refentry version="5.0-subset Scilab" xml:id="accsum_dblcompsum" xml:lang="en"
xmlns="http://docbook.org/ns/docbook"
xmlns:xlink="http://www.w3.org/1999/xlink"
xmlns:svg="http://www.w3.org/2000/svg"
xmlns:ns3="http://www.w3.org/1999/xhtml"
xmlns:mml="http://www.w3.org/1998/Math/MathML"
xmlns:db="http://docbook.org/ns/docbook">
<refnamediv>
<refname>accsum_dblcompsum</refname><refpurpose>The doubly compensated sum of a matrix.</refpurpose>
</refnamediv>
<refsynopsisdiv>
<title>Calling Sequence</title>
<synopsis>
s = accsum_dblcompsum ( x )
[s,e] = accsum_dblcompsum ( x )
</synopsis>
</refsynopsisdiv>
<refsection>
<title>Parameters</title>
<variablelist>
<varlistentry><term>x :</term>
<listitem><para> a n-by-m matrix of doubles</para></listitem></varlistentry>
<varlistentry><term>s :</term>
<listitem><para> a 1-by-1 matrix of doubles, the highest significant digits of the sum</para></listitem></varlistentry>
<varlistentry><term>e :</term>
<listitem><para> a 1-by-1, matrix of doubles, the lowest significant digits of the sum</para></listitem></varlistentry>
</variablelist>
</refsection>
<refsection>
<title>Description</title>
<para>
Returns the sum of x.
Algorithm 4.2 in SNAA (Doubly Compensated Summation)
Due to Priest (and Kahan)
Assumes that Scilab uses round-to-nearest.
</para>
<para>
</para>
</refsection>
<refsection>
<title>Examples</title>
<programlisting role="example"><![CDATA[
[s,e] = accsum_dblcompsum ( [2 1] ) // 3
[s,e] = accsum_dblcompsum ( [1 2] ) // 3
// is 1+(%eps/2) but is 1 without algo
[s,e] = accsum_dblcompsum ( [1 %eps/2 ] )
// is 1+(%eps/2) but is 1 without algo
[s,e] = accsum_dblcompsum ( [%eps/2 1] )
]]></programlisting>
</refsection>
<refsection>
<title>Authors</title>
<simplelist type="vert">
<member>Michael Baudin, 2010</member>
</simplelist>
</refsection>
<refsection>
<title>Bibliography</title>
<para>"Stability and numerical accuracy of algorithms", Nicolas Higham</para>
<para>"Handbook of Floating Point Computations", Muller et al</para>
</refsection>
</refentry>
sci_gateway/c/sci_accsum_fcompsum.c
1
1
2
23
34
45
......
1819
1920
2021
21
22
23
24
2225
2326
2427
......
2730
2831
2932
30
3133
34
35
3236
33
37
3438
3539
3640
......
6872
6973
7074
71
72
73
74
75
76
77
78
79
8075
8176
82
83
84
85
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
86104
87
88105
89106
// Copyright (C) 2010 - Michael Baudin
// Copyright (C) 2010-2011 - Michael Baudin
// Copyright (C) 2011 - DIGITEO - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
#include "accsum.h"
// s = accsum_fcompsum(x) returns the sum of x, with a compensated summation.
// s = accsum_fcompsum(x)
// [s,e] = accsum_fcompsum(x)
// returns the sum of x, with a compensated summation.
int sci_accsum_fcompsum (char * fname) {
SciErr sciErr;
int rowsX,colsX;
double* lX = NULL;
int iComplex = 0;
double* lS = NULL;
int i;
double s;
double e;
int minlhs=1, minrhs=1, maxlhs=1, maxrhs=1;
int minlhs=1, minrhs=1, maxlhs=2, maxrhs=1;
CheckRhs(minrhs,maxrhs) ;
CheckLhs(minlhs,maxlhs) ;
return 0;
}
//
// Create s
//
sciErr = allocMatrixOfDouble(pvApiCtx,Rhs+1, 1, 1, &lS);
if(sciErr.iErr)
{
printError(&sciErr, 0);
return 0;
}
//
// Perform the sum
//
accsum_fcompsum(lX, rowsX*colsX, lS);
LhsVar(1) = Rhs+1;
accsum_fcompsum(lX, rowsX*colsX, &s, &e);
//
// Set LHS
if ( Lhs>=1 )
{
//
// Create s
sciErr = createMatrixOfDouble(pvApiCtx,Rhs+1, 1, 1, &s);
if(sciErr.iErr)
{
printError(&sciErr, 0);
return 0;
}
LhsVar(1) = Rhs+1;
}
if ( Lhs==2 )
{
//
// Create e
sciErr = createMatrixOfDouble(pvApiCtx,Rhs+2, 1, 1, &e);
if(sciErr.iErr)
{
printError(&sciErr, 0);
return 0;
}
LhsVar(2) = Rhs+2;
}
return(0);
}
sci_gateway/c/sci_accsum_fscs.c
1
1
2
23
34
45
......
1819
1920
2021
21
22
23
24
2225
2326
2427
......
2730
2831
2932
30
3133
34
35
3236
33
37
3438
3539
3640
......
6872
6973
7074
71
72
73
74
75
76
77
78
79
8075
81
82
76
8377
84
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
85104
86105
87106
// Copyright (C) 2010 - Michael Baudin
// Copyright (C) 2010-2011 - Michael Baudin
// Copyright (C) 2011 - DIGITEO - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
#include "accsum.h"
// s = accsum_fscs(x) returns the sum of x, with a simply compensated summation.
// s = accsum_fdcs(x)
// [s,e] = accsum_fdcs(x)
// returns the sum of x, with a simply compensated summation.
int sci_accsum_fscs (char * fname) {
SciErr sciErr;
int rowsX,colsX;
double* lX = NULL;
int iComplex = 0;
double* lS = NULL;
int i;
double s;
double e;
int minlhs=1, minrhs=1, maxlhs=1, maxrhs=1;
int minlhs=1, minrhs=1, maxlhs=2, maxrhs=1;
CheckRhs(minrhs,maxrhs) ;
CheckLhs(minlhs,maxlhs) ;
return 0;
}
//
// Create s
//
sciErr = allocMatrixOfDouble(pvApiCtx,Rhs+1, 1, 1, &lS);
if(sciErr.iErr)
{
printError(&sciErr, 0);
return 0;
}
//
// Perform the sum
//
accsum_fscs(lX, rowsX*colsX, lS);
accsum_fscs(lX, rowsX*colsX, &s, &e);
LhsVar(1) = Rhs+1;
//
// Set LHS
if ( Lhs>=1 )
{
//
// Create s
sciErr = createMatrixOfDouble(pvApiCtx,Rhs+1, 1, 1, &s);
if(sciErr.iErr)
{
printError(&sciErr, 0);
return 0;
}
LhsVar(1) = Rhs+1;
}
if ( Lhs==2 )
{
//
// Create e
sciErr = createMatrixOfDouble(pvApiCtx,Rhs+2, 1, 1, &e);
if(sciErr.iErr)
{
printError(&sciErr, 0);
return 0;
}
LhsVar(2) = Rhs+2;
}
return(0);
sci_gateway/c/sci_accsum_fdcs.c
1
1
2
23
34
45
......
1819
1920
2021
21
22
23
24
2225
2326
2427
......
2730
2831
2932
30
3133
34
35
3236
33
37
3438
3539
3640
......
6872
6973
7074
71
72
73
74
75
76
77
78
79
8075
81
82
76
8377
84
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
85104
86105
87106
// Copyright (C) 2010 - Michael Baudin
// Copyright (C) 2010-2011 - Michael Baudin
// Copyright (C) 2011 - DIGITEO - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
#include "accsum.h"
// s = accsum_fdcs(x) returns the sum of x, with a doubly compensated summation.
// s = accsum_fdcs(x)
// [s,e] = accsum_fdcs(x)
// returns the sum of x, with a doubly compensated summation.
int sci_accsum_fdcs (char * fname) {
SciErr sciErr;
int rowsX,colsX;
double* lX = NULL;
int iComplex = 0;
double* lS = NULL;
int i;
double s;
double e;
int minlhs=1, minrhs=1, maxlhs=1, maxrhs=1;
int minlhs=1, minrhs=1, maxlhs=2, maxrhs=1;
CheckRhs(minrhs,maxrhs) ;
CheckLhs(minlhs,maxlhs) ;
return 0;
}
//
// Create s
//
sciErr = allocMatrixOfDouble(pvApiCtx,Rhs+1, 1, 1, &lS);
if(sciErr.iErr)
{
printError(&sciErr, 0);
return 0;
}
//
// Perform the sum
//
accsum_fdcs(lX, rowsX*colsX, lS);
accsum_fdcs(lX, rowsX*colsX, &s, &e);
LhsVar(1) = Rhs+1;
//
// Set LHS
if ( Lhs>=1 )
{
//
// Create s
sciErr = createMatrixOfDouble(pvApiCtx,Rhs+1, 1, 1, &s);
if(sciErr.iErr)
{
printError(&sciErr, 0);
return 0;
}
LhsVar(1) = Rhs+1;
}
if ( Lhs==2 )
{
//
// Create e
sciErr = createMatrixOfDouble(pvApiCtx,Rhs+2, 1, 1, &e);
if(sciErr.iErr)
{
printError(&sciErr, 0);
return 0;
}
LhsVar(2) = Rhs+2;
}
return(0);

Archive Download the corresponding diff file

Revision: 15