stacksize(80000000)

cd "C:\Program Files\scilab-4.1.2\intfftw";
exec "loader.sce";

disp('The following benchmark program will print the average timings');
disp('to calculate the functions by 3 runs.');
disp(' ');
disp(' ');
disp('!!! SCILAB - Benchmarkprogram !!!');
disp('=================================');
disp(' ');
 
//* Test of IO function and descriptive statistics *

result = 0;
a = 0;
b = 0;
runs = 3;
for i=1:runs
   tic();
   [fd,err]=mopen('D:\Mathematik\Ncrunch5\currency2.txt','r',1);
   header=mgetstr(154,fd);
   datmat=fscanfMat('D:\Mathematik\Ncrunch5\currency2.txt');
   LimitsArea={1,261,522,784,1045,1305,1565,1827,2088,2349,2610,2871,3131,3158};
  for k=1:2000
    for j=1:13
      Year_Data=datmat(LimitsArea(j):LimitsArea(j+1)-1,5:38);
      Min_Year=min(Year_Data,1);
      Max_Year=max(Year_Data,1);
      Mean_Year=mean(Year_Data,1);
      GainLoss_Year=100-(Year_Data(1,:)+0.00000001)./(Year_Data($,:)+0.00000001)*100;
    end
  end
   zeit=toc();
   result=result+zeit;
end
result=result/runs;
disp('IO test and descriptive statistics__________________: '+string(result)+' sec.');


//* Misc. operation *
 
result = 0;
a=1;
for i = 1:runs
  tic();
  for x = 1:15000
    for y = 1:15000
      a=a+x+y;
    end
  end
  zeit = toc();
  result = result+zeit;
end
result = result/runs;
disp('Loop test____________________________________________: '+string(result)+' sec.');
clear('a');
clear('b');
for i = 1:runs
  for j=1:200
    b = 1+(rand(2000,2000,'uniform')/100);
    tic();
    a = b.^1000;
    result = result+toc();
  end
end
result = result/runs;
disp('2000x2000 normal distributed random matrix^1000______: '+string(result)+' sec.');
clear('a');
clear('b');
result = 0;
for i = 1:runs
  for j=1:100
    a = rand(1000000,1,'uniform');
    tic();
    b = sort(a);
    result = result+toc();
  end
end
result = result/runs;
disp('1000000 values sorted ascending______________________: '+string(result)+' sec.');
 
//* Analysis *
 
clear('a');
clear('b');
result = 0;
for i = 1:runs
  for j=1:100
    a = rand(1048576,1,'uniform');
    tic();
    b = fftw(a);
    result = result+toc();    
  end
end
result = result/runs;
disp('FFT over 1048576 values______________________________: '+string(result)+' sec.');
 
//* Algebra *
 
clear('a');
clear('b');
result = 0;
for i = 1:runs
  for j=1:100
    a = rand(1500,1500,'uniform');
    tic();
    b = det(a);
    result = result+toc();
  end
end
result = result/runs;
disp('Determinant of a 1500x1500 random matrix_____________: '+string(result)+' sec.');
clear('a');
clear('b');
result = 0;
for i = 1:runs
  for j=1:100
    a = rand(1500,1500,'uniform');
    tic();
    b = inv(a);
    result = result+toc();
  end
end
result = result/runs;
disp('Inverse of a 1500x1500 uniform distr. random matrix__: '+string(result)+' sec.');
clear('a');
clear('b');
result = 0;
for i = 1:runs
  a = rand(1200,1200,'uniform');
  tic();
  c = spec(a);
  result = result+toc();
end
result = result/runs;
disp('Eigenval. of a normal distr. 1200x1200 randommatrix__: '+string(result)+' sec.');
clear('a');
clear('b');
result = 0;
for i = 1:runs
  for j=1:100
    a = rand(1500,1500,'uniform');
    a = a'*a;
    tic();
    b = chol(a);
    result = result+toc();
  end
end
result = result/runs;
disp('Cholesky decomposition of a 1500x1500-matrix_________: '+string(result)+' sec.');
clear('a');
clear('b');
result = 0;
for i = 1:runs
  for j=1:100
    a = rand(1500,1500,'uniform');
    tic();
    b = a'*a;
    result = result+toc();    
  end
end
result = result/runs;
disp('1500x1500 cross-product matrix_______________________: '+string(result)+' sec.');
 
//* Number theory *
 
clear('a');
clear('b');
phi = 1.61803398874989;
result = 0;
for i = 1:runs
  a = floor(1000*rand(10000000,1,'uniform'));
  tic();
  b = (phi.^a-(-phi).^(-a))/sqrt(5);
  result = result+toc();
end
result = result/runs;
disp('Calculation of 10000000 fibonacci numbers_____________: '+string(result)+' sec.');

clear('a');
clear('b');
// * pca without output */
function [lambda,facpr,comprinc]=pca2(varargin)

//This  function performs  several  computations known  as
//"principal component  analysis". It includes  drawing of
//"correlations circle",  i.e. in the  horizontal axis the
//correlation   values  r(c1;xj)   and  in   the  vertical
//r(c2;xj). It is an extension of the pca function.
//
//The  idea  behind this  method  is  to  represent in  an
//approximative  manner a  cluster of  n individuals  in a
//smaller  dimensional subspace.  In order  to do  that it
//projects the cluster onto a subspace.  The choice of the
//k-dimensional projection subspace is  made in such a way
//that  the distances  in  the projection  have a  minimal
//deformation: we are looking for a k-dimensional subspace
//such that the squares of the distances in the projection
//is  as  big  as  possible  (in  fact  in  a  projection,
//distances can only stretch).  In other words, inertia of
//the projection  onto the k dimensional  subspace must be
//maximal.
//
//x is a nxp (n  individuals, p variables) real matrix.
//
//lambda is a px2 numerical matrix. In the first column we
//find the  eigenvalues of V,  where V is  the correlation
//pxp matrix  and in the  second column are the  ratios of
//the   corresponding   eigenvalue   over   the   sum   of
//eigenvalues.
//
//facpr  are  the principal  factors:  eigenvectors of  V.
//Each  column is an  eigenvector element  of the  dual of
//R^p. Is an orthogonal matrix.
//
//comprinc  are  the  principal components.   Each  column
//(c_i=Xu_i)  of  this  nxn  matrix  is  the  M-orthogonal
//projection of individuals onto principal axis.  Each one
//of this columns is a linear combination of the variables
//x1,  ...,xp   with  maximum  variance   under  condition
//u'_iM^(-1)u_i=1.
//
//Verification: comprinc*facpr=x
//
//References: Saporta, Gilbert, Probabilites,  Analyse des
//Donnees et Statistique, Editions Technip, Paris, 1990.
//
//author: carlos klimann
//
//date: 2002-02-05
//commentary fixed 2003-19-24 ??
// update disable graphics output Allan CORNET 2008

  [lhs,rhs]=argn(0)
 
  x = [];
  N = [];
  no_graphics = %f;
  
  select rhs,
  	case 1 then
  	  x = varargin(1);
  	  N=[1 2];
  	case 2 then 
  	  x = varargin(1);
  	  N = varargin(2);
  	case 3 then
  	  x = varargin(1);
  	  N = varargin(2);
  	  no_graphics = varargin(3);
  	else
  	error(sprintf("%s: Wrong number of input argument(s).\n","pca"));
  	end

  if ( type(no_graphics) <> 4 ) then
    error(sprintf("%s: Wrong type for third input argument: Boolean expected.\n" ,"pca"));
  end

  if x==[] then
    lambda=%nan;
    facpr=%nan;
    comprinc=%nan; 
    return; 
  end

  [rowx colx]=size(x)
  if size(N,'*')<>2 then error('Second parameter has bad dimensions'), end,
  if max(N)>colx then disp('Graph demand out of bounds'), return, end
  
  xbar=sum(x,'r')/rowx
  y=x-ones(rowx,1)*xbar
  std=(sum(y .^2,'r')) .^ .5 
  V=(y'*y) ./ (std'*std)
  [lambda facpr]=bdiag(V,(1/%eps))
  [lambda order]=sort(diag(lambda))
  lambda(:,2)=lambda(:,1)/sum(lambda(:,1))
  facpr=facpr(:,order)
  comprinc=x*facpr
  
  if ~no_graphics then
    w = winsid();
    if (w == []) then 
      w = 0; 
    else
      w = max(w)+1;
    end
    fig1 = scf(w);
    
    rc = (ones(colx,1)* sqrt((lambda(N,1))')) .* facpr(:,order(N)) ;
    rango = rank(V);
    ra = [1:rango]';
    if ( rango <= 1 ) then return, end
    plot2d(-rc(ra,1),rc(ra,2),style=-10);
    legend('(r(c1,xj),r(c2,xj)');
    ax=gca();
    ax.x_location="middle";
    ax.y_location = "middle";
    blue=color('blue')
    for k=1:rango,
      xstring(rc(k,1),rc(k,2),'X'+string(k));
      e=gce();
      e.foreground=blue;
    end
    title(' -Correlations Circle- ');
    fig2 = scf(w+1);
    plot2d3([0;ra;rango+1]',[0; lambda(ra,2);0]);
    plot2d(ra,lambda(ra,2),style=9);
    ax=gca();
    ax.grid=[31 31];
    plot2d3([0;ra;rango+1]',[0; lambda(ra,2);0])
    plot2d(ra,lambda(ra,2),style=9)
    for k=1:rango,
      xstring(k,0,'l'+string(k)),
      e=gce();e.font_style=1
    end
    title(' -Eigenvalues (in percent)- ')
    ylabel('%')
  end  
endfunction
// * END pca without output */


result=0; 

for i=1:runs;
    a=rand(10000,1000,'uniform');
    tic();
    b=pca2(a,[1 2],%t);
    result=result+toc();
end;


result=result/runs;
disp(['Calc. of the principal components of a 10000x1000 matrix : '+string(result)+' sec.'])
clear pca_without_graphics;
 
//* Stochastic-statistic *
 
clear('a');
clear('b');
result = 0;
for i = 1:runs
  for j=1:100
    a = rand(1500,1500,'uniform');
    tic();
    b = gamma(a(:));
    result = result+toc();
  end
end
result = result/runs;
disp('Gamma function over a 1500x1500 matrix_______________: '+string(result)+' sec.');
clear('a');
clear('b');
result = 0;
for i = 1:runs
  for j=1:100
    a = rand(1500,1500,'uniform');
    tic();
    b = erf(a(:));
    result = result+toc();    
  end
end
result = result/runs;
disp('Gaussian error function over a 1500x1500 matrix______: '+string(result)+' sec.');
clear('a');
clear('b');
result = 0;
for i = 1:runs
  for j=1:100
    a = rand(1000,1000,'uniform');
    b = 1:1000;
    tic();
    b = a\b';
    result = result+toc();
  end
end
result = result/runs;
disp('Linear regression over a 1000x1000 matrix____________: '+string(result)+' sec.');