From 15aeca78f18db6211e0259dca7e745a611af1477 Mon Sep 17 00:00:00 2001 From: Mantas Mikaitis Date: Wed, 12 Jan 2022 11:17:35 +0000 Subject: [PATCH] Add code for reproducing experiments in the survey paper --- .gitmodules | 6 ++ LICENSE | 26 +++++++ ODE_tests.m | 163 +++++++++++++++++++++++++++++++++++++++++++ README.md | 22 ++++++ deps/chop | 1 + deps/cpfloat | 1 + run_tests.m | 38 ++++++++++ test_matvec.m | 68 ++++++++++++++++++ test_sum.m | 94 +++++++++++++++++++++++++ unit_circle_ODE.m | 172 ++++++++++++++++++++++++++++++++++++++++++++++ 10 files changed, 591 insertions(+) create mode 100644 .gitmodules create mode 100644 LICENSE create mode 100644 ODE_tests.m create mode 100644 README.md create mode 160000 deps/chop create mode 160000 deps/cpfloat create mode 100644 run_tests.m create mode 100644 test_matvec.m create mode 100644 test_sum.m create mode 100644 unit_circle_ODE.m diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000..00db3ed --- /dev/null +++ b/.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 diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..63538d2 --- /dev/null +++ b/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. diff --git a/ODE_tests.m b/ODE_tests.m new file mode 100644 index 0000000..aff6481 --- /dev/null +++ b/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 diff --git a/README.md b/README.md new file mode 100644 index 0000000..36ed23e --- /dev/null +++ b/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). diff --git a/deps/chop b/deps/chop new file mode 160000 index 0000000..74ab966 --- /dev/null +++ b/deps/chop @@ -0,0 +1 @@ +Subproject commit 74ab966a5cb1b1f4461012a464ef4d54e52797e0 diff --git a/deps/cpfloat b/deps/cpfloat new file mode 160000 index 0000000..3c38e61 --- /dev/null +++ b/deps/cpfloat @@ -0,0 +1 @@ +Subproject commit 3c38e615b267af2cf01d8a16b02f439ba0a27fe1 diff --git a/run_tests.m b/run_tests.m new file mode 100644 index 0000000..d70f2ae --- /dev/null +++ b/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 diff --git a/test_matvec.m b/test_matvec.m new file mode 100644 index 0000000..2991147 --- /dev/null +++ b/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 diff --git a/test_sum.m b/test_sum.m new file mode 100644 index 0000000..ca6d90b --- /dev/null +++ b/test_sum.m @@ -0,0 +1,94 @@ +% test_sum.m Test accuracy of partial sums of harmonic series. +% This script evaluates the harmonic sum \sum_{i=1}^{\infty} 1/i +% in different arithmetics and with different rounding modes. +% This code uses the function chop at https://github.com/higham/chop. +% +% References: +% [1] M. Fasi, M. Mikaitis. Algorithms for stochastically rounded +% elementary arithmetic operations in IEEE 754 floating-point +% arithmetic. IEEE Trans. Emerg. Topics Comput., 9(3), pp. 1451–1466. +% July 2021. +% http://dx.doi.org/10.1109/TETC.2021.3069165 +% +% [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. + +rng(4) + +formats = ['h', 'b']; +t = [11, 8]; % Number of precision bits. + +for global_index=1:2 + options.format = formats(global_index); + + % nvals = [10 50 100 5e2 1000 5e3 1e4 5e4 1e5 5e5 1e6]; + nmax = 1e6; + npoints = 10; + nvals = round( exp( linspace(log(10), log(nmax), npoints ) ) ); + + clear err + err = zeros(length(nvals), 4); %RN, RS avg, RS worst, RS best + + % alpha = 0.9; + alpha = 1; + + for i = 1:length(nvals) + n = nvals(i); + x = chop(ones(1,n) ./ (1:n).^alpha); + sx = sum(x); + for k = [1 5] + + options.round = k; + % Initialize: subsequent calls chop(x) reuse options. + chop([],options) + + if (k==5) + rep = 10; + else + rep = 1; + end + errors = zeros(1,rep); + for repeat = 1:rep + s = 0; + for j = 1:n + s = chop(s + x(j)); + end + errors(repeat) = abs((sx-s)/sx); + end + eavg = sum(errors)/rep; + eworst = max(errors); + ebest = min(errors); + err(i,round(k/5)+1) = eavg; + err(i,3) = eworst; + err(i,4) = ebest; + end + + end + nearmax = err(:,1); + stochmax = err(:,2); + stochworst = err(:,3); + stochbest = err(:,4); + + fs=12; ms=7; lw=1; + + dim = nvals; + lam = 1; + u = 2^(-11); % unit roundoff for fp16 + pbound = exp(lam*sqrt(dim).*(2*u) + dim.*(2*u)^2/(1-2*u)) - 1; + pbound(pbound < 0) = 1; pbound(pbound > 1) = 1; + + dbound = (dim.*u)./(1-dim.*u); + dbound(dbound < 0) = 1; dbound(dbound > 1) = 1; + + filename = sprintf('test_sum_t%d.dat', t(global_index)); + fid = fopen(filename, 'w'); + for i=1:length(dim) + fprintf(fid, "%d %f %f %f %f %f %f\n",... + dim(i), nearmax(i), stochmax(i), stochworst(i), ... + stochbest(i), dbound(i), pbound(i)); + end + fclose(fid); +end \ No newline at end of file diff --git a/unit_circle_ODE.m b/unit_circle_ODE.m new file mode 100644 index 0000000..153c1ea --- /dev/null +++ b/unit_circle_ODE.m @@ -0,0 +1,172 @@ +% unit_circle_ODE.m Tests for Euler integration of unit circle ODE. +% This script integrates the unit circle ODE system +% u'(t) = v(t), u(0) = 1, +% v'(t) = u(t), v(0) = 0, +% using the forward Euler method with various floating-point arithmetics. +% This code uses the function chop at https://github.com/higham/chop. +% +% References: +% [1] N. J. Higham (ed.). Princeton Companion to Applied Mathematics. +% p. 51. 2015. +% +% [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), pp. 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 2021. Revised January 2022. + +% Clear chop options and reset the PRNG seed. +clear coordinates_exact coordinates coordinates_temp +rng(500) + +% Set the number of times to repeat the SR experiments. +rep = 10; + +% Steps of ODE integration. +if ~exist('N', 'var') + N = 65536; +end + +% Target format. +if ~exist('format', 'var') + format = 'fp16'; +end + +% Set up number of coordinates to sample. +points = min(100, N); + +% Arrays for coordinates. +coordinates_exact = zeros(2, points+2); +coordinates = zeros(4, points+2, 2); +coordinates_temp = zeros(2, points+2, 2, rep); + +% Get the coordinates of the exact solution. +h_dp = 2*pi/N; +j = 1; +coordinates_exact(:, 1) = [1,0]; +for i=1:N + if (mod(i, floor(N/points)) == 0) + j = j+1; + coordinates_exact(:, j) = [cos((i-1)*h_dp), -sin((i-1)*h_dp)]; + end +end +coordinates_exact(:, end) = [1,0]; + +options.round = 1; % RN +% Run two cases: with RN and with SR. +for k = 1:2 + switch k + case 1, options.format = format; options.subnormal = 1; sr = 0; + case 2, options.format = format; options.subnormal = 1; sr = 1; + end + + fprintf('k = %1.0f, prec = %s, subnormal = %1.0f\n',... + k,options.format,options.subnormal) + chop([],options) + + % Timestep size + h_fp = chop(2*pi/N, options); + closest_avg = 999; + closest_avg_i = 0; + furthest_avg_outside = 0; + furthest_avg_outside_i = 0; + furthest_avg_inside = 0; + furthest_avg_inside_i = 0; + avg_err = 0; + avg_err_signed = 0; + if sr + repeat = rep; + else + repeat = 1; + end + for l=1:repeat + j=1; avg_err = 0; + % Initial values. + u_fp = chop(1); + v_fp = chop(0); + coordinates_temp(k, j, :, l) = [u_fp, v_fp]; + for i=1:N + [u_fp, v_fp] = Euler(0, h_fp, u_fp, v_fp, sr, options); + if (mod(i, floor(N/points)) == 0 && j < points+1) + j = j+1; + coordinates_temp(k, j, :, l) = [u_fp, v_fp]; + end + end + coordinates_temp(k, end, :, l) = [u_fp, v_fp]; + end + if (sr) + % Compute closest and farthest points. + u_pos = squeeze(coordinates_temp(k,:,1,:)); + v_pos = squeeze(coordinates_temp(k,:,2,:)); + norms = sqrt(u_pos.^2 + v_pos.^2); + [~, minpos] = min(norms, [], 2); + [~, maxpos] = max(norms, [], 2); + for i = 1:points+2 + coordinates(2, i, 1) = u_pos(i, maxpos(i)); + coordinates(2, i, 2) = v_pos(i, maxpos(i)); + coordinates(3, i, 1) = u_pos(i, minpos(i)); + coordinates(3, i, 2) = v_pos(i, minpos(i)); + end + coordinates(4, :, :) = squeeze( ... + sum(coordinates_temp(k, :, :, :), 4)) / repeat; + else + coordinates(1, :, :) = squeeze( ... + sum(coordinates_temp(1, :, :, :), 4)) / repeat; + end +end + +fileNamePrefix = sprintf('unit_circle_%s', options.format); +fileNameSuffix = sprintf('_%d.dat', N); + +fileID = fopen('unit_circle_exact.dat', 'w'); +fprintf(fileID, 'u v \n'); +fprintf(fileID, '%e %e \n', [coordinates_exact(1,:); ... + coordinates_exact(2,:)]); +fclose(fileID); + +filename = sprintf('%s_rn%s', fileNamePrefix, fileNameSuffix); +fileID = fopen(filename, 'w'); +fprintf(fileID, 'u v \n'); +fprintf(fileID, '%e %e \n', [coordinates(1,:,1); coordinates(1,:,2)]); +fclose(fileID); + +filename = sprintf('%s_sr_worst_outside%s', fileNamePrefix, ... + fileNameSuffix); +fileID = fopen(filename, 'w'); +fprintf(fileID, 'u v \n'); +fprintf(fileID, '%e %e \n', [coordinates(2,:,1); coordinates(2,:,2)]); +fclose(fileID); + +filename = sprintf('%s_sr_worst_inside%s', fileNamePrefix, fileNameSuffix); +fileID = fopen(filename, 'w'); +fprintf(fileID, 'u v \n'); +fprintf(fileID, '%e %e \n', [coordinates(3,:,1); coordinates(3,:,2)]); +fclose(fileID); + +filename = sprintf('%s_sr_best%s', fileNamePrefix, fileNameSuffix); +fileID = fopen(filename, 'w'); +fprintf(fileID, 'u v \n'); +fprintf(fileID, '%e %e \n', [coordinates(4,:,1); coordinates(4,:,2)]); +fclose(fileID); + +function [u, v] = Euler(accurate, h, u_prev, v_prev, SR, options) + if (accurate) + u = u_prev + h*v_prev; + v = v_prev - h*u_prev; + elseif (SR) + temp = options.round; + options.round = 5; + u = chop(u_prev + chop(h * v_prev, options), options); + v = chop(v_prev - chop(h * u_prev, options), options); + options.round = temp; + else + u = chop(u_prev + chop(h * v_prev, options), options); + v = chop(v_prev - chop(h * u_prev, options), options); + end +end \ No newline at end of file