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.