# 2007-11-27: Comparison of numerical results on various Matlab
# versions, CPUs and operating systems. The test-case uses pdepe
# and fft, besides some basic array manipulation.
#
# OS=gentoo means linux kernel >2.6.20, glibc 2.5, gcc 4.1.x
# OS=rhel5 means RedHat Enterprise 5: kernel 2.6.18.x, glibc 2.5.
# 
# 64|32 bit refers to the kernel's support for 64bit programs,
# and if these are supported then the OS is a 64bit-capable one
# and the 64bit Matlab is used.
#
# A line of shell, used to run the example file:
# for n in 5 4 3 1 ; do 
#      matlab-7.$n -nojvm -nosplash -r 'ml_numres_example; exit'
# done


# First, ordered by system and versions.
# CPU	  64|32bit    OS     MLvers     output

Athlon-XP   32       gentoo   7.5       2.095668e-11, 4.166408e-15

Xeon-Core2  64       gentoo   7.5       2.095668e-11, 4.166408e-15
Xeon-Core2  64       gentoo   7.3       2.095661e-11, 4.100776e-15
Xeon-Core2  64       gentoo   7.4       2.095661e-11, 4.100776e-15
Xeon-Core2  64       gentoo   7.1       2.095664e-11, 4.154064e-15

Xeon-P4     64       gentoo   7.5       2.095668e-11, 4.166408e-15
Xeon-P4     64       gentoo   7.4       2.095661e-11, 4.100776e-15
Xeon-P4     64       gentoo   7.3       2.095661e-11, 4.100776e-15
Xeon-P4     64       gentoo   7.1       2.095664e-11, 4.154064e-15

Athlon64-X2 64       gentoo   7.5       2.095624e-11, 3.829428e-15
Athlon64-X2 64       gentoo   7.4       2.095661e-11, 4.100776e-15
Athlon64-X2 64       gentoo   7.3       2.095661e-11, 4.100776e-15
Athlon64-X2 64       gentoo   7.1       2.095664e-11, 4.154064e-15

Opteron     64       rhel5    7.5       2.095624e-11, 3.829428e-15
Opteron     64       rhel5    7.4       2.095661e-11, 4.100776e-15
Opteron     64       rhel5    7.3       2.095661e-11, 4.100776e-15
Opteron     64       rhel5    7.1       2.095664e-11, 4.154064e-15

Pentium4    32       gentoo   7.5       2.095668e-11, 4.166408e-15
Pentium4    32       gentoo   7.4       2.095661e-11, 4.100776e-15
Pentium4    32       gentoo   7.3       2.095661e-11, 4.100776e-15
Pentium4    32       gentoo   7.1       2.095664e-11, 4.154064e-15

Athlon64    64       rhel5    7.5       2.095624e-11, 3.829428e-15
Athlon64    64       rhel5    7.4       2.095661e-11, 4.100776e-15
Athlon64    64       rhel5    7.3       2.095661e-11, 4.100776e-15
Athlon64    64       rhel5    7.1       2.095664e-11, 4.154064e-15


# Now, the same lines, sorted on the last field to group similar
# numerical behaviours.
# CPU	  64|32bit    OS     MLvers     output

Athlon64    64       rhel5    7.5       2.095624e-11, 3.829428e-15
Athlon64-X2 64       gentoo   7.5       2.095624e-11, 3.829428e-15
Opteron     64       rhel5    7.5       2.095624e-11, 3.829428e-15

Athlon64    64       rhel5    7.3       2.095661e-11, 4.100776e-15
Athlon64    64       rhel5    7.4       2.095661e-11, 4.100776e-15
Athlon64-X2 64       gentoo   7.3       2.095661e-11, 4.100776e-15
Athlon64-X2 64       gentoo   7.4       2.095661e-11, 4.100776e-15
Opteron     64       rhel5    7.3       2.095661e-11, 4.100776e-15
Opteron     64       rhel5    7.4       2.095661e-11, 4.100776e-15
Pentium4    32       gentoo   7.3       2.095661e-11, 4.100776e-15
Pentium4    32       gentoo   7.4       2.095661e-11, 4.100776e-15
Xeon-Core2  64       gentoo   7.3       2.095661e-11, 4.100776e-15
Xeon-Core2  64       gentoo   7.4       2.095661e-11, 4.100776e-15
Xeon-P4     64       gentoo   7.3       2.095661e-11, 4.100776e-15
Xeon-P4     64       gentoo   7.4       2.095661e-11, 4.100776e-15

Athlon64    64       rhel5    7.1       2.095664e-11, 4.154064e-15
Athlon64-X2 64       gentoo   7.1       2.095664e-11, 4.154064e-15
Opteron     64       rhel5    7.1       2.095664e-11, 4.154064e-15
Pentium4    32       gentoo   7.1       2.095664e-11, 4.154064e-15
Xeon-Core2  64       gentoo   7.1       2.095664e-11, 4.154064e-15
Xeon-P4     64       gentoo   7.1       2.095664e-11, 4.154064e-15

Athlon-XP   32       gentoo   7.5       2.095668e-11, 4.166408e-15
Pentium4    32       gentoo   7.5       2.095668e-11, 4.166408e-15
Xeon-Core2  64       gentoo   7.5       2.095668e-11, 4.166408e-15
Xeon-P4     64       gentoo   7.5       2.095668e-11, 4.166408e-15


# `We' already knew that numerical differences *can* arise due 
# to different compiler options, or different code being chosen
# for different CPUs.  It was assumed that Matlab would be 
# designed to provide similar results for its own functions, 
# independent of the CPU type, even at expense of some optimisation.
#
# Matlab uses FFTW: this is known for doing clever things to increase
# speed but lose some accuracy. Is it *just* the FFT part that's 
# making the difference? 
# Tried looking at just the time-domain values (i_td) that are an 
# intermediate step in the calculation. These also revealed 
# differences, with the same pair being similar for the one host
# tried (`Opteron').
# 
# The groupings of results show that there is no consistent effect of
# a single parameter of CPU, 64|32 or OS, on the results.
# Only Matlab version 7.5 showed variations between different systems; 
# this divided the systems into AMD's 64bit Opteron and Athlon64, 
# against any Intel or the older 32bit Athlon-XP.
# Matlab versions 7.3 and 7.4 gave the same results as each other, 
# consistent over all systems tried.
# Matlab version 7.1 was consistent between different systems, but 
# differed from all of the other versions.
#
# This particular case isn't a practical problem. More, I'm puzzled
# by whether one really should expect such variation.

# end

# The example script.
cat <<EOT
% Example of difference between calculation on Intel and AMD hardware.
% First two columns of output are the same for any system (determined
% by input).
%        40000,  3.0000e-05,     2.095624e-11,   3.829428e-15
%        40000,  3.0000e-05,     2.095668e-11,   4.166408e-15
function Ccomp = ml_numres_example()
tstart = cputime;
j = sqrt(-1);
Nt = 3 * 32;
f = 0.000188495559215388 / ( 2 * pi );
Tc = 1/f;
Ts = Tc/32;
maxstep = 1/32 * Ts;
options = odeset( 'MaxStep', maxstep, ...
    'RelTol', 10e-6, ...
    'AbsTol', 10e-6  );
t_td = NaN*zeros( Nt, 1 );
v_td = t_td;
i_td = t_td;
f_fd = NaN*zeros( 8, 1 );
V_fd = f_fd;
I_fd = f_fd;
dashedline = '--------------------------------------------------------------';
disp(dashedline)
tic; cpustart = cputime;
x = linspace( 0, 0.08, 32 );
t_all = ((-2*32):(Nt-1)) * Ts;
soln_all = pdepe( 0, @nlc_pde, @nlc_ic, @nlc_bc, x, t_all, options );
t_td = t_all( (2*32+1):end );
soln = soln_all( (2*32+1):end, : );
fprintf( '    PDE solution took %5.2fs real time, %5.2fs cpu time\n', ...
    toc, cputime-cpustart );
tic; cpustart = cputime;
v_td = zeros( Nt, 1 );
dudx = zeros( Nt, 1 );
for tidx=1:Nt,
    [ v_td(tidx), dudx(tidx) ] = pdeval( 0, x, soln(tidx,:), x(1) );
end
i_td = - dudx .* 1/1.8182e+12 .* exp( 0.00115 * abs(dudx).^0.66666 );
fft_i = fft( i_td );
I_fd = 2/length(i_td) * fft_i( 1 + (1:8)*3 );
I_fd_dc = fft_i(1)/length(i_td);
fft_v = fft( v_td );
V_fd = 2/length(v_td) * fft_v( 1 + (1:8)*3 );
V_fd_dc = fft_v(1)/length(v_td);
f_fd = f*(1:8).';
Ccomp =  I_fd ./ ( j*(0.000188495559215388*(1:8).') .* V_fd );
fprintf( '    Postprocessing took %5.3fs real time, %5.3fs cpu time\n', ...
    toc, cputime-cpustart );
fprintf('    Total CPU time: %5.3fs\n', cputime-tstart );
disp(dashedline)
fprintf( '\t%5d,\t%1.4e, \t%1.6e,\t%1.6e\n', ...
    round(40000), f, real(Ccomp(1)), -imag(Ccomp(1)) );
fout = fopen('out_itd.m','a');
fprintf(fout, 'i_td = [ ');
fprintf(fout, '%1.6e ', i_td); 
fprintf(fout, '] \n');
disp(dashedline)
function [ c, f, s ] = nlc_pde( x, t, u, dudx )
s = 0;
c = 2.6625e-10;
f = dudx .* 1/1.8182e+12 .* exp( 0.00115 * abs(dudx).^0.66666 );
function u0 = nlc_ic(x);
u0 = 0*x;
function [ pl, ql, pr, qr ] = nlc_bc( xl, ul, xr, ur, t )
pl = ul - v_supply(t);
ql = 0;
pr = 0*ur;
qr = 1;
function vs = v_supply( t )
vs = 40000 * sin( 0.000188495559215388*t + 0 );
EOT