Eigenvalues and Eigenvectors in GSLL

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.