# accsum Commit Details

Date: 2011-03-23 21:59:53 (7 years 9 months ago) Michael Baudin 18 17 Used consistently apifun. Created missing .ref Created higham example.

## File differences

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
 ␊ ␊ ␊ ␊ ␊ ␊ ␊ accsum_orderdynamicTry to find an order which makes the sum perform bad/good.␊ ␊ ␊ ␊ Calling Sequence␊ s = accsum_orderdynamic ( x , order )␊ ␊ ␊ ␊ Parametersx : a m-by-n, matrix of doublesorder : 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 accuracys : a 1-by-1 matrix of doubles, the sum␊ ␊ Description␊ Try to find an order which makes the sum perform bad/good.␊ ␊ ␊ Examples␊ ␊ AuthorsMichael Baudin, 2010␊ ␊ Bibliography"Stability and numerical accuracy of algorithms", Nicolas Higham"Handbook of Floating Point Computations", Muller et alhttps://hpcrd.lbl.gov/SCG/ocean/NRS/ECMWF/img16.htm␍␊ ␍␊ ␍␊ ␍␊ ␍␊ ␍␊ ␍␊ ␍␊ accsum_orderdynamicCompute 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␍␊ ␍␊ ␍␊ ␍␊ ␍␊ 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␍␊ ␍␊ ␍␊
help/en_US/support/accsum_order.xml
 35 35 36 36 37 37 38 38 39 39 40 40 41 41 ... ... 49 49 50 50 51 51 52 53 54 55 56 57 58 59 60 61 62 63 52 64 53 65 54 66
 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␍␊ 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.␍␊ ␍␊ 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.␍␊ ␍␊ ␍␊ ␍␊ ␍␊ ␍␊
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
 ␍␊ ␍␊ ␍␊ ␍␊ ␍␊ ␍␊ ␍␊ ␍␊ accsum_highamReturns 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␍␊ ␍␊ ␍␊ ␍␊ ␍␊ 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␍␊ ␍␊ ␍␊
help/en_US/support/accsum_priestx.xml
 59 59 60 60 61 61 62 63 62 63 64 64 65 65 66 66
 Examples␍␊ ␍␊ ␍␊ ␍␊
help/en_US/support/accsum_sumcond.xml
 24 24 25 25 26 26 27 27 28 28 29 29 30 30 31 ... ... 60 61 61 62 62 63 63 64 64 65 65 66 66 67
 ␍␊ Calling Sequence␍␊ ␍␊ x = accsum_sumcond ( )␍␊ s = accsum_sumcond ( x )␍␊ [s,c] = accsum_sumcond ( x )␍␊ ␍␊ ␍␊ ␍␊ [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/accsum_straight.xml
 42 42 43 43 44 44 45 45 46 46 47 47 48 48 49
 ␍␊ Description␍␊ ␍␊ Computes the sum.␍␊ Computes the sum by a straightforward sum.␍␊ This is called "recursive summation" in SNAA.␍␊ ␍␊ ␍␊ Uses an unvectorized loop.␍␊
help/en_US/update_help.sce
 16 16 17 17 18 18 19 20 19 21 20 22 21 23 22 24 23 25 24 26 25 27 26 28 27 29 28 30 29 31 30 32 31 32 33 33 34 34 35 35 ... ... 37 37 38 38 39 39 40 40 41 41 42 42 43 43 ... ... 52 52 53 53 54 54 55 55 56 56 57 57 58 58
 "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 );␍␊
etc/accsum.start
 32 32 33 33 34 34 35 35 36 36 37 37 38 38
 // 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));␊
 13 13 14 14 15 15 16 17 16 17 18 18 19 19 20 20 21 22 23 24 25 26 27 28 29 21 30 22 31 23 24 32 25 33 26 27 34 28 29 35 36 30 37 38 39 40 31 41 32 42 33 43 34 44 35 45 46 36 47 37 48 38 49 39 50 40 51 52 41 53 42 54 43 55
 --------␊ ␊ 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␊ ------␊
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
 9 9 10 10 11 11 12 12 13 13 14 14 15
 // <-- 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_orderdynamic.sci
 19 19 20 20 21 21 22 22 23 24 25 26 27 28 29 30 31 32 33 34 23 35 24 36 25 37 ... ... 47 59 48 60 49 61 62 63 64 65 50 66 51 67 52 68
 // 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
 41 41 42 42 43 43 44 45 46 47 44 48 45 49 46 50
 // "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
 18 18 19 19 20 20 21 21 22 22 23 23 24 24 25 ... ... 41 42 42 43 43 44 45 46 47 48 44 49 45 50 46 51
 // 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
 42 42 43 43 44 44 45 46 47 48 45 49 46 50 47 51
 // 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
 21 21 22 22 23 23 24 25 26 27 24 28 25 29 26 30
 // 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
 39 39 40 40 41 41 42 43 44 45 42 46 43 47 44 48
 // "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_twosum.sci
 45 45 46 46 47 47 48 49 50 51 48 52 49 53 50 54
 // "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_priestx.sci
 30 30 31 31 32 32 33 34 33 34 35 35 36 36 37 37 ... ... 42 42 43 43 44 44 45 46 47 48 45 49 46 50 47 51
 // ␊ // 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_dblcompsum.sci
 40 40 41 41 42 42 43 44 45 46 43 47 44 48 45 49
 // "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
 30 30 31 31 32 32 33 33 34 34 35 35 36 36 ... ... 41 41 42 42 43 43 44 45 46 47 48 44 49 45 50 46 51
 // [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
 42 42 43 43 44 44 45 46 47 48 45 49 46 50 47 51
 // 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,"*")␊
macros/accsum_order.sci
 15 15 16 16 17 17 18 18 19 19 20 20 21 21 22 22 23 23 24 24 25 26 27 28 29 30 31 32 33 34 25 35 26 36 27 37 ... ... 38 48 39 49 40 50 51 52 53 54 41 55 42 56 43 57
 // ␊ // 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
 44 44 45 45 46 46 47 48 49 50 47 51 48 52 49 53
 // "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
 40 40 41 41 42 42 43 44 45 46 43 47 44 48 45 49
 // "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␊

Revision: 18