// Ox Benchmark 2 (8 March 2003) // version 2, scaled to get 1 +/- 0.1 sec with R 1.6.2 // using the standard ATLAS library (Rblas.dll) // on a Pentium IV 1.6 Ghz with 1 Gb Ram on Win XP pro // Author : Philippe Grosjean // eMail : phgrosjean@sciviews.org // Web : http://www.sciviews.org // License: GPL 2 or above at your convenience (see: http://www.gnu.org) // // Several tests are adapted from: //*************************************************************************** //* Matlab Benchmark program version 2.0 * //* Author : Stefan Steinhaus * //* EMAIL : stst@informatik.uni-frankfurt.de * //* This program is public domain. Feel free to copy it freely. * //*************************************************************************** // and from the corresponding Ox Benchmark V.2 from S. Steinhaus // Escoufier's equivalents vectors (III.5) is adapted from Planque & Fromentin, 1996 // Ref: Escoufier Y., 1970. Echantillonnage dans une population de variables // aleatoires réelles. Publ. Inst. Statis. Univ. Paris 19 Fasc 4, 1-47. // // On the dos prompt, type: // cd c:\ // oxl ox2 // to start the test #include // This is for test III.3 gcd2(a, b) { if (b <= 1.0E-4) return a; else { b[vecindex(b, 0)] = a[vecindex(b, 0)]; return gcd2(b, fmod(a, b)); } } main() { decl runs = 3; // Number of times the tests are executed decl times, cumul, timeini; decl a, b, c, phi, i, j, k; decl x, x2, p, vt, vr, vrt, rvt; decl RV, R, Rxx, Ryy, Rxy, Ryx, Rvmax; times = zeros(5, 3); print(" Ox Benchmark 2\n"); print(" ==============\n"); print("Number of times each test is run__________________________: ", runs, "\n\n"); print(" I. Matrix calculation\n"); print(" ---------------------\n"); // (1) cumul = 0; a = 0; b = 0; for (i = 0; i < runs; i++) { timeini = timer(); a = fabs(rann(1500, 1500)/10); b = a'; a = reshape(b, 750, 3000); b = a'; cumul += (timer() - timeini) / 100; } cumul /= runs; times[0][0] = cumul; print("Creation, transp., deformation of a 1500x1500 matrix (sec): " , cumul, "\n"); delete a; delete b; // (2) cumul = 0; b = 0; for (i = 0; i < runs; i++) { a = fabs(rann(800, 800)/2); timeini = timer(); b = (a).^1000; cumul += (timer() - timeini) / 100; } cumul /= runs; times[1][0] = cumul; print("800x800 normal distributed random matrix ^1000______ (sec): ", cumul, "\n"); delete a; delete b; // (3) cumul = 0; b = 0; for (i = 0; i < runs; i++) { a = rann(2000000, 1); timeini = timer(); b = sortc(a); cumul += (timer() - timeini) / 100; } cumul /= runs; times[2][0] = cumul; print("Sorting of 2,000,000 random values__________________ (sec): ", cumul, "\n"); delete a; delete b; //(4) cumul = 0; b = 0; for (i = 0; i < runs; i++) { a = rann(700, 700); timeini = timer(); b = a'a; cumul += (timer() - timeini) / 100; } cumul /= runs; times[3][0] = cumul; print("700x700 cross-product matrix (b = a' * a)___________ (sec): ", cumul, "\n"); delete a; delete b; // (5) cumul = 0; c = 0; for (i = 0; i < runs; i++) { a = rann(600, 600); b = range(1, 600); timeini = timer(); olsr(b, a, &c); cumul += (timer() - timeini) / 100; } cumul /= runs; times[4][0] = cumul; print("Linear regression over a 600x600 matrix (c = a \ b') (sec): ", cumul, "\n"); delete a; delete b; delete c; times = sortc(times); print(" --------------------------------------------\n"); print(" Trimmed geom. mean (2 extremes eliminated): ", exp(meanc(log(times[1:3][0]))), "\n\n"); print(" II. Matrix functions\n"); print(" --------------------\n"); // (1) cumul = 0; b = 0; for (i = 0; i < runs; i++) { a = rann(1, 800000); timeini = timer(); b = fft(a); cumul += (timer() - timeini) / 100; } cumul /= runs; times[0][1] = cumul; print("FFT over 800,000 random values______________________ (sec): ", cumul, "\n"); delete a; delete b; // (2) cumul = 0; b = 0; for (i = 0; i < runs; i++) { a = rann(320, 320); timeini = timer(); eigen(a, &b); cumul += (timer() - timeini) / 100; } cumul /= runs; times[1][1] = cumul; print("Eigenvalues of a 320x320 random matrix______________ (sec): ", cumul, "\n"); delete a; delete b; // (3) cumul = 0; b = 0; for (i = 0; i < runs; i++) { a = rann(650, 650)/10; timeini = timer(); b = determinant(a); cumul += (timer() - timeini) / 100; } cumul /= runs; times[2][1] = cumul; print("Determinant of a 650x650 random matrix______________ (sec): ", cumul, "\n"); delete a; delete b; // (4) cumul = 0; b = 0; for (i = 0; i < runs; i++) { a = rann(900, 900); a = a'a; timeini = timer(); b = choleski(a); cumul += (timer() - timeini) / 100; } cumul /= runs; times[3][1] = cumul; print("Cholesky decomposition of a 900x900 matrix__________ (sec): ", cumul, "\n"); delete a; delete b; // (5) cumul = 0; b = 0; for (i = 0; i < runs; i++) { a = rann(400, 400); timeini = timer(); b = invert(a); cumul += (timer() - timeini) / 100; } cumul /= runs; times[4][1] = cumul; print("Inverse of a 400x400 random matrix__________________ (sec): ", cumul, "\n"); delete a; delete b; times = sortc(times); print(" --------------------------------------------\n"); print(" Trimmed geom. mean (2 extremes eliminated): ", exp(meanc(log(times[1:3][1]))), "\n\n"); print(" III. Programmation\n"); print(" ------------------\n"); // (1) cumul = 0; a = 0; b = 0; phi = 1.6180339887498949; for (i = 0; i < runs; i++) { a = floor(1000 * ranu(750000, 1)); timeini = timer(); b = (phi^a - (-phi)^(-a)) / sqrt(5); cumul += (timer() - timeini) / 100; } cumul /= runs; times[0][2] = cumul; print("750,000 Fibonacci numbers calculation (vector calc)_ (sec): ", cumul, "\n"); delete a; delete b; delete phi; // (2) cumul = 0; a = 2250; b = 0; for (i = 0; i < runs; i++) { timeini = timer(); b = ones(a, a)./((range(1, a))' * ones(1, a) + ones(a, 1) * (range(0, a-1))); cumul += (timer() - timeini) / 100; } cumul /= runs; times[1][2] = cumul; print("Creation of a 2250x2250 Hilbert matrix (matrix calc) (sec): ", cumul, "\n"); delete a; delete b; // (3) cumul = 0; c = 0; for (i = 0; i < runs; i++) { a = ceil(1000 * ranu(70000, 1)); b = ceil(1000 * ranu(70000, 1)); timeini = timer(); c = gcd2(a, b); // gcd2 is a recursive function cumul += (timer() - timeini) / 100; } cumul /= runs; times[2][2] = cumul; print("Grand common divisors of 70,000 pairs (recursion)___ (sec): ", cumul, "\n"); delete a; delete b; delete c; // (4) cumul = 0; b = 0; for (i = 0; i < runs; i++) { b = zeros(220, 220); timeini = timer(); for (j = 0; j < 220; j++) { for (k = 0; k < 220; k++) { b[k][j] = fabs(j - k) + 1; } } cumul += (timer() - timeini) / 100; } cumul /= runs; times[3][2] = cumul; print("Creation of a 220x220 Toeplitz matrix (loops)_______ (sec): ", cumul, "\n"); delete b; delete j; delete k; // (5) cumul = 0; p = 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; for (i = 0; i < runs; i++) { x = fabs(rann(37, 37)); timeini = timer(); // Calculation of Escoufier's equivalent vectors p = sizec(x) - 1; vt = range(0, p); // Variables to test vr = ones(p+1, 1); // Result: ordered variables RV = ones(p+1, 1); // Result: correlations for (j = 0; j < p+1; j++) { // loop on the variable number Rvmax = 0; for (k = 0; k < p-j+1; k++) { // loop on the variables x2 = x~x[][vt[k]]; // New table to test R = correlation(x2); // Correlations table Ryy = R[0:p][0:p]; Rxx = R[p+1:p+j+1][p+1:p+j+1]; Rxy = R[p+1:p+j+1][0: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 = vt[k]; // temporary held variable } } vr[j][0] = vrt + 1; // +1 because we index variables from 1! RV[j][0] = Rvmax; // Result: correlation x = x~x[][vrt]; // add the held variable to x vt = deletec(vt, vrt); // reidentify variables to test } cumul += (timer() - timeini) / 100; } cumul /= runs; times[4][2] = cumul; print("Escoufier's method on a 37x37 matrix (mixed)________ (sec): ", cumul, "\n"); delete x; delete p; delete vt; delete vr; delete vrt; delete rvt; delete RV; delete j; delete k; delete x2; delete R; delete Rxx; delete Ryy; delete Rxy; delete Ryx; delete Rvmax; times = sortc(times); print(" --------------------------------------------\n"); print(" Trimmed geom. mean (2 extremes eliminated): ", exp(meanc(log(times[1:3][2]))), "\n\n"); print("\n"); print("Total time for all 15 tests_________________________ (sec): ", sumr(sumc(times)), "\n"); print("Overall mean (sum of I, II and III trimmed means/3)_ (sec): ", exp(meanr(meanc(log(times[1:3][0:2])))), "\n"); delete cumul; delete timeini; delete times; delete runs; delete i; print(" --- End of test ---\n"); }