bm_version = ["bm ", "2.3"];

# Benchmark for octave, released under the GNU GPL, version 2 or later
# Francesco Potort́ <Potorti@isti.cnr.it>  1996, 1997, 2005
# 2005/11/24 14:44:53
#
# latest version at <ftp://fly.isti/pub/software/octave/benchmark.m>
# old results list at <ftp://fly.isti.cnr.it/pub/software/octave/bm_results>

printf ("Octave benchmark version %s\n", bm_version);

# To add reference times for your machine run the benchmark and
# add the values contained in the bm_mytime vector.
#
bm_refname = "Pentium II 350 MHz"; # with Debian's lapack3 and refblas3
bm_reftime = [1.231  3.54  0.885  0.792  1.135];
# Using Debian's atlas3-base: [1.77 1.38 .857 .778 1.12]

# Use clock() if cputime() does not work on this particular port of octave.
# In this case, time will be computed on a wall clock, and will make sense
# only on a machine where no other processes are consuming significant cpu
# time while the benchmark is running.
global bm_uses_cputime = (cputime() != 0);
if (!bm_uses_cputime)
  disp ...
    ("WARNING: if other processes are running the figures will be inaccurate");
endif
function t = bm_start ()
  global bm_uses_cputime
  if (bm_uses_cputime)
    t = cputime();
  else
    t  = clock();
  endif
endfunction
function et = bm_stop (t);
  global bm_uses_cputime
  if (bm_uses_cputime)
    et = cputime()-t;
  else
    et = etime(clock(),t);
  endif
endfunction

# Used for the lsode test.
clear xdot
function xdot = xdot (x, t)
  r = 0.25; k = 1.4; a = 1.5; b = 0.16; c = 0.9; d = 0.8;
  xdot(1) = r*x(1)*(1 - x(1)/k) - a*x(1)*x(2)/(1 + b*x(1));
  xdot(2) = c*a*x(1)*x(2)/(1 + b*x(1)) - d*x(2);
endfunction

#
# Do benchmark
#
function [name, time] = bm_test(f,rep)		# Actual test functions
  global t;
  start = bm_start();
  for i = 1:rep
    if     (f==1) name="Matrix inversion (LAPACK)";
                  bm_x=inv(hadamard(512));
    elseif (f==2) name="Schur decomposition (LAPACK)";
                  bm_x=schur(hadamard(256));
    elseif (f==3) name="Differential equation (LSODE)";
                  bm_x=lsode("xdot",[1;2],(t=linspace(0,50,400)'));
    elseif (f==4) name="Fourier transforms (FFTPACK)";
                  bm_x=ifft2(fft2(hadamard(512)));
    elseif (f==5) name="for loop";
                  for i=1:20000;bm_x=i^2;endfor
    endif
  endfor
  time = bm_stop(start)/rep;
endfunction

bm_targetaccuracy = 0.025;		# target accuracy of mean of times
bm_minrepetitions = 7;			# min number of repetitions per test
bm_maxtime = 60;			# max runtime per test [seconds]
bm_mintime = 0.3;			# min runtime per test [seconds]
bm_runtime = 3;				# target runtime per test [seconds]

printf ("Speed of octave %s on %s\n...
  %s has score 1.0, higher numbers are better\n", ...
        version(), computer(), bm_refname); fflush(stdout);
bm_mytime = zeros(size(bm_reftime));
for f = 1:length(bm_reftime)
  res = [];
  bm_test(f,1);				# increase the RSS, load things
  rep = 1;				# number of repetitions per run
  while (1)				# we would need a do..while really
    [name,time] = bm_test(f,rep);	# evaluate name and time
    if (time*rep > bm_mintime)		# run for at least bm_mintime
      break;				# found approximate time
    endif
    rep = 2*rep;			# approaching min run time
  endwhile
  printf("%-33s", name); fflush(stdout);# print name
  rep = round(bm_runtime/time);		# no. of repetitions per run
  rep = max(1,rep);			# slow machines need this
  for runs = 1:bm_maxtime/bm_runtime	# do runs
    [name,time] = bm_test(f,rep);	# run
    res(runs) = bm_reftime(f)/time;	# store relative performance
    if (runs < bm_minrepetitions)	# jump rest of for loop
      continue
    endif
    res = sort(res);
    bm_mean = mean(res(2:runs-1));	# remove min and max results
    if (std(res)/bm_mean < bm_targetaccuracy)
      break
    endif
  endfor				# end of repetitions loop
  bm_mytime(f) = bm_reftime(f)/bm_mean;
  # print 95% confidence interval
  printf("%5.2f +/- %.1f%% (%d runs)\n", ...
         bm_mean, 200*std(res)/bm_mean, runs*rep); fflush(stdout);
endfor
clear bm_x
# Display the geometric mean of the results
printf("-- Performance index (%s): %.2g\n\n", bm_version, ...
       prod(bm_reftime./bm_mytime)^(1/length(bm_reftime)));