I was trying to learn how to use gsll (the lisp bindings of GNU Scientific Library), specifically calculating eigenvalues and eigenvectors in gsll. I wrote the following code to calculate the eigenvalues and eigenvectors of a 400x400 randomly generated real matrix.
(let* ((N 400)
(rng (gsll:make-random-number-generator gsll::+mt19937+))
(aa (grid:make-foreign-array
'double-float
:initial-contents (loop repeat N collect
(loop repeat N collect
(gsll:sample rng :uniform))))))
(gsll:eigenvalues-eigenvectors-nonsymm aa))
I ran this code 10 times and the average came out to be 1.98 seconds on SBCL on my debian box (Intel Core2 2.8GHz with 2Gb RAM) at work.
Then I thought, the code and how fast it runs does not mean anything. To put this in a context, I needed to write the code for the same problem in C and see how fast it runs. So, I wrote the following code in C:
#include <malloc.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_rng.h>
void main(void)
{
double *data;
int i, j, n=400;
const gsl_rng_type *T;
gsl_rng *r;
gsl_vector_complex *eval;
gsl_matrix_complex *evec;
gsl_eigen_nonsymmv_workspace *w;
gsl_matrix_view m;
gsl_rng_env_setup();
T = gsl_rng_default;
r = gsl_rng_alloc (T);
data = (double *)malloc(n*n*sizeof(double));
for (i = 0; i < n; i++)
for (j = 0; j < n; j++)
data[n*i+j] = gsl_rng_uniform (r);
gsl_rng_free(r);
m = gsl_matrix_view_array(data, n, n);
eval = gsl_vector_complex_alloc(n);
evec = gsl_matrix_complex_alloc(n,n);
w = gsl_eigen_nonsymmv_alloc(n);
gsl_eigen_nonsymmv(&m.matrix, eval, evec, w);
gsl_eigen_nonsymmv_free(w);
gsl_vector_complex_free(eval);
gsl_matrix_complex_free(evec);
free(data);
}
I ran this 10 times. Guess what: the average was 1.88 seconds. I wasn’t expecting a close result. Then I did an experiment with different values of N.
The results are averages of 3 runs per each N. The C code is modified so that I provided N on the command line. Lisp A results come from running the code within the SBCL repl, and Lisp B results are obtained by running SBCL on the command line with the specific N provided.
I was surprised by the fact that virtually there is no difference between the C results and the lisp results obtained within the repl. Notice also the notch in Lisp A between N=50 to N=100. My guess is that the code is internally compiled once for the Lisp A runs. So, Lisp A and Lisp B results are identical for N=50, and there is a divergence for later N because SBCL compiles the code internally each time it runs from the command line.