--- /dev/null
+
+//###########################################################################
+// Benchmarks
+//
+//
+// Auteur : Pierre MARECHAL
+// Date : 03 Fev 2005
+//
+//###########################################################################
+
+lines(0);
+
+stacksize(30000000)
+runs = 3; // Number of times the tests are executed
+times = zeros(5,3);
+
+
+u=file('open','benchmark.ini','new');
+write(u,'[Benchmark]');
+write(u,'');
+
+disp(' Scilab Benchmark 2')
+disp('==================')
+disp('Number of times each test is run__________________________: ' + string(runs))
+disp(' ')
+
+//###########################################################################
+disp(' I. Matrix calculation')
+disp(' ---------------------')
+//###########################################################################
+
+
+//***************************************************************************
+// (1)
+//***************************************************************************
+
+cumulate = 0; a = 0; b = 0;
+for i = 1:runs
+ timer();
+ a = abs(rand(1500, 1500, 'n')/10);
+ b = a';
+ a = matrix(b, 750, 3000);
+ b = a';
+ timing = timer();
+ cumulate = cumulate + timing;
+end;
+timing = cumulate/runs;
+times(1, 1) = timing;
+
+write(u,'test1 = ' + string(timing));
+
+disp('Creation, transp., deformation of a 1500x1500 matrix (sec): ' + string(timing));
+clear('a'); clear('b');
+
+
+//***************************************************************************
+// (2)
+//***************************************************************************
+
+cumulate = 0; b = 0;
+for i = 1:runs
+ a = abs(rand(800, 800, 'n')/2);
+ timer();
+ b = a.^1000;
+ timing = timer();
+ cumulate = cumulate + timing;
+end
+timing = cumulate/runs;
+times(2, 1) = timing;
+
+write(u,'test2 = ' + string(timing));
+
+disp('800x800 normal distributed random matrix ^1000______ (sec): ' + string(timing))
+clear('a'); clear('b');
+
+
+//***************************************************************************
+// (3)
+//***************************************************************************
+
+cumulate = 0; b = 0;
+for i = 1:runs
+ a = rand(2000000, 1, 'n');
+ timer();
+ b = -sort(-a);
+ timing = timer();
+ cumulate = cumulate + timing;
+end
+timing = cumulate/runs;
+times(3, 1) = timing;
+
+write(u,'test3 = ' + string(timing));
+
+disp('Sorting of 2,000,000 random values__________________ (sec): ' + string(timing))
+clear('a'); clear('b');
+
+
+//***************************************************************************
+// (4)
+//***************************************************************************
+
+cumulate = 0; b = 0;
+for i = 1:runs
+ a = rand(700, 700, 'n');
+ timer();
+ b = a'*a;
+ timing = timer();
+ cumulate = cumulate + timing;
+end
+timing = cumulate/runs;
+times(4, 1) = timing;
+
+write(u,'test4 = ' + string(timing));
+
+disp('700x700 cross-product matrix (b = a'' * a)___________ (sec): ' + string(timing))
+clear('a'); clear('b');
+
+
+//***************************************************************************
+// (5)
+//***************************************************************************
+
+cumulate = 0; c = 0;
+for i = 1:runs
+ a = rand(600, 600, 'n');
+ b = 1:600;
+ timer();
+ c = a\b';
+ timing = timer();
+ cumulate = cumulate + timing;
+end
+timing = cumulate/runs;
+times(5, 1) = timing;
+
+write(u,'test5 = ' + string(timing));
+
+disp('Linear regression over a 600x600 matrix (c = a \ b'') (sec): ' + string(timing))
+clear('a'); clear('b'); clear('c');
+
+//***************************************************************************
+// Moyenne
+//***************************************************************************
+
+times = sort(times, 'r');
+disp(' --------------------------------------------------')
+disp(' Trimmed geom. mean (2 extremes eliminated): ' + string(exp(mean(log(times(2:4,1))))))
+disp(' ')
+
+write(u,'');
+write(u,'mean1 = ' + string(exp(mean(log(times(2:4,1))))));
+write(u,'');
+
+
+//###########################################################################
+disp(' II. Matrix functions')
+disp(' --------------------')
+//###########################################################################
+
+//***************************************************************************
+// (1)
+//***************************************************************************
+
+cumulate = 0; b = 0;
+for i = 1:runs
+ a = rand(800000, 1, 'n');
+ timer();
+ b = fft(a, -1);
+ timing = timer();
+ cumulate = cumulate + timing;
+end
+timing = cumulate/runs;
+times(1, 2) = timing;
+
+write(u,'test6 = ' + string(timing));
+
+disp('FFT over 800,000 random values______________________ (sec): ' + string(timing))
+clear('a'); clear('b');
+
+//***************************************************************************
+// (2)
+//***************************************************************************
+
+cumulate = 0; b = 0;
+for i = 1:runs
+ a = rand(320, 320, 'n');
+ timer();
+ b = spec(a);
+ timing = timer();
+ cumulate = cumulate + timing;
+end
+timing = cumulate/runs;
+times(2, 2) = timing;
+
+write(u,'test7 = ' + string(timing));
+
+disp('Eigenvalues of a 320x320 random matrix______________ (sec): ' + string(timing))
+clear('a'); clear('b');
+
+//***************************************************************************
+// (3)
+//***************************************************************************
+
+cumulate = 0; b = 0;
+for i = 1:runs
+ a = rand(650, 650, 'n');
+ timer();
+ b = det(a);
+ timing = timer();
+ cumulate = cumulate + timing;
+end
+timing = cumulate/runs;
+times(3, 2) = timing;
+
+write(u,'test8 = ' + string(timing));
+
+disp('Determinant of a 650x650 random matrix______________ (sec): ' + string(timing))
+clear('a'); clear('b');
+
+//***************************************************************************
+// (4)
+//***************************************************************************
+
+cumulate = 0; b = 0;
+for i = 1:runs
+ a = rand(900, 900, 'n');
+ a = a'*a;
+ timer();
+ b = chol(a);
+ timing = timer();
+ cumulate = cumulate + timing;
+end
+timing = cumulate/runs;
+times(4, 2) = timing;
+
+write(u,'test9 = ' + string(timing));
+
+disp('Cholesky decomposition of a 900x900 matrix__________ (sec): ' + string(timing))
+clear('a'); clear('b');
+
+//***************************************************************************
+// (5)
+//***************************************************************************
+
+cumulate = 0; b = 0;
+for i = 1:runs
+ a = rand(400, 400, 'n');
+ timer();
+ b = inv(a);
+ timing = timer();
+ cumulate = cumulate + timing;
+end
+timing = cumulate/runs;
+times(5, 2) = timing;
+
+write(u,'test10 = ' + string(timing));
+
+disp('Inverse of a 400x400 random matrix__________________ (sec): ' + string(timing))
+clear('a'); clear('b');
+
+
+//***************************************************************************
+// Moyenne
+//***************************************************************************
+
+times = sort(times, 'r');
+disp(' --------------------------------------------------')
+disp(' Trimmed geom. mean (2 extremes eliminated): ' + string(exp(mean(log(times(2:4,2))))))
+disp(' ')
+
+write(u,'');
+write(u,'mean2 = ' + string(exp(mean(log(times(2:4,2))))));
+write(u,'');
+
+
+//###########################################################################
+disp(' III. Programmation')
+disp(' ------------------')
+//###########################################################################
+
+
+//***************************************************************************
+// (1)
+//***************************************************************************
+
+cumulate = 0; a = 0; b = 0; phi = 1.6180339887498949;
+for i = 1:runs
+ a = floor(1000 * rand(750000, 1, 'u'));
+ timer();
+ b = (phi.^a - (-phi).^(-a)) / sqrt(5);
+ timing = timer();
+ cumulate = cumulate + timing;
+end
+timing = cumulate/runs;
+times(1, 3) = timing;
+
+write(u,'test11 = ' + string(timing));
+
+disp('750,000 Fibonacci numbers calculation (vector calc)_ (sec): ' + string(timing))
+clear('a'); clear('b'); clear('phi');
+
+
+//***************************************************************************
+// (2)
+//***************************************************************************
+
+cumulate = 0; a = 2250; b = 0;
+for i = 1:runs
+ timer();
+ b = ones(a, a)./((1:a)' * ones(1, a) + ones(a, 1) * (0:(a-1)));
+ timing = timer();
+ cumulate = cumulate + timing;
+end
+timing = cumulate/runs;
+times(2, 3) = timing;
+
+write(u,'test12 = ' + string(timing));
+
+disp('Creation of a 2250x2250 Hilbert matrix (matrix calc) (sec): ' + string(timing))
+clear('a'); clear('b');
+
+
+//***************************************************************************
+// (3)
+//***************************************************************************
+
+cumulate = 0; c = 0;
+deff('[c]=gcd2(a, b)',['if b <= 1.0E-4'; 'c = a;'; 'else'; 'b(b == 0) = a(b == 0);'; 'c = gcd2(b, modulo(a, b));'; 'end'])
+for i = 1:runs
+ a = ceil(1000 * rand(70000, 1, 'u'));
+ b = ceil(1000 * rand(70000, 1, 'u'));
+ timer();
+ c = gcd2(a, b); // gcd2 is a recursive function
+ timing = timer();
+ cumulate = cumulate + timing;
+end
+timing = cumulate/runs;
+times(3, 3) = timing;
+
+write(u,'test13 = ' + string(timing));
+
+disp('Grand common divisors of 70,000 pairs (recursion)___ (sec): ' + string(timing))
+clear('a'); clear('b'); clear('c');
+
+
+//***************************************************************************
+// (4)
+//***************************************************************************
+
+cumulate = 0; b = 0;
+for i = 1:runs
+ b = zeros(220, 220);
+ timer();
+ for j = 1:220
+ for k = 1:220
+ b(k,j) = abs(j - k) + 1;
+ end
+ end
+ timing = timer();
+ cumulate = cumulate + timing;
+end
+timing = cumulate/runs;
+times(4, 3) = timing;
+
+write(u,'test14 = ' + string(timing));
+
+disp('Creation of a 220x220 Toeplitz matrix (loops)_______ (sec): ' + string(timing))
+clear('b'); clear('j'); clear('k');
+
+
+//***************************************************************************
+// (5)
+//***************************************************************************
+
+cumulate = 0; p = 0; q = 0; vt = 0; vr = 0; vrt = 0; rvt = 0; RV = 0; j = 0; k = 0;
+x2 = 0; R = 0; Rxx = 0; Ryy = 0; Rxy = 0; Ryx = 0; Rvmax = 0; f = 0;, c= 0; d = 0;
+for i = 1:runs
+ x = abs(rand(37, 37, 'n'));
+ timer();
+ // Calculation of Escoufier's equivalent vectors
+ p = size(x, 2);
+ vt = [1:p]; // Variables to test
+ vr = []; // Result: ordered variables
+ RV = [1:p]; // Result: correlations
+ for j = 1:p // loop on the variable number
+ Rvmax = 0;
+ for k = 1:(p-j+1) // loop on the variables
+ x2 = [x, x(:, vr), x(:, vt(k))]; // New table to test
+ // R = corrcoef(x2); // Correlations table
+ // Not in scilab, so the 5 following lines do it
+ q = size(x2, 1);
+ x2 = x2 - ones (q, 1) * sum (x2, 'r') / q;
+ c = conj(x2' * x2 / (q - 1));
+ d=diag(c);
+ R = c ./ sqrt(d * d');
+ Ryy = R(1:p, 1:p);
+ Rxx = R(p+1:p+j, p+1:p+j);
+ Rxy = R(p+1:p+j, 1:p);
+ Ryx = Rxy';
+ rvt = trace(Ryx*Rxy)/((trace(Ryy^2)*trace(Rxx^2))^0.5); // RV calculation
+ if rvt > Rvmax
+ Rvmax = rvt; // test of RV
+ vrt(j) = vt(k); // temporary held variable
+ end
+ end
+ vr(j) = vrt(j); // Result: variable
+ RV(j) = Rvmax; // Result: correlation
+ f = find(vt~=vr(j)); // identify the held variable
+ vt = vt(f); // reidentify variables to test
+ end
+ timing = timer();
+ cumulate = cumulate + timing;
+end
+times(5, 3) = timing;
+
+write(u,'test15 = ' + string(timing));
+
+disp('Escoufier''s method on a 37x37 matrix (mixed)________ (sec): ' + string(timing))
+clear('x'); clear('p'); clear('q'); clear('vt'); clear('vr'); clear('vrt'); clear('rvt'); clear('RV'); clear('j'); clear('k');
+clear('x2'); clear('R'); clear('Rxx'); clear('Ryy'); clear('Rxy'); clear('Ryx'); clear('Rvmax'); clear('f'); clear('c'); clear('d');
+
+
+//***************************************************************************
+// Moyenne
+//***************************************************************************
+
+times = sort(times, 'r');
+disp(' --------------------------------------------------')
+disp(' Trimmed geom. mean (2 extremes eliminated): ' + string(exp(mean(log(times(2:4,3))))))
+disp(' ')
+
+write(u,'');
+write(u,'mean3 = ' + string(exp(mean(log(times(2:4,3))))));
+write(u,'');
+
+//###########################################################################
+//###########################################################################
+
+disp(' ')
+disp('Total time for all 15 tests_________________________ (sec): ' + string(sum(times)))
+disp('Overall mean (sum of I, II and III trimmed means/3)_ (sec): ' + string(exp(mean(log(times(2:4,:))))))
+disp(' --- End of test ---')
+
+write(u,'');
+write(u,'total = ' + string(sum(times)));
+write(u,'mean4 = ' + string(exp(mean(log(times(2:4,:))))));
+file('close',u);
+
+clear('cumulate'); clear('timing'); clear('times'); clear('runs'); clear('i');
+
+exit