accsum

accsum Commit Details

Date:2011-03-23 21:59:53 (7 years 3 months ago)
Author:Michael Baudin
Commit:18
Parents: 17
Message:Used consistently apifun. Created missing .ref Created higham example.
Changes:
A/help/en_US/support/accsum_higham.xml
A/tests/unit_tests/shuffle.dia.ref
A/tests/unit_tests/higham.dia.ref
A/tests/unit_tests/higham.tst
A/tests/unit_tests/wilkinson.dia.ref
M/macros/accsum_dblcompsum.sci
M/tests/unit_tests/scs.dia.ref
M/help/en_US/update_help.sce
M/macros/accsum_sumcond.sci
M/macros/accsum_dcs.sci
M/help/en_US/support/accsum_priestx.xml
M/macros/accsum_order.sci
M/macros/accsum_fasttwosum.sci
M/tests/unit_tests/shuffle.tst
M/macros/accsum_wilkinson.sci
M/readme.txt
M/help/en_US/accsum_orderdynamic.xml
M/macros/accsum_orderdynamic.sci
M/help/en_US/support/accsum_sumcond.xml
M/macros/accsum_compsum.sci
M/help/en_US/accsum_straight.xml
M/macros/accsum_straight.sci
M/etc/accsum.start
M/help/en_US/support/accsum_order.xml
M/tests/unit_tests/dcs.dia.ref
M/macros/accsum_scs.sci
M/macros/accsum_shuffle.sci
M/macros/accsum_getpath.sci
M/macros/accsum_priestx.sci
M/macros/accsum_twosum.sci

File differences

tests/unit_tests/shuffle.dia.ref
1
2
3
4
5
6
7
8
9
10
11
12
13
14
// Copyright (C) 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
//
// <-- JVM NOT MANDATORY -->
// <-- ENGLISH IMPOSED -->
x = (0:10);
y1 = accsum_shuffle ( x );
y2 = accsum_shuffle ( x );
assert_checktrue(~(and(y1==y2)));
tests/unit_tests/higham.dia.ref
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
// Copyright (C) 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
//
// <-- JVM NOT MANDATORY -->
// <-- ENGLISH IMPOSED -->
x = accsum_higham ( );
// With direct input x: expected = 0
s = accsum_straight(x);
assert_checkequal(s,0);
// With increasing magnitude ordering: expected = 0
s = accsum_straight(accsum_order ( x , 4 ));
// With dynamic magnitude ordering: expected = 0
assert_checkequal(s,0);
s = accsum_orderdynamic ( x , 1 );
assert_checkequal(s,0);
// With decreasing magnitude ordering: expected = 1
s = accsum_straight(accsum_order ( x , 5 ));
assert_checkequal(s,1);
tests/unit_tests/shuffle.tst
99
1010
1111
12
1213
1314
1415
// <-- JVM NOT MANDATORY -->
// <-- ENGLISH IMPOSED -->
x = (0:10);
y1 = accsum_shuffle ( x );
y2 = accsum_shuffle ( x );
assert_checktrue(~(and(y1==y2)));
tests/unit_tests/higham.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
// Copyright (C) 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
//
// <-- JVM NOT MANDATORY -->
// <-- ENGLISH IMPOSED -->
x = accsum_higham ( );
// With direct input x: expected = 0
s = accsum_straight(x);
assert_checkequal(s,0);
// With increasing magnitude ordering: expected = 0
s = accsum_straight(accsum_order ( x , 4 ));
// With dynamic magnitude ordering: expected = 0
assert_checkequal(s,0);
s = accsum_orderdynamic ( x , 1 );
assert_checkequal(s,0);
// With decreasing magnitude ordering: expected = 1
s = accsum_straight(accsum_order ( x , 5 ));
assert_checkequal(s,1);
tests/unit_tests/dcs.dia.ref
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
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
// Copyright (C) 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
//
// <-- JVM NOT MANDATORY -->
// <-- ENGLISH IMPOSED -->
[s,e] = accsum_DCS ( [2 1] );
assert_checkequal(s,3);
assert_checkequal(e,0);
//
[s,e] = accsum_DCS ( [1 2] );
assert_checkequal(s,3);
assert_checkequal(e,0);
//
[s,e] = accsum_DCS ( [2 1] );
assert_checkequal(s,3);
assert_checkequal(e,0);
//
[s,e] = accsum_DCS ( [1 2] );
assert_checkequal(s,3);
assert_checkequal(e,0);
//
x = accsum_wilkinson(10);
[s,e] = accsum_DCS ( x );
assert_checkalmostequal(s,1024,1.e-12);
assert_checkalmostequal(e,-3.786e-14,[],1.e-10);
// Copyright (C) 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
//
// <-- JVM NOT MANDATORY -->
// <-- ENGLISH IMPOSED -->
[s,e] = accsum_dcs ( [2 1] );
assert_checkequal(s,3);
assert_checkequal(e,0);
//
[s,e] = accsum_dcs ( [1 2] );
assert_checkequal(s,3);
assert_checkequal(e,0);
//
[s,e] = accsum_dcs ( [2 1] );
assert_checkequal(s,3);
assert_checkequal(e,0);
//
[s,e] = accsum_dcs ( [1 2] );
assert_checkequal(s,3);
assert_checkequal(e,0);
//
x = accsum_wilkinson(10);
[s,e] = accsum_dcs ( x );
assert_checkalmostequal(s,1024,1.e-12);
assert_checkalmostequal(e,-3.786e-14,[],1.e-10);
tests/unit_tests/wilkinson.dia.ref
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
// Copyright (C) 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
//
// <-- JVM NOT MANDATORY -->
// <-- ENGLISH IMPOSED -->
// 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]';
assert_checkalmostequal(x,e,2*%eps);
tests/unit_tests/scs.dia.ref
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
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
// Copyright (C) 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
//
// <-- JVM NOT MANDATORY -->
// <-- ENGLISH IMPOSED -->
[s,e] = accsum_SCS ( [2 1] );
assert_checkequal(s,3);
assert_checkequal(e,0);
//
[s,e] = accsum_SCS ( [1 2] );
assert_checkequal(s,3);
assert_checkequal(e,0);
//
[s,e] = accsum_SCS ( [2 1] );
assert_checkequal(s,3);
assert_checkequal(e,0);
//
[s,e] = accsum_SCS ( [1 2] );
assert_checkequal(s,3);
assert_checkequal(e,0);
//
x = accsum_wilkinson(10);
[s,e] = accsum_SCS ( x );
assert_checkalmostequal(s,1024,1.e-12);
assert_checkalmostequal(e,-3.786e-14,[],1.e-10);
// Copyright (C) 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
//
// <-- JVM NOT MANDATORY -->
// <-- ENGLISH IMPOSED -->
[s,e] = accsum_scs ( [2 1] );
assert_checkequal(s,3);
assert_checkequal(e,0);
//
[s,e] = accsum_scs ( [1 2] );
assert_checkequal(s,3);
assert_checkequal(e,0);
//
[s,e] = accsum_scs ( [2 1] );
assert_checkequal(s,3);
assert_checkequal(e,0);
//
[s,e] = accsum_scs ( [1 2] );
assert_checkequal(s,3);
assert_checkequal(e,0);
//
x = accsum_wilkinson(10);
[s,e] = accsum_scs ( x );
assert_checkalmostequal(s,1024,1.e-12);
assert_checkalmostequal(e,-3.786e-14,[],1.e-10);
macros/accsum_order.sci
1515
1616
1717
18
18
1919
2020
2121
2222
2323
2424
25
26
27
28
29
30
31
32
33
34
2535
2636
2737
......
3848
3949
4050
51
52
53
54
4155
4256
4357
//
// 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
// 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]
// "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
macros/accsum_fasttwosum.sci
4444
4545
4646
47
48
49
50
4751
4852
4953
// "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
macros/accsum_wilkinson.sci
4040
4141
4242
43
44
45
46
4347
4448
4549
// "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
macros/accsum_orderdynamic.sci
1919
2020
2121
22
22
23
24
25
26
27
28
29
30
31
32
33
34
2335
2436
2537
......
4759
4860
4961
62
63
64
65
5066
5167
5268
// s : a 1-by-1 matrix of doubles, the sum
//
// Description
// Try to find an order which makes the sum perform bad/good.
// 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
// "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,"*")
macros/accsum_compsum.sci
4141
4242
4343
44
45
46
47
4448
4549
4650
// "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,"*")
macros/accsum_straight.sci
1818
1919
2020
21
21
22
2223
2324
2425
......
4142
4243
4344
45
46
47
48
4449
4550
4651
// s : a 1-by-1, matrix of doubles, the sum
//
// Description
// Computes the sum.
// 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
// "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
macros/accsum_scs.sci
4242
4343
4444
45
46
47
48
4549
4650
4751
// 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,"*")
macros/accsum_getpath.sci
2121
2222
2323
24
25
26
27
2428
2529
2630
// 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),".."))
endfunction
macros/accsum_shuffle.sci
3939
4040
4141
42
43
44
45
4246
4347
4448
// "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)
endfunction
macros/accsum_priestx.sci
3030
3131
3232
33
34
33
34
3535
3636
3737
......
4242
4343
4444
45
46
47
48
4549
4650
4751
//
// Examples
// x = accsum_priestx ( )
// [s,e] = accsum_SCS ( x ) // = 3 + 0
// [s,e] = accsum_DCS ( x ) // = 2 + 0
// [s,e] = accsum_scs ( x ) // = 3 + 0
// [s,e] = accsum_dcs ( x ) // = 2 + 0
//
// Authors
// "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
macros/accsum_twosum.sci
4545
4646
4747
48
49
50
51
4852
4953
5054
// "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
macros/accsum_dblcompsum.sci
4040
4141
4242
43
44
45
46
4347
4448
4549
// "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)
macros/accsum_sumcond.sci
3030
3131
3232
33
33
3434
3535
3636
......
4141
4242
4343
44
45
46
47
48
4449
4550
4651
// [s,c] = accsum_sumcond(x)
// for o = 1 : 8
// xo = accsum_order ( x , o );
// [s,e] = accsum_DCS ( xo );
// [s,e] = accsum_dcs ( xo );
// disp([o s e])
// end
//
// "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 )
//
s = sum(x)
varf = sum(abs(x))
c = varf / abs(s)
macros/accsum_dcs.sci
4242
4343
4444
45
46
47
48
4549
4650
4751
// 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,"*")
help/en_US/update_help.sce
1616
1717
1818
19
2019
2120
2221
2322
24
23
2524
2625
2726
2827
2928
3029
31
30
3231
32
3333
3434
3535
......
3737
3838
3939
40
40
4141
4242
4343
......
5252
5353
5454
55
55
5656
5757
5858
"accsum_dcs"
"accsum_fasttwosum"
"accsum_orderdynamic"
"accsum_priestx"
"accsum_scs"
"accsum_straight"
"accsum_twosum"
];
macrosdir = helpdir +"../../macros";
macrosdir = cwd +"../../macros";
demosdir = [];
modulename = "accsum";
helptbx_helpupdate ( funmat , helpdir , macrosdir , demosdir , modulename , %t );
//
// Generate the library help
mprintf("Updating accsum/support\n");
helpdir = cwd;
helpdir = fullfile(cwd,"support");
funmat = [
"accsum_higham"
"accsum_order"
"accsum_priestx"
"accsum_shuffle"
"accsum_wilkinson"
"accsum_getpath"
];
macrosdir = helpdir +"../../macros";
macrosdir = cwd +"../../macros";
demosdir = [];
modulename = "accsum";
helptbx_helpupdate ( funmat , helpdir , macrosdir , demosdir , modulename , %t );
"accsum_fscs"
"accsum_fcompsum"
];
macrosdir = fullfile(helpdir,"pseudomacros");
macrosdir = fullfile(cwd,"pseudomacros");
demosdir = [];
modulename = "accsum";
helptbx_helpupdate ( funmat , helpdir , macrosdir , demosdir , modulename , %t );
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
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
<?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>
<?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>Compute the sum with a dynamic re-ordering.</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.
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.
</para>
<para>
This algorithm is called "Psum" in the Chapter 4 "Summation" of
SNAA.
</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-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/img16.htm</para>
</refsection>
</refentry>
help/en_US/support/accsum_sumcond.xml
2424
2525
2626
27
27
28
2829
2930
3031
......
6061
6162
6263
63
64
6465
6566
6667
<refsynopsisdiv>
<title>Calling Sequence</title>
<synopsis>
x = accsum_sumcond ( )
s = accsum_sumcond ( x )
[s,c] = accsum_sumcond ( x )
</synopsis>
</refsynopsisdiv>
[s,c] = accsum_sumcond(x)
for o = 1 : 8
xo = accsum_order ( x , o );
[s,e] = accsum_DCS ( xo );
[s,e] = accsum_dcs ( xo );
disp([o s e])
end
help/en_US/support/accsum_order.xml
3535
3636
3737
38
38
3939
4040
4141
......
4949
5050
5151
52
53
54
55
56
57
58
59
60
61
62
63
5264
5365
5466
<varlistentry><term>x :</term>
<listitem><para> a m-by-n, matrix of doubles</para></listitem></varlistentry>
<varlistentry><term>o :</term>
<listitem><para> 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</para></listitem></varlistentry>
<listitem><para> 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</para></listitem></varlistentry>
<varlistentry><term>y :</term>
<listitem><para> a m-by-n, matrix of doubles, the same content as x, but in different order.</para></listitem></varlistentry>
</variablelist>
modifies its state.
</para>
<para>
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.
</para>
<para>
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.
</para>
<para>
</para>
</refsection>
help/en_US/support/accsum_higham.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
<?xml version="1.0" encoding="UTF-8"?>
<!--
*
* This help file was generated from accsum_higham.sci using help_from_sci().
*
-->
<refentry version="5.0-subset Scilab" xml:id="accsum_higham" 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_higham</refname><refpurpose>Returns an example designed by Higham.</refpurpose>
</refnamediv>
<refsynopsisdiv>
<title>Calling Sequence</title>
<synopsis>
x = accsum_higham ( )
</synopsis>
</refsynopsisdiv>
<refsection>
<title>Parameters</title>
<variablelist>
<varlistentry><term>x :</term>
<listitem><para> a 4-by-1 matrix of doubles</para></listitem></varlistentry>
</variablelist>
</refsection>
<refsection>
<title>Description</title>
<para>
This example is so that the orderings produce
different results.
</para>
<para>
The exact sum is 1.
</para>
<para>
Increasing magnitude and Psum (dynamic ordering, with order=1)
produce 0, while decreasing magnitude produce 1.
</para>
<para>
</para>
</refsection>
<refsection>
<title>Examples</title>
<programlisting role="example"><![CDATA[
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 ))
]]></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/img16.htm</para>
</refsection>
</refentry>
help/en_US/support/accsum_priestx.xml
5959
6060
6161
62
63
62
63
6464
6565
6666
<title>Examples</title>
<programlisting role="example"><![CDATA[
x = accsum_priestx ( )
[s,e] = accsum_SCS ( x ) // = 3 + 0
[s,e] = accsum_DCS ( x ) // = 2 + 0
[s,e] = accsum_scs ( x ) // = 3 + 0
[s,e] = accsum_dcs ( x ) // = 2 + 0
]]></programlisting>
</refsection>
</refentry>
help/en_US/accsum_straight.xml
4242
4343
4444
45
45
46
4647
4748
4849
<refsection>
<title>Description</title>
<para>
Computes the sum.
Computes the sum by a straightforward sum.
This is called "recursive summation" in SNAA.
</para>
<para>
Uses an unvectorized loop.
etc/accsum.start
3232
3333
3434
35
35
3636
3737
3838
// load gateways
// =============================================================================
if ( %t ) then
if ( %f ) then
mprintf("\tLoad gateways\n");
ilib_verbose(0);
exec( pathconvert(root_tlbx+"/sci_gateway/loader_gateway.sce",%f));
readme.txt
1313
1414
1515
16
17
16
17
1818
1919
2020
21
22
23
24
25
26
27
28
29
2130
2231
23
2432
2533
26
2734
28
2935
36
3037
38
39
40
3141
3242
3343
3444
3545
46
3647
3748
3849
3950
40
51
52
4153
4254
4355
--------
The following is a list of the current accsum functions :
* accsum_DCS : A Doubly Self Compensated Sum algorithm
* accsum_SCS : A Self Compensated Sum algorithm
* accsum_dcs : A Doubly Self Compensated Sum algorithm
* accsum_scs : A Self Compensated Sum algorithm
* accsum_compsum : The compensated sum of a matrix.
* accsum_dblcompsum : The doubly compensated sum of a matrix.
* accsum_fasttwosum : The fast2sum sum of a and b.
* accsum_orderdynamic : Returns the sum with a dynamic re-ordering.
* accsum_straight : The straightforward sum of a matrix.
* accsum_twosum : The twosum sum of a and b.
* accsum_fcompsum : The compensated sum of a matrix.
* accsum_fdcs : A Doubly Self Compensated Sum algorithm
* accsum_fscs : A Self Compensated Sum algorithm
and support functions:
* accsum_getpath : Returns the path to the current module.
* accsum_order : Re-order the matrix.
* accsum_orderdynamic : Try to find an order which makes the sum perform bad/good.
* accsum_priestx : A difficult example for SCS by Priest.
* accsum_shuffle : Randomly shuffles the input.
* accsum_straight : The straightforward sum of a matrix.
* accsum_sumcond : Condition number of the sum function.
* accsum_twosum : The twosum sum of a and b.
* accsum_wilkinson : A test vector by Wilkinson.
* accsum_higham : Returns an example designed by Higham.
The accsum_fcompsum, accsum_fdcs and accsum_fscs functions
are based on compiled source code and are faster than the other.
Dependencies
----
* This modules depends on the assert module.
* This modules depends on the helptbx module.
* This modules depends on the apifun module.
TODO
----
* Fix the bug in accsum_fdcs.
* Implement pairwise summation.
* Implement optimal summation.
Author
------

Archive Download the corresponding diff file

Revision: 18