// 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 the matrix-matrix product

// References
// "Programming in Scilab", Michael Baudin, 2010, http://forge.scilab.org/index.php/p/docprogscilab/downloads/

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

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

stacksize("max");
rand( "normal" );
lines(0);

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("Matrix-Matrix multiply","Matrix Order (n)","Megaflops");
n = 1;
k = 1;
perftable = [];
while ( %t )
  A = rand(n,n);
  B = rand(n,n);
  tic();
  ierr = execstr("C = A * B","errcatch");
  t = toc();
  if ( ierr <> 0 ) then
    laerr = lasterror();
    disp("Error:");
    disp(laerr);
    break
  end
  if ( t > timemin ) then
    mflops = 2*n^3/t/1.e6;
    perftable(k,:) = [n t mflops];
    plot(n,mflops,"bo-")
    mprintf("Run #%d: n=%6d, T=%.3f (s), Mflops=%6d\n",k,perftable(k,1),perftable(k,2),perftable(k,3))
    k = k+1;
  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 (s), MFLOPS=%d\n",perftable(k,1),perftable(k,2),perftable(k,3));