// Copyright (C) 2010 - DIGITEO - 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

// Benchmarking backslash

// References
// Cleve Moler, "Benchmarks: LINPACK and MATLAB - Fame and fortune from megaflops", 1994, http://www.mathworks.com/company/newsletters/news_notes/pdf/sumfall94cleve.pdf
// "Programming in Scilab", Michael Baudin, 2010, http://forge.scilab.org/index.php/p/docprogscilab/downloads/

/////////////////////////////////////////////////////////////////
//
// A basic benchmark
if ( %f ) then
stacksize("max")
rand( "normal" );
n = 1000;
A = rand(n,n);
b = rand(n,1);
tic();
x = A\b;
t = toc();
mflops = (2/3*n^3 + 2*n^2)/t/1.e6;
disp([n t mflops])
end

/////////////////////////////////////////////////////////////////
//
// A subtle benchmark
scf();

stacksize("max")
rand( "normal" );
s = stacksize();
MB = round(s(1)*8/10^6);
nmax = sqrt(s(1));
mprintf("Memory: %d (MB)\n",MB);
mprintf("Maximum n: %d\n",nmax);

timemin = 0.1;
timemax = 8.0;
nfact = 1.2;

//
// Make a loop over n
xtitle("Backslash","Matrix Order (n)","Megaflops");
n = 1;
k=1;
perftable = [];
while ( %t )
  A = rand(n,n);
  b = rand(n,1);
  c = ceil(log10(1/rcond(A)));
  if ( c > 7 ) then
    // Avoid least square computation.
    // See : http://bugzilla.scilab.org/show_bug.cgi?id=7487
    mprintf("Skipping ill-conditionned matrix: c=%d\n",c);
    continue;
  end
  tic();
  ierr = execstr("x = A\b;","errcatch");
  t = toc();
  if ( ierr <> 0 ) then
    laerr = lasterror();
    disp("Error:");
    disp(laerr);
    break
  end
  if ( t > timemin ) then
    mflops = (2/3*n^3 + 2*n^2)/t/1.e6;
    perftable(k,:) = [n t mflops];
    mprintf("Run #%d: n=%6d, T=%.3f (s), Mflops=%6d, Log(cond)=%2d\n",k,perftable(k,1),perftable(k,2),perftable(k,3),c)
    k = k+1;
    plot(n,mflops,"bo-")
  end
  if ( t > timemax ) then
    break
  end
  n = ceil(nfact * n);
end

// Search for best performance
[M,k] = max(perftable(:,3));
mprintf("Best performance:")
mprintf(" N=%d, T=%.3f, MFLOPS=%d\n",perftable(k,1),perftable(k,2),perftable(k,3));