Add code for reproducing experiments in the survey paper

.gitmodules

@@ -0,0 +1,6 @@
+[submodule "deps/chop"]
+	path = deps/chop
+	url = git@github.com:higham/chop.git
+[submodule "deps/cpfloat"]
+	path = deps/cpfloat
+	url = git@github.com:mfasi/cpfloat.git

LICENSE

@@ -0,0 +1,26 @@
+BSD 2-Clause License
+
+Copyright (c) 2022, Matteo Croci, Massimiliano Fasi, Nicholas J. Higham,
+Theo Mary, and Mantas Mikaitis
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+1. Redistributions of source code must retain the above copyright notice, this
+   list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright notice,
+   this list of conditions and the following disclaimer in the documentation
+   and/or other materials provided with the distribution.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

ODE_tests.m

@@ -0,0 +1,163 @@
+% ODE_tests.m   Test accuracy of ODE solvers in different rounding modes.
+%  This code uses the function chop at https://github.com/higham/chop.
+%
+% References:
+%   [1] N. J. Higham, S. Pranesh. Simulating Low Precision Floating-Point
+%       Arithmetic. SIAM J. Sci. Comput., 41(5), pp. 585–602. October 2019.
+%       http://dx.doi.org/10.1137/19M1251308
+%
+%   [2] M. Fasi, M. Mikaitis. Algorithms for stochastically rounded
+%       elementary arithmetic operations in IEEE 754 floating-point
+%       arithmetic. IEEE Trans. Emerg. Topics Comput., 9(3): 1451–1466.
+%       July 2021.
+%       http://dx.doi.org/10.1109/TETC.2021.3069165
+%
+%   [3] M. Croci, M. Fasi, N. J. Higham, T. Mary, M. Mikaitis.
+%       Stochastic Rounding: Implementation, Error Analysis, and
+%       Applications. Tech. Report 2021.17, Manchester Institute for
+%       Mathematical Sciences, The University of Manchester, UK.
+%       October 2022. Revised January 2022.
+
+% Clear chop options and reset the PRNG seed.
+clear options
+rng(1)
+
+% Set the number of times to repeat the SR experiments.
+rep = 10;
+
+% Choose a testcase (only 0 and 1 are supported)
+if ~exist('testcase', 'var')
+    testcase = 1;
+end
+
+% Set up initial ODE conditions.
+a = 0;
+if (testcase == 0)
+    b = 1.0;
+    y0 = 0.015625;
+elseif (testcase == 1)
+    b = 0.015625;
+    y0 = 1.0;
+else
+    error('This value of testcase is not supported.');
+end
+
+% Exact solution to the exponential decay ODE.
+if (testcase == 0)
+    yexact = exp(-b)*y0;
+else
+    yexact = exp(-b/20)*y0;
+end
+
+% Decay function.
+if (testcase == 0)
+    decay_ODE_acc = @(y, options)(-y);
+    decay_ODE_inacc = @(y, options)(-chop(y, options));
+else
+    decay_ODE_acc = @(y, options)(-y/20);
+    decay_ODE_inacc = @(y, options)(-chop(y/20, options));
+end
+
+nrange = round(10.^linspace(1, 6, 16));
+m = length(nrange);
+
+% Solution in binary64.
+for j = 1:m
+    n = nrange(j);
+    x_dp = a;
+    h_dp = (b-a)/n;
+    y_dp = y0;
+
+    for i=1:n
+        y_dp = Euler(true, decay_ODE_acc, h_dp, y_dp, 0, []);
+    end
+    efp(j, 1) = abs(y_dp - yexact);
+end
+
+options.round = 1; % RN
+
+% All chop formats.
+for k = 1:6
+    switch k
+      case 1, options.format = 'b'; ...
+            options.subnormal = 1; sr = 0;
+      case 2, options.format = 'b'; ...
+            options.subnormal = 1; sr = 1;
+      case 3, options.format = 'h'; ...
+            options.subnormal = 1; sr = 0;
+      case 4, options.format = 'h'; ...
+            options.subnormal = 1; sr = 1;
+      case 5, options.format = 's'; ...
+            options.subnormal = 1; sr = 0;
+      case 6, options.format = 's'; ...
+            options.subnormal = 1; sr = 1;
+    end
+
+    fprintf('k = %1.0f, prec = %s, subnormal = %1.0f\n', ...
+            k, options.format, options.subnormal)
+    chop([],options)
+
+    a = chop(a); b = chop(b); y0 = chop(y0);
+
+    for j = 1:m
+        n = nrange(j);
+        h_fp = chop((b-a)/n);
+        y_fp_init = chop(y0);
+
+        if sr
+            repeat = rep;
+        else
+            repeat = 1;
+        end
+
+        avg_err = 0;
+        max_err(j,k+1) = 0;
+        min_err(j,k+1) = Inf;
+
+        for l=1:repeat
+            y_fp = y_fp_init;
+            for i=1:n
+                y_fp = Euler(false, decay_ODE_inacc, h_fp, y_fp, ...
+                             sr, options);
+            end
+            err = abs(y_fp - yexact);
+            if err > max_err(j, k+1)
+                max_err(j, k+1) = err;
+            end
+            if err < min_err(j, k+1)
+                min_err(j, k+1) = err;
+            end
+            avg_err = avg_err + err;
+        end
+        efp(j, k+1) = avg_err/repeat;
+    end
+end
+
+fileName = sprintf('euler%d.dat', testcase);
+fileID = fopen(fileName, 'w');
+fprintf(fileID, '%s %s %s %s %s %s %s %s %s %s %s %s %s %s\n', 'fp64', ...
+        'bf16-rn','bf16-sr-avg', 'bf16-sr-worst', 'bf16-sr-best', ...
+        'fp16-rn','fp16-sr-avg', 'fp16-sr-worst', 'fp16-sr-best', ...
+        'fp32-rn', 'fp32-sr-avg', 'fp32-sr-worst', 'fp32-sr-best', 'n');
+fprintf(fileID, '%e %e %e %e %e %e %e %e %e %e %e %e %e %d\n', ...
+        [efp(:, 1:3), max_err(:,3), min_err(:,3), ...
+         efp(:, 4:5), max_err(:,5), min_err(:,5), ...
+         efp(:, 6:7), max_err(:,7), min_err(:,7), nrange']');
+fclose(fileID);
+
+function f = Euler(accurate, decay_ODE, h, y, SR, options)
+    if (accurate)
+        f = y + h*decay_ODE(y, options);
+    elseif (SR)
+        temp = options.round;
+        options.round = 5;
+        f = chop(y + ...
+                 chop(h*chop(decay_ODE(y, options), options), ...
+                      options), options);
+        options.round = temp;
+    else
+        f = chop(y + ...
+                 chop(h*chop(decay_ODE(y, options), options), ...
+                      options), options);
+    end
+end

README.md

@@ -0,0 +1,22 @@
+# Numerical MATLAB experiments with stochastic rounding
+This repository contains the source code for reproducing the results in [Sec. 8, 1]. The following scripts for generating the data presented in the survey are made available:
+   * [`test_sum.m`](./test_sum.m) - Experiment with the harmonic sum [Sec. 8a, 1].
+   * [`test_matvec.m`](./test_matvec.m) - Experiment with matrix-vector multiplication [Sec. 8a, 1].
+   * [`ODE_test.m`](./ODE_tests.m) - Solution of an exponential decay ODE with Euler's method [Sec. 8d, 1].
+   * [`unit_circle_ODE.m`](./unit_circle_ODE.m) - Solution of a unit circle ODE with Euler's method [Sec. 8d, 1].
+
+The script [`run_tests.m`](./run_tests.m) performs all of the tests from the paper [1] in one run.
+
+These scripts rely on the [chop](https://github.com/higham/chop) library for implementing low precision arithmetics with stochastic rounding and therefore it should be downloaded and placed on the MATLAB search path. These experiments were developed and run on MATLAB version 2021b.
+
+### References
+
+ [1] M. Croci, M. Fasi, N. J. Higham, T. Mary, and M. Mikaitis. [*Stochastic Rounding: Implementation, Error Analysis, and Applications*](http://eprints.maths.manchester.ac.uk/2836/). Technical Report 2021.17, Manchester Institute for Mathematical Sciences, The University of Manchester, UK, October 20201. Revised January 2022. To appear in R. Soc. Open Sci.
+
+### License
+
+This software is distributed under the terms of the 2-clause BSD software license (see [LICENSE](./LICENSE)).
+
+The MATLAB function `chop` is distributed under the terms of the [BSD 2-Clause "Simplified" License](https://raw.githubusercontent.com/higham/chop/master/license.txt).
+
+The CPFloat C library is distributed under the [GNU Lesser General Public License, Version 2.1 or later](https://raw.githubusercontent.com/mfasi/cpfloat/master/LICENSES/LGPL-2.1-or-later.txt).

run_tests.m

@@ -0,0 +1,38 @@
+% run_tests.m   Scripts to reproduce the data used in [1, Fig. 8.1--8.4].
+%
+% References:
+%   [1] M. Croci, M. Fasi, N. J. Higham, T. Mary, M. Mikaitis.
+%       Stochastic Rounding: Implementation, Error Analysis, and
+%       Applications. Tech. Report 2021.17, Manchester Institute for
+%       Mathematical Sciences, The University of Manchester, UK.
+%       October 2022. Revised January 2022.
+
+% Add dependencies
+addpath('./deps/chop');
+addpath('./deps/cpfloat/mex');
+currdir = pwd();
+cd('./deps/cpfloat/mex');
+cpfloat_compile_nomake;
+cd(currdir);
+
+% Run tests that produce Figure 8.1
+test_sum
+
+% Run tests that produce Figure 8.2
+test_matvec
+
+% Run tests that produce Figure 8.3
+for testcase = [0, 1]
+    ODE_tests
+end
+
+% Run tests that produce Figure 8.4
+format = 'bfloat16';
+for N = [32 512 2048 8192]
+    unit_circle_ODE
+end
+
+format = 'fp16';
+for N = [32 512 16384 65536]
+    unit_circle_ODE
+end

test_matvec.m

@@ -0,0 +1,68 @@
+% test_matvec.m   Test accuracy of matrix-vector products.
+%  This code uses the CPFloat library at https://github.com/mfasi/cpfloat.
+%
+% References:
+%   [1] M. P. Connolly, N. J. Higham, T. Mary. Stochastic rounding and its
+%       probabilistic backward error analysis. SIAM J. Sci. Comput., 43(1),
+%       pp. 566–585. February 2021. http://dx.doi.org/10.1137/20m1334796
+%
+%   [2] M. Croci, M. Fasi, N. J. Higham, T. Mary, M. Mikaitis.
+%       Stochastic Rounding: Implementation, Error Analysis, and
+%       Applications. Tech. Report 2021.17, Manchester Institute for
+%       Mathematical Sciences, The University of Manchester, UK.
+%       October 2022. Revised January 2022.
+
+clear all
+rng(1)
+
+fs = 14; ms = 7;
+m = 100;
+nlist = round(logspace(1,6,20));
+rep = 10;
+
+formats = ['b', 'h'];
+precisions = [8, 11];
+options.explim = 1;
+
+for k = 1:2
+    t = precisions(k);
+    u = 2^-t;
+    options.format = formats(k);
+    i = 0;
+    for n = nlist
+        i = i + 1;
+        fprintf('i = %2d, n = %7d\n',i,n);
+        A = 1e-3*rand(m,n);
+        x = rand(n,1);
+        ye = A*x;
+        absAx = abs(A)*abs(x);
+        options.round = 1;
+        y1 = matvec(A,x,options);
+        berr1(i) = max(abs(ye-y1)./absAx);
+        for j=1:rep
+            options.round = 5;
+            y2 = matvec(A,x,options);
+            berr2(i,j) = max(abs(ye-y2)./absAx);
+        end
+    end
+    berr2avg = sum(berr2,2)/rep;
+    berr2max = max(berr2,[],2);
+    berr2min = min(berr2,[],2);
+
+    filename = sprintf('RTN_vs_SR_t%d.dat', t);
+    fid = fopen(filename, 'w');
+    for i=1:length(nlist)
+        fprintf(fid, "%d %f %f %f %f %f %f\n", nlist(i), berr1(i), ...
+                berr2avg(i), berr2max(i), berr2min(i), min(1,nlist(i)*u), ...
+                min(1,sqrt(nlist(i))*u));
+    end
+    fclose(fid);
+end
+
+function y = matvec(A,x,options)
+    [m,n] = size(A);
+    y = zeros(m,1);
+    for i=1:n
+        y = cpfloat(y + cpfloat(A(:,i)*x(i), options), options);
+    end
+end