# R 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 # revised and optimized for R v. 1.5.x, 8 June 2002 # Requires additionnal libraries: Matrix, SuppDists # 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 the Splus Benchmark Test V. 2 # by Stephan Steinhaus (stst@informatik.uni-frankfurt.de) # # Reference for Escoufier's equivalents vectors (test III.5): # Escoufier Y., 1970. Echantillonnage dans une population de variables # aleatoires réelles. Publ. Inst. Statis. Univ. Paris 19 Fasc 4, 1-47. # # type source("c://R2.R") to start the test runs <- 3 # Number of times the tests are executed times <- rep(0, 15); dim(times) <- c(5,3) require(Matrix) # Optimized matrix operations require(SuppDists) # Optimized random number generators Runif <- rMWC1019 # The fast uniform number generator # If you don't have SuppDists, you can use: Runif <- runif a <- rMWC1019(10, new.start=TRUE, seed=492166) # Init. the generator Rnorm <- rziggurat # The fast normal number generator # If you don't have SuppDists, you can use: Rnorm <- rnorm b <- rziggurat(10, new.start=TRUE) # Init. the generator remove("a", "b") options(object.size=100000000) cat("\n\n R Benchmark 2\n") cat(" =============\n") cat(c("Number of times each test is run__________________________: ", runs)) cat("\n\n") cat(" I. Matrix calculation\n") cat(" ---------------------\n") if (R.Version()$os == "Win32" || R.Version()$os == "mingw32") flush.console() # (1) cumulate <- 0; a <- 0; b <- 0 for (i in 1:runs) { timing <- system.time({ a <- abs(matrix(Rnorm(1500*1500)/10, ncol=1500, nrow=1500)); b <- t(a); dim(b) <- c(750, 3000); a <- t(b) })[3] cumulate <- cumulate + timing } timing <- cumulate/runs times[1, 1] <- timing cat(c("Creation, transp., deformation of a 1500x1500 matrix (sec): ", timing, "\n")) remove("a", "b") if (R.Version()$os == "Win32" || R.Version()$os == "mingw32") flush.console() # (2) cumulate <- 0; b <- 0 for (i in 1:runs) { a <- abs(matrix(Rnorm(800*800)/2, ncol=800, nrow=800)); timing <- system.time({ b <- a^1000 })[3] cumulate <- cumulate + timing } timing <- cumulate/runs times[2, 1] <- timing cat(c("800x800 normal distributed random matrix ^1000______ (sec): ", timing, "\n")) remove("a", "b") if (R.Version()$os == "Win32" || R.Version()$os == "mingw32") flush.console() # (3) cumulate <- 0; b <- 0 for (i in 1:runs) { a <- Rnorm(2000000) timing <- system.time({ b <- sort(a, method="quick") # Sort is modified in v. 1.5.x # And there is now a quick method that better competes with other packages!!! })[3] cumulate <- cumulate + timing } timing <- cumulate/runs times[3, 1] <- timing cat(c("Sorting of 2,000,000 random values__________________ (sec): ", timing, "\n")) remove("a", "b") if (R.Version()$os == "Win32" || R.Version()$os == "mingw32") flush.console() # (4) cumulate <- 0; b <- 0 for (i in 1:runs) { a <- Rnorm(700*700); dim(a) <- c(700, 700) timing <- system.time({ b <- crossprod(a) # equivalent to: b <- t(a) %*% a })[3] cumulate <- cumulate + timing } timing <- cumulate/runs times[4, 1] <- timing cat(c("700x700 cross-product matrix (b = a' * a)___________ (sec): ", timing, "\n")) remove("a", "b") if (R.Version()$os == "Win32" || R.Version()$os == "mingw32") flush.console() # (5) cumulate <- 0; c <- 0; qra <-0 for (i in 1:runs) { a <- Rnorm(600*600); dim(a) <- c(600,600) b <- 1:600 timing <- system.time({ qra <- qr(a, tol = 1e-7); c <- qr.coef(qra, b) #Rem: a little faster than c <- lsfit(a, b, inter=F)$coefficients })[3] cumulate <- cumulate + timing } timing <- cumulate/runs times[5, 1] <- timing cat(c("Linear regression over a 600x600 matrix (c = a \\ b') (sec): ", timing, "\n")) remove("a", "b", "c", "qra") if (R.Version()$os == "Win32" || R.Version()$os == "mingw32") flush.console() times[ , 1] <- sort(times[ , 1]) cat(" --------------------------------------------\n") cat(c(" Trimmed geom. mean (2 extremes eliminated): ", exp(mean(log(times[2:4, 1]))), "\n\n")) cat(" II. Matrix functions\n") cat(" --------------------\n") if (R.Version()$os == "Win32") flush.console() # (1) cumulate <- 0; b <- 0 for (i in 1:runs) { a <- Rnorm(800000) timing <- system.time({ b <- fft(a) })[3] cumulate <- cumulate + timing } timing <- cumulate/runs times[1, 2] <- timing cat(c("FFT over 800,000 random values______________________ (sec): ", timing, "\n")) remove("a", "b") if (R.Version()$os == "Win32" || R.Version()$os == "mingw32") flush.console() # (2) cumulate <- 0; b <- 0 for (i in 1:runs) { a <- Rnorm(320*320); dim(a) <- c(320, 320) # Only needed if using eigen.Matrix(): Matrix.class(a) timing <- system.time({ b <- eigen.default(a, symmetric=F, only.values=T)$Value # Rem: on my machine, it is faster than: # b <- La.eigen(a, symmetric=F, only.values=T, method="dsyevr")$Value # b <- La.eigen(a, symmetric=F, only.values=T, method="dsyev")$Value # b <- eigen.Matrix(a, vectors = F)$Value })[3] cumulate <- cumulate + timing } timing <- cumulate/runs times[2, 2] <- timing cat(c("Eigenvalues of a 320x320 random matrix______________ (sec): ", timing, "\n")) remove("a", "b") if (R.Version()$os == "Win32" || R.Version()$os == "mingw32") flush.console() # (3) cumulate <- 0; b <- 0 for (i in 1:runs) { a <- Rnorm(650*650); dim(a) <- c(650, 650) Matrix.class(a) timing <- system.time({ b <- det.Matrix(a, logarithm=F) # Rem: the following is slower on my computer! # b <- det.default(a) })[3] cumulate <- cumulate + timing } timing <- cumulate/runs times[3, 2] <- timing cat(c("Determinant of a 650x650 random matrix______________ (sec): ", timing, "\n")) remove("a", "b") if (R.Version()$os == "Win32" || R.Version()$os == "mingw32") flush.console() # (4) cumulate <- 0; b <- 0 for (i in 1:runs) { a <- Rnorm(900*900); dim(a) <- c(900, 900) a <- crossprod(a, a) timing <- system.time({ b <- chol(a) })[3] cumulate <- cumulate + timing } timing <- cumulate/runs times[4, 2] <- timing cat(c("Cholesky decomposition of a 900x900 matrix__________ (sec): ", timing, "\n")) remove("a", "b") if (R.Version()$os == "Win32" || R.Version()$os == "mingw32") flush.console() # (5) cumulate <- 0; b <- 0 for (i in 1:runs) { a <- Rnorm(400*400); dim(a) <- c(400, 400) timing <- system.time({ # b <- qr.solve(a) # Rem: a little faster than b <- solve(a) })[3] cumulate <- cumulate + timing } timing <- cumulate/runs times[5, 2] <- timing cat(c("Inverse of a 400x400 random matrix__________________ (sec): ", timing, "\n")) remove("a", "b") if (R.Version()$os == "Win32" || R.Version()$os == "mingw32") flush.console() times[ , 2] <- sort(times[ , 2]) cat(" --------------------------------------------\n") cat(c(" Trimmed geom. mean (2 extremes eliminated): ", exp(mean(log(times[2:4, 2]))), "\n\n")) cat(" III. Programmation\n") cat(" ------------------\n") if (R.Version()$os == "Win32" || R.Version()$os == "mingw32") flush.console() # (1) cumulate <- 0; a <- 0; b <- 0; phi <- 1.6180339887498949 for (i in 1:runs) { a <- floor(Runif(750000)*1000) timing <- system.time({ b <- (phi^a - (-phi)^(-a))/sqrt(5) })[3] cumulate <- cumulate + timing } timing <- cumulate/runs times[1, 3] <- timing cat(c("750,000 Fibonacci numbers calculation (vector calc)_ (sec): ", timing, "\n")) remove("a", "b", "phi") if (R.Version()$os == "Win32" || R.Version()$os == "mingw32") flush.console() # (2) cumulate <- 0; a <- 2250; b <- 0 for (i in 1:runs) { timing <- system.time({ b <- rep(1:a, a); dim(b) <- c(a, a); b <- 1 / (t(b) + 0:(a-1)) # Rem: this is twice as fast as the following code proposed by R programmers # a <- 1:a; b <- 1 / outer(a - 1, a, "+") })[3] cumulate <- cumulate + timing } timing <- cumulate/runs times[2, 3] <- timing cat(c("Creation of a 2250x2250 Hilbert matrix (matrix calc) (sec): ", timing, "\n")) remove("a", "b") if (R.Version()$os == "Win32" || R.Version()$os == "mingw32") flush.console() # (3) cumulate <- 0; c <- 0 gcd2 <- function(x, y) {if (sum(y > 1.0E-4) == 0) x else {y[y == 0] <- x[y == 0]; Recall(y, x %% y)}} for (i in 1:runs) { a <- ceiling(Runif(70000)*1000) b <- ceiling(Runif(70000)*1000) timing <- system.time({ c <- gcd2(a, b) # gcd2 is a recursive function })[3] cumulate <- cumulate + timing } timing <- cumulate/runs times[3, 3] <- timing cat(c("Grand common divisors of 70,000 pairs (recursion)___ (sec): ", timing, "\n")) remove("a", "b", "c", "gcd2") if (R.Version()$os == "Win32" || R.Version()$os == "mingw32") flush.console() # (4) cumulate <- 0; b <- 0 for (i in 1:runs) { b <- rep(0, 220*220); dim(b) <- c(220, 220) timing <- system.time({ # Rem: there are faster ways to do this # but here we want to time loops (220*220 'for' loops)! for (j in 1:220) { for (k in 1:220) { b[k,j] <- abs(j - k) + 1 } } })[3] cumulate <- cumulate + timing } timing <- cumulate/runs times[4, 3] <- timing cat(c("Creation of a 220x220 Toeplitz matrix (loops)_______ (sec): ", timing, "\n")) remove("b", "j", "k") if (R.Version()$os == "Win32" || R.Version()$os == "mingw32") flush.console() # (5) cumulate <- 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 # Calculate the trace of a matrix (sum of its diagonal elements) Trace <- function(y) {sum(c(y)[1 + 0:(min(dim(y)) - 1) * (dim(y)[1] + 1)], na.rm=FALSE)} for (i in 1:runs) { x <- abs(Rnorm(37*37)); dim(x) <- c(37, 37) timing <- system.time({ # Calculation of Escoufier's equivalent vectors p <- ncol(x) vt <- 1:p # Variables to test vr <- NULL # Result: ordered variables RV <- 1:p # Result: correlations vrt <- NULL for (j in 1:p) { # loop on the variable number Rvmax <- 0 for (k in 1:(p-j+1)) { # loop on the variables x2 <- cbind(x, x[,vr], x[,vt[k]]) R <- cor(x2) # Correlations table 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 <- t(Rxy) rvt <- Trace(Ryx %*% Rxy) / sqrt(Trace(Ryy %*% Ryy) * Trace(Rxx %*% Rxx)) # RV calculation if (rvt > Rvmax) { Rvmax <- rvt # test of RV vrt <- vt[k] # temporary held variable } } vr[j] <- vrt # Result: variable RV[j] <- Rvmax # Result: correlation vt <- vt[vt!=vr[j]] # reidentify variables to test } })[3] cumulate <- cumulate + timing } times[5, 3] <- timing cat(c("Escoufier's method on a 37x37 matrix (mixed)________ (sec): ", timing, "\n")) remove("x", "p", "vt", "vr", "vrt", "rvt", "RV", "j", "k") remove("x2", "R", "Rxx", "Ryy", "Rxy", "Ryx", "Rvmax", "Trace") if (R.Version()$os == "Win32" || R.Version()$os == "mingw32") flush.console() times[ , 3] <- sort(times[ , 3]) cat(" --------------------------------------------\n") cat(c(" Trimmed geom. mean (2 extremes eliminated): ", exp(mean(log(times[2:4, 3]))), "\n\n\n")) cat(c("Total time for all 15 tests_________________________ (sec): ", sum(times), "\n")) cat(c("Overall mean (sum of I, II and III trimmed means/3)_ (sec): ", exp(mean(log(times[2:4, ]))), "\n")) remove("cumulate", "timing", "times", "runs", "i") cat(" --- End of test ---\n\n")