accsum

accsum Commit Details

Date:2015-09-02 15:30:19 (2 years 10 months ago)
Author:John Gliksberg
Commit:33
Parents: 32
Message:more cosmetics in src/ and macros/
Changes:
M/macros/accsum_orderdynamic.sci
M/macros/accsum_compsum.sci
M/macros/accsum_straight.sci
M/macros/accsum_scs.sci
M/macros/accsum_getpath.sci
M/macros/accsum_shuffle.sci
M/macros/accsum_higham.sci
M/macros/accsum_priestx.sci
M/macros/accsum_twosum.sci
M/macros/accsum_dblcompsum.sci
M/macros/accsum_sumcond.sci
M/src/c/accsum.c
M/macros/accsum_dcs.sci
M/macros/accsum_order.sci
M/macros/accsum_fasttwosum.sci
M/macros/accsum_wilkinson.sci
M/src/c/builder_c.sce
M/src/c/accsum.h

File differences

macros/accsum_fasttwosum.sci
77
88
99
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
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
4648
47
48
49
50
51
52
53
54
55
56
57
58
59
49
50
51
52
53
54
55
56
57
58
59
60
61
62
6063
6164
62
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
function [s,t] = accsum_fasttwosum ( a , b )
// The fast2sum sum of a and b.
//
// Calling Sequence
// s = accsum_fasttwosum ( a , b )
// [s,t] = accsum_fasttwosum ( 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.
//
// Assumes that |a|>=|b|. If not, the arguments are switched.
//
// Algorithm 4.3 in HFCC
// Due to Dekker (and Kahan)
// Assumes that Scilab uses round-to-nearest.
//
// Examples
// [s,t] = accsum_fasttwosum ( 2 , 1 ) // 3
// // is 1+(%eps/2) but is 1 without algo
// [s,t] = accsum_fasttwosum ( 1 , %eps/2 )
//
// Authors
// Michael Baudin, 2010
//
// Bibliography
// "Stability and numerical accuracy of algorithms", Nicolas Higham
// "Handbook of Floating Point Computations", Muller et al
function [s, t] = accsum_fasttwosum(a, b)
// The fast2sum sum of a and b.
//
// Calling Sequence
// s = accsum_fasttwosum ( a , b )
// [s,t] = accsum_fasttwosum ( 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.
//
// Enforces |a| >= |b|
//
// Algorithm 4.3 in HFCC
// Due to Dekker (and Kahan)
// Assumes that Scilab uses round-to-nearest.
//
// Examples
// [s,t] = accsum_fasttwosum ( 2 , 1 ) // 3
// // is 1+(%eps/2) but is 1 without algo
// [s,t] = accsum_fasttwosum ( 1 , %eps/2 )
//
// Authors
// Michael Baudin, 2010
//
// Bibliography
// "Stability and numerical accuracy of algorithms", Nicolas Higham
// "Handbook of Floating Point Computations", Muller et al
[lhs, rhs] = argn()
apifun_checkrhs ( "accsum_fasttwosum" , rhs , 2:2 )
apifun_checklhs ( "accsum_fasttwosum" , lhs , 0:2 )
//
if ( abs(a) < abs(b) ) then
// Switch a and b
tmp = a
a = b
b = tmp
end
s = a + b
z = s - a
t = b - z
[lhs, rhs] = argn();
apifun_checkrhs("accsum_fasttwosum", rhs, 2:2);
apifun_checklhs("accsum_fasttwosum", lhs, 0:2);
if abs(a) < abs(b) then
// Switch a and b
tmp = a;
a = b;
b = tmp;
end
s = a + b;
z = s - a;
t = b - z;
endfunction
macros/accsum_wilkinson.sci
77
88
99
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
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
4249
43
44
45
46
47
48
49
50
51
52
53
54
50
51
52
53
54
55
56
57
58
59
60
61
5562
5663
// 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
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
[lhs, rhs] = argn()
apifun_checkrhs ( "accsum_wilkinson" , rhs , 1:1 )
apifun_checklhs ( "accsum_wilkinson" , lhs , 0:1 )
//
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
[lhs, rhs] = argn();
apifun_checkrhs("accsum_wilkinson", rhs, 1:1);
apifun_checklhs("accsum_wilkinson", lhs, 0:1);
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
77
88
99
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
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
6165
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
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
8091
81
82
83
84
85
86
87
8892
8993
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
function s = accsum_orderdynamic ( x , order )
// Compute the sum with a dynamic re-ordering.
//
// 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.
// Finding the optimal ordering of a sum, which minimizes the
// intermediate sums, is a combinatorial problem which is too
// expensive to solve in general.
// A compromise is used here.
// For k from 1 to n,
// we search for i which minimizes (if order=1) or maximizes
// (if order=2) the term |s+x(i)|, among the entries from
// 1 to n-k+1.
// Then the entries #i and #n-k+1 are switched.
//
// This algorithm is called "Psum" in the Chapter 4 "Summation" of
// SNAA.
//
// 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-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/img16.htm
function s = accsum_orderdynamic(x, order)
// Compute the sum with a dynamic re-ordering.
//
// 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.
// Finding the optimal ordering of a sum, which minimizes the
// intermediate sums, is a combinatorial problem which is too
// expensive to solve in general.
// A compromise is used here.
// For k from 1 to n,
// we search for i which minimizes (if order=1) or maximizes
// (if order=2) the term |s+x(i)|, among the entries from
// 1 to n-k+1.
// Then the entries #i and #n-k+1 are switched.
//
// This algorithm is called "Psum" in the Chapter 4 "Summation" of
// SNAA.
//
// 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-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/img16.htm
[lhs, rhs] = argn()
apifun_checkrhs ( "accsum_orderdynamic" , rhs , 2:2 )
apifun_checklhs ( "accsum_orderdynamic" , lhs , 0:1 )
//
s = 0
z = x
n = size(x,"*")
// We remove an entry of z at each iteration of the loop:
// in the end, the array z is empty.
for k = 1 : n
// Sort into increasing order, from 1 to n-k+1.
[B,i] = gsort(abs(s+z(1:n-k+1)),"g","i")
if ( order == 1 ) then
// Set in i the indice which makes |s+z(i)| minimum
i = i(1)
else
// Set in i the indice which makes |s+z(i)| maximum
i = i($)
[lhs, rhs] = argn();
apifun_checkrhs("accsum_orderdynamic", rhs, 2:2);
apifun_checklhs("accsum_orderdynamic", lhs, 0:1);
s = 0;
z = x;
n = size(x, "*");
// We remove an entry of z at each iteration of the loop:
// in the end, the array z is empty.
for k = 1 : n
// Sort into increasing order, from 1 to n-k+1.
[B, i] = gsort(abs(s+z(1:n-k+1)), "g", "i");
if order == 1 then
// Set in i the indice which makes |s+z(i)| minimum
i = i(1);
else
// Set in i the indice which makes |s+z(i)| maximum
i = i($);
end
// Take into account for entry #i
s = s + z(i);
// Switch entries #n-k+1 and #i
t = z(i);
z(i) = z(n-k+1);
z(n-k+1) = t;
end
// Take into account for entry #i
s = s + z(i)
// Switch entries #n-k+1 and #i
t = z(i)
z(i) = z(n-k+1)
z(n-k+1) = t
end
endfunction
macros/accsum_compsum.sci
77
88
99
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
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
4345
44
45
46
47
48
49
50
51
52
53
54
55
56
46
47
48
49
50
51
52
53
54
55
56
57
58
5759
5860
// 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
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
[lhs, rhs] = argn()
apifun_checkrhs ( "accsum_compsum" , rhs , 1:1 )
apifun_checklhs ( "accsum_compsum" , lhs , 0:2 )
//
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
[lhs, rhs] = argn();
apifun_checkrhs("accsum_compsum", rhs, 1:1);
apifun_checklhs("accsum_compsum", lhs, 0:2);
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
77
88
99
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
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
4444
45
46
47
48
49
50
51
52
53
45
46
47
48
49
50
51
52
53
5454
5555
// 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 by a straightforward sum.
// This is called "recursive summation" in SNAA.
//
// 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
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 by a straightforward sum.
// This is called "recursive summation" in SNAA.
//
// 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
[lhs, rhs] = argn()
apifun_checkrhs ( "accsum_straight" , rhs , 1:1 )
apifun_checklhs ( "accsum_straight" , lhs , 0:1 )
//
n = size(x,"*")
s = 0
for i = 1 : n
s = s + x(i)
end
[lhs, rhs] = argn();
apifun_checkrhs("accsum_straight", rhs, 1:1);
apifun_checklhs("accsum_straight", lhs, 0:1);
//
n = size(x, "*");
s = 0;
for i = 1 : n
s = s + x(i);
end
endfunction
macros/accsum_scs.sci
77
88
99
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
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
4446
45
46
47
48
49
50
51
52
53
54
47
48
49
50
51
52
53
54
55
56
5557
5658
// 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
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
[lhs, rhs] = argn()
apifun_checkrhs ( "accsum_scs" , rhs , 1:1 )
apifun_checklhs ( "accsum_scs" , lhs , 0:2 )
//
s = 0
e = 0
n = size(x,"*")
for i = 1 : n
[s,e] = accsum_fasttwosum ( s , x(i) + e )
end
[lhs, rhs] = argn();
apifun_checkrhs("accsum_scs", rhs, 1:1);
apifun_checklhs("accsum_scs", lhs, 0:2);
s = 0;
e = 0;
n = size(x, "*");
for i = 1 : n
[s, e] = accsum_fasttwosum(s, x(i) + e);
end
endfunction
macros/accsum_getpath.sci
66
77
88
9
9
1010
11
11
1212
13
13
1414
1515
1616
1717
1818
19
19
2020
2121
2222
2323
24
25
26
27
28
29
24
25
26
27
28
29
3030
3131
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
function path = accsum_getpath ( )
function path = accsum_getpath()
// Returns the path to the current module.
//
//
// Calling Sequence
// path = accsum_getpath ( )
// path = accsum_getpath()
//
// Parameters
// path : a 1-by-1 matrix of strings, the path to the current module.
//
// Examples
// path = accsum_getpath ( )
// path = accsum_getpath()
//
// Authors
// 2010 - DIGITEO - Michael Baudin
[lhs, rhs] = argn()
apifun_checkrhs ( "accsum_getpath" , rhs , 0:0 )
apifun_checklhs ( "accsum_getpath" , lhs , 0:1 )
//
path = get_function_path("accsum_getpath")
path = fullpath(fullfile(fileparts(path),".."))
[lhs, rhs] = argn();
apifun_checkrhs("accsum_getpath", rhs, 0:0);
apifun_checklhs("accsum_getpath", lhs, 0:1);
path = get_function_path("accsum_getpath");
path = fullpath(fullfile(fileparts(path), ".."));
endfunction
macros/accsum_shuffle.sci
77
88
99
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
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
4141
42
43
44
45
46
47
42
43
44
45
46
47
4848
4949
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
function y = accsum_shuffle ( x )
// Randomly shuffles the input.
//
// Calling Sequence
// y = accsum_shuffle ( x )
//
// Parameters
// x : a m-by-n, matrix of doubles
// y : a m-by-n, matrix of doubles, the randomly permuted values.
//
// Description
// Randomly shuffles the input.
// This function makes use of the grand function and
// modifies its state.
//
// Examples
// // Use with larger r
// x = accsum_wilkinson(10);
// for i = 1 : 10
// x = accsum_shuffle ( x );
// [s1,e1] = accsum_compsum ( x );
// [s2,e2] = accsum_dblcompsum ( x );
// mprintf("#%5d CS=%e + %e DCS=%e + %e\n",i,s1,e1,s2,e2);
// end
//
// Authors
// Michael Baudin, 2010
//
// Bibliography
// "Stability and numerical accuracy of algorithms", Nicolas Higham
// "Handbook of Floating Point Computations", Muller et al
function y = accsum_shuffle(x)
// Randomly shuffles the input.
//
// Calling Sequence
// y = accsum_shuffle(x)
//
// Parameters
// x : a m-by-n, matrix of doubles
// y : a m-by-n, matrix of doubles, the randomly permuted values.
//
// Description
// Randomly shuffles the input.
// This function makes use of the grand function and
// modifies its state.
//
// Examples
// // Use with larger r
// x = accsum_wilkinson(10);
// for i = 1 : 10
// x = accsum_shuffle(x);
// [s1,e1] = accsum_compsum(x);
// [s2,e2] = accsum_dblcompsum(x);
// mprintf("#%5d CS=%e + %e DCS=%e + %e\n",i,s1,e1,s2,e2);
// end
//
// Authors
// Michael Baudin, 2010
//
// Bibliography
// "Stability and numerical accuracy of algorithms", Nicolas Higham
// "Handbook of Floating Point Computations", Muller et al
[lhs, rhs] = argn()
apifun_checkrhs ( "accsum_shuffle" , rhs , 1:1 )
apifun_checklhs ( "accsum_shuffle" , lhs , 0:1 )
//
x = x(:)
y = grand(1,"prm",x)
[lhs, rhs] = argn();
apifun_checkrhs("accsum_shuffle", rhs, 1:1);
apifun_checklhs("accsum_shuffle", lhs, 0:1);
x = x(:);
y = grand(1, "prm", x);
endfunction
macros/accsum_higham.sci
77
88
99
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
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
4747
48
49
50
51
52
53
48
49
50
51
52
53
54
55
56
5457
5558
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
function x = accsum_higham ( )
// Returns an example designed by Higham.
//
// Calling Sequence
// x = accsum_higham ( )
//
// Parameters
// x : a 4-by-1 matrix of doubles
//
// Description
// This example is so that the orderings produce
// different results.
//
// The exact sum is 1.
//
// Increasing magnitude and Psum (dynamic ordering, with order=1)
// produce 0, while decreasing magnitude produce 1.
//
// Examples
// x = accsum_higham ( )
// sum(x)
// // With direct input x: expected = 0
// s = accsum_straight(x)
// // With increasing magnitude ordering: expected = 0
// s = accsum_straight(accsum_order ( x , 4 ))
// // With dynamic magnitude ordering: expected = 0
// s = accsum_orderdynamic ( x , 1 )
// // With decreasing magnitude ordering: expected = 1
// s = accsum_straight(accsum_order ( x , 5 ))
//
// 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/img16.htm
function x = accsum_higham()
// Returns an example designed by Higham.
//
// Calling Sequence
// x = accsum_higham()
//
// Parameters
// x : a 4-by-1 matrix of doubles
//
// Description
// This example is so that the orderings produce
// different results.
//
// The exact sum is 1.
//
// Increasing magnitude and Psum (dynamic ordering, with order=1)
// produce 0, while decreasing magnitude produce 1.
//
// Examples
// x = accsum_higham()
// sum(x)
// // With direct input x: expected = 0
// s = accsum_straight(x)
// // With increasing magnitude ordering: expected = 0
// s = accsum_straight(accsum_order ( x , 4 ))
// // With dynamic magnitude ordering: expected = 0
// s = accsum_orderdynamic ( x , 1 )
// // With decreasing magnitude ordering: expected = 1
// s = accsum_straight(accsum_order ( x , 5 ))
//
// 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/img16.htm
[lhs, rhs] = argn()
apifun_checkrhs ( "accsum_higham" , rhs , 0:0 )
apifun_checklhs ( "accsum_higham" , lhs , 0:1 )
//
M = 2^53
x = [1;M;2*M;-3*M]
[lhs, rhs] = argn();
apifun_checkrhs("accsum_higham", rhs, 0:0);
apifun_checklhs("accsum_higham", lhs, 0:1);
M = 2 ^ 53;
x = [1
M
2 * M
-3 * M];
endfunction
macros/accsum_priestx.sci
77
88
99
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
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
3535
36
37
38
39
40
41
42
43
36
37
38
39
40
41
42
43
44
4445
45
46
47
48
49
50
51
52
46
47
48
49
50
51
52
53
5354
5455
// 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
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
//
// 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
[lhs, rhs] = argn()
apifun_checkrhs ( "accsum_priestx" , rhs , 0:0 )
apifun_checklhs ( "accsum_priestx" , lhs , 0:1 )
//
t = 53
x(1) = 2^(t+1)
x(2) = 2^(t+1) - 2
x(3:6) = -(2^t - 1)
[lhs, rhs] = argn();
apifun_checkrhs("accsum_priestx", rhs, 0:0);
apifun_checklhs("accsum_priestx", lhs, 0:1);
t = 53;
x(1) = 2 ^ (t+1);
x(2) = 2 ^ (t+1) - 2;
x(3:6) = - (2 ^ t - 1);
endfunction
macros/accsum_twosum.sci
77
88
99
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
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
4749
48
49
50
51
52
53
54
55
56
57
50
51
52
53
54
55
56
57
58
59
5860
5961
// 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
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
[lhs, rhs] = argn()
apifun_checkrhs ( "accsum_twosum" , rhs , 2:2 )
apifun_checklhs ( "accsum_twosum" , lhs , 0:2 )
//
s = a + b
ap = s - b
bp = s - ap
da = a - ap
db = b - bp
t = da + db
[lhs, rhs] = argn();
apifun_checkrhs("accsum_twosum", rhs, 2:2);
apifun_checklhs("accsum_twosum", lhs, 0:2);
s = a + b;
ap = s - b;
bp = s - ap;
da = a - ap;
db = b - bp;
t = da + db;
endfunction
macros/accsum_dblcompsum.sci
77
88
99
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
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
4244
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
6365
6466
// 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
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
[lhs, rhs] = argn()
apifun_checkrhs ( "accsum_dblcompsum" , rhs , 1:1 )
apifun_checklhs ( "accsum_dblcompsum" , lhs , 0:2 )
//
// 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
[lhs, rhs] = argn();
apifun_checkrhs("accsum_dblcompsum", rhs, 1:1);
apifun_checklhs("accsum_dblcompsum", lhs, 0:2);
// 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
77
88
99
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
10
11
4712
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
4843
49
50
51
44
45
46
47
48
49
50
5251
5352
// 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
[lhs, rhs] = argn()
apifun_checkrhs ( "accsum_sumcond" , rhs , 1:1 )
apifun_checklhs ( "accsum_sumcond" , lhs , 0:2 )
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)
[lhs, rhs] = argn();
apifun_checkrhs("accsum_sumcond", rhs, 1:1);
apifun_checklhs("accsum_sumcond", lhs, 0:2);
s = sum(x);
varf = sum(abs(x));
c = varf / abs(s);
endfunction
macros/accsum_dcs.sci
77
88
99
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
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
4850
49
50
51
52
53
54
55
56
57
58
59
60
51
52
53
54
55
56
57
58
59
60
61
62
6163
6264
// 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.
// The input data x must be ordered in decreasing magnitude.
// To do this, we may use the accsum_order
// function with order=5.
//
// 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)
// x = accsum_order ( x , 5 );
// [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
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.
// The input data x must be ordered in decreasing magnitude.
// To do this, we may use the accsum_order
// function with order=5.
//
// 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)
// x = accsum_order ( x , 5 );
// [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
[lhs, rhs] = argn()
apifun_checkrhs ( "accsum_dcs" , rhs , 1:1 )
apifun_checklhs ( "accsum_dcs" , lhs , 0:2 )
//
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
[lhs, rhs] = argn();
apifun_checkrhs("accsum_dcs", rhs, 1:1);
apifun_checklhs("accsum_dcs", lhs, 0:2);
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
77
88
99
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
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
5059
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
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
8695
8796
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
function y = accsum_order ( x , o )
// Re-order the matrix.
//
// Calling Sequence
// y = accsum_order ( x , o )
//
// Parameters
// x : a m-by-n, matrix of doubles
// o : a 1-by-1 matrix of doubles, a positive integer, the order. o = 1: no-op i.e. y=x, o = 2: increasing order, o = 3: decreasing order, o = 4: increasing magnitude order, o = 5: decreasing magnitude order (good order for DCS), o = 6: positive reverse from order 2, o = 7: negative reverse from order 2, o = 8: random order
// y : a m-by-n, matrix of doubles, the same content as x, but in different order.
//
// Description
// Re-order the matrix.
// If o=8, then this function uses the grand function and
// modifies its state.
//
// The increasing magnitude order (o=4) is good when summing nonnegative
// numbers by recursive summation (i.e. with accsum_straight),
// in the sense of having the smallest a priori forward error bound.
//
// The decreasing magnitude ordering order (o=5) is good when
// there is heavy cancellation, i.e. when the condition number
// returned by accsum_sumcond is large.
// In this case, the recursive summation (i.e. accsum_straight)
// can be used with decreasing magnitude ordering.
//
// Examples
// x = [-9.9 5.5 2.2 -3.3 -6.6 0 1.1]
// for o = 1 : 8
// y = accsum_order ( x , o );
// disp([o y'])
// end
//
// 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
function y = accsum_order(x, o)
// Re-order the matrix.
//
// Calling Sequence
// y = accsum_order(x, o)
//
// Parameters
// x : a m-by-n, matrix of doubles
// o : a 1-by-1 matrix of doubles, a positive integer, the order.
// o = 1 : no-op i.e. y = x,
// o = 2 : increasing order,
// o = 3 : decreasing order,
// o = 4 : increasing magnitude order,
// o = 5 : decreasing magnitude order (good order for DCS),
// o = 6 : positive reverse from order 2,
// o = 7 : negative reverse from order 2,
// o = 8 : random order
// y : a m-by-n, matrix of doubles, the same content as x,
// but in different order.
//
// Description
// Re-order the matrix.
// If o=8, then this function uses the grand function and
// modifies its state.
//
// The increasing magnitude order (o=4) is good when summing nonnegative
// numbers by recursive summation (i.e. with accsum_straight),
// in the sense of having the smallest a priori forward error bound.
//
// The decreasing magnitude ordering order (o=5) is good when
// there is heavy cancellation, i.e. when the condition number
// returned by accsum_sumcond is large.
// In this case, the recursive summation (i.e. accsum_straight)
// can be used with decreasing magnitude ordering.
//
// Examples
// x = [-9.9 5.5 2.2 -3.3 -6.6 0 1.1]
// for o = 1 : 8
// y = accsum_order(x, o);
// disp([o y'])
// end
//
// 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
[lhs, rhs] = argn()
apifun_checkrhs ( "accsum_order" , rhs , 2:2 )
apifun_checklhs ( "accsum_order" , lhs , 0:1 )
//
x = x(:)
n = size(x,"*")
select o
case 1
k = 1:n
case 2
[B,k] = gsort(x,"g","i")
case 3
[B,k] = gsort(x,"g","d")
case 4
[B,k] = gsort(abs(x),"g","i")
case 5
[B,k] = gsort(abs(x),"g","d")
case 6
k1 = find(x<=0)
[B1,k11] = gsort(x(k1),"g","i")
k2 = find(x>0)
[B,k22] = gsort(x(k2),"g","d")
k = [k1(k11) k2(k22)]
case 7
k1 = find(x<0)
[B,k11] = gsort(x(k1),"g","d")
k2 = find(x>=0)
[B,k22] = gsort(x(k2),"g","i")
k = [k1(k11) k2(k22)]
case 8
k = grand(1,"prm",(1:n)')
else
error ( msprintf("Unkown order %d",o) )
end
y = x(k)
[lhs, rhs] = argn();
apifun_checkrhs("accsum_order", rhs, 2:2);
apifun_checklhs("accsum_order", lhs, 0:1);
x = x(:);
n = size(x, "*");
select o
case 1
k = 1:n;
case 2
[B, k] = gsort(x, "g", "i");
case 3
[B, k] = gsort(x, "g", "d");
case 4
[B, k] = gsort(abs(x), "g", "i");
case 5
[B, k] = gsort(abs(x), "g", "d");
case 6
k1 = find(x<=0);
[B1, k11] = gsort(x(k1), "g", "i");
k2 = find(x>0);
[B, k22] = gsort(x(k2), "g", "d");
k = [k1(k11) k2(k22)];
case 7
k1 = find(x<0);
[B, k11] = gsort(x(k1), "g", "d");
k2 = find(x>=0);
[B, k22] = gsort(x(k2), "g", "i");
k = [k1(k11) k2(k22)];
case 8
k = grand(1, "prm", (1:n)');
else
error(msprintf("Unkown order %d", o));
end
y = x(k);
endfunction
src/c/accsum.h
1010
1111
1212
13
14
15
16
17
13
14
15
16
17
1818
19
19
2020
2121
2222
......
3232
3333
3434
35
35
36
3637
37
38
39
3840
39
41
42
4043
41
44
4245
4346
4447
45
4648
4749
#define _ACCSUM_H_
#ifdef _MSC_VER
#if LIBACCSUM_EXPORTS
#define ACCSUM_IMPORTEXPORT __declspec (dllexport)
#else
#define ACCSUM_IMPORTEXPORT __declspec (dllimport)
#endif
#if LIBACCSUM_EXPORTS
#define ACCSUM_IMPORTEXPORT __declspec (dllexport)
#else
#define ACCSUM_IMPORTEXPORT __declspec (dllimport)
#endif
#else
#define ACCSUM_IMPORTEXPORT
#define ACCSUM_IMPORTEXPORT
#endif
#undef __BEGIN_DECLS
__BEGIN_DECLS
// Sets s as the sum of x entries, with a double compensated summation.
void accsum_fdcs(double *x, int n, double *s, double *e);
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);
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, double *e);
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 , double *e);
void accsum_fcompsum(double * x, int n, double * s, double * e);
__END_DECLS
#endif /* _ACCSUM_H_ */
src/c/builder_c.sce
1313
1414
1515
16
17
16
17
1818
1919
2020
21
22
21
22
2323
2424
2525
files = "accsum.c";
ldflags = "";
if getos() == "Windows" then
// Configure the floating point unit so as to use
// SSE2 units, potentially instead of x87 registers.
// Configure the floating point unit so as to use
// SSE2 units, potentially instead of x87 registers.
cflags = "-DWIN32 -DACCSUM_EXPORTS /arch:SSE2";
libs = [];
else
// Configure the floating point unit so as to use
// SSE units, instead of registers.
// Configure the floating point unit so as to use
// SSE units, instead of registers.
// The x86 registers may use extended precision.
// http://gcc.gnu.org/wiki/FloatingPointMath
include1 = src_c_path;
src/c/accsum.c
2020
2121
2222
23
24
25
26
23
24
25
26
27
28
2729
28
2930
3031
3132
......
3637
3738
3839
39
40
4041
41
42
43
44
45
42
43
44
45
46
4647
47
48
49
50
51
52
53
54
48
49
50
51
52
53
54
55
5556
5657
5758
......
6465
6566
6667
67
68
6869
69
70
71
70
71
72
7273
73
74
75
76
77
78
79
80
74
75
76
77
78
79
80
81
8182
8283
8384
8485
8586
8687
87
88
8889
8990
9091
9192
92
93
9394
94
95
95
96
9697
97
98
99
100
101
102
103
104
105
106
98
99
100
101
102
103
104
105
106
107
107108
108109
109
110
110111
111
112
113
114
115
116
117
118
119
120
121
122
123
112
113
114
115
116
117
118
119
120
121
122
123
124
124125
125
126
127
128
129
130
131
132
133
134
135
126
127
128
129
130
131
132
133
134
135
136
136137
137138
138
139
139140
140
141
142
143
144
145
146
147
148
149
150
151
141
142
143
144
145
146
147
148
149
150
151
152
152153
153
154
155
154
155
156
156157
157
158
159
160
161
162
163
164
165
158
159
160
161
162
163
164
165
166
166167
167168
// https://hpcrd.lbl.gov/SCG/ocean/NRS/ECMWF/img14.htm
// https://hpcrd.lbl.gov/SCG/ocean/NRS/SCSsum.F
//
// "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.
// "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.
// 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
// "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, double *e)
void accsum_fdcs(double * x, int n, double * s, double * e)
{
double s1;
double e1;
double s2;
double e2;
int i;
double s1;
double e1;
double s2;
double e2;
int i;
*s = 0;
*e = 0;
for (i = 0; i< n; i++)
{
accsum_fasttwosum (*e, x[i], &s1, &e1 );
accsum_fasttwosum (*s, s1, &s2, &e2);
accsum_fasttwosum (s2, e1 + e2, s, e);
}
*s = 0;
*e = 0;
for (i = 0; i < n; i++)
{
accsum_fasttwosum(*e, x[i], &s1, &e1);
accsum_fasttwosum(*s, s1, &s2, &e2);
accsum_fasttwosum(s2, e1 + e2, s, e);
}
}
// Returns s and e such that s+e is the sum of x.
// "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, double *e)
void accsum_fscs(double * x, int n, double * s, double * e)
{
double s1;
double e1;
int i;
double s1;
double e1;
int i;
*s = 0;
*e = 0;
for (i = 0; i< n; i++)
{
accsum_fasttwosum ( *s , x[i] + *e , &s1, &e1);
*s = s1;
*e = e1;
}
*s = 0;
*e = 0;
for (i = 0; i < n; i++)
{
accsum_fasttwosum(*s, x[i] + *e, &s1, &e1);
*s = s1;
*e = e1;
}
}
// Returns s and t so that a + b = s + t exactly.
// where s contains the highest significant digits
// and t the lowest significant digits.
//
// Assumes that |a|>=|b|. If not, the arguments are switched.
// Assumes that |a|>=|b|. If not, the arguments are switched.
//
// Algorithm 4.3 in HFCC
// Due to Dekker (and Kahan)
// Assumes that we use round-to-nearest.
void accsum_fasttwosum(double a, double b, double *s, double *t)
void accsum_fasttwosum(double a, double b, double * s, double * t)
{
double tmp;
double z;
double tmp;
double z;
if ( fabs(a) < fabs(b) )
{
// Switch a and b
tmp = a;
a = b;
b = tmp;
}
*s = a + b;
z = *s - a;
*t = b - z;
if (fabs(a) < fabs(b))
{
// Switch a and b
tmp = a;
a = b;
b = tmp;
}
*s = a + b;
z = *s - a;
*t = b - z;
}
void accsum_twosum ( double a , double b , double * s , double * t )
void accsum_twosum(double a, double b, double * s, double * t)
{
// The twosum sum of a and b.
// 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 we use round-to-nearest.
// Bibliography
// "Stability and numerical accuracy of algorithms", Nicolas Higham
// "Handbook of Floating Point Computations", Muller et al
// The twosum sum of a and b.
// 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 we use round-to-nearest.
// Bibliography
// "Stability and numerical accuracy of algorithms", Nicolas Higham
// "Handbook of Floating Point Computations", Muller et al
double ap;
double bp;
double da;
double db;
*s = a + b;
ap = *s - b;
bp = *s - ap;
da = a - ap;
db = b - bp;
*t = da + db;
double ap;
double bp;
double da;
double db;
*s = a + b;
ap = *s - b;
bp = *s - ap;
da = a - ap;
db = b - bp;
*t = da + db;
}
void accsum_fcompsum ( double * x , int n , double * s , double * e )
void accsum_fcompsum(double * x, int n, double * s, double * e)
{
// The compensated sum of a matrix.
// Description
// 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
// Assumes that Scilab uses round-to-nearest.
// Bibliography
// "Stability and numerical accuracy of algorithms", Nicolas Higham
// "Handbook of Floating Point Computations", Muller et al
// The compensated sum of a matrix.
// Description
// 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
// Assumes that Scilab uses round-to-nearest.
// Bibliography
// "Stability and numerical accuracy of algorithms", Nicolas Higham
// "Handbook of Floating Point Computations", Muller et al
int i;
double t;
double y;
int i;
double t;
double y;
*s = 0;
*e = 0;
for (i = 0; i< n; i++)
{
t = *s;
y = x[i] + *e;
*s = t + y;
*e = (t - *s) + y;
}
*s = 0;
*e = 0;
for (i = 0; i < n; i++)
{
t = *s;
y = x[i] + *e;
*s = t + y;
*e = (t - *s) + y;
}
}

Archive Download the corresponding diff file

Revision: 33