An introduction to the GNU Scientific Library (GSL)

The GNU Scientific Library (GSL) is a numerical library for C and C++. It is free software. It provides lots of data structures and methods for mathematics, including vectors and matrices, BLAS routines, advanced operations of linear algebra and optimization, and much more.

Personally, I mainly use GSL to work with matrices and BLAS routines. With GSL, you can play with matrices and vectors in a high-level object-oriented way, a bit like in Matlab or Python+Numpy. No more headaches with double pointers and hand-coded indexing!

In this post, I assume you already know how to code in C.

Install

To install it on a GNU/Linux distribution with apt, just type in a terminal apt install libgsl23 libgsl-dev (for other systems, see the project home page).

Example

Let’s begin with a simple example of C code using GSL. Create a file called test.c, and copy-paste the following code. It creates a random matrix, and then normalizes its columns:

/* To compile: gcc test.c -lgsl -lgslcblas -o test */

#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>

/* Function to print a gsl_matrix in a human-readable way */
void
print_gslmat (const gsl_matrix * m)
{
  for (uint i = 0; i < m->size1; ++i)
    {
      for (uint j = 0; j < m->size2; ++j)
	    printf ("%g ", gsl_matrix_get (m, i, j));
      printf ("\n");
    }
}

int
main ()
{
  /* Create a matrix with 3 rows and 4 columns */
  gsl_matrix * mat = gsl_matrix_alloc (3, 4);

  /* Fill it with random numbers between 0 and 1 */
  for (uint i = 0; i < mat->size1; ++i)
    for (uint j = 0; j < mat->size2; ++j)
      gsl_matrix_set (mat, i, j, (double)rand()/RAND_MAX);

  /* Display the matrix */
  printf ("Original matrix\n");
  print_gslmat (mat);

  for (uint j = 0; j < mat->size2; ++j)
    {
      /* Create a view of the j-th column of mat */
      gsl_vector_view col = gsl_matrix_column (mat, j);
      /* Compute the l1-norm of this column (sum of entries) */
      double d = gsl_blas_dasum (&col.vector);
      /* Divide all entries of this column by its norm */
      gsl_blas_dscal (1 / d, &col.vector);
    }

  /* Display the matrix */
  printf ("Normalized matrix\n");
  print_gslmat (mat);

  /* Exit */
  gsl_matrix_free (mat);
  return 0;
}

You can compile it with gcc test.c -lgsl -lgslcblas -o test, and then run it with ./test. If it works, then GSL is well installed on your computer! Let’s now examine this code.

Details

First, you need to include the library headers. There is one header for every area covered by the GSL. In our case, we juste need the headers matrix and blas.

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>

To create a matrix, you have to use alloc functions similar to the standard C malloc, taking in parameters the dimensions (number of rows and number of columns) of the matrix.

gsl_matrix * mat = gsl_matrix_alloc (3, 4);

Note that alloc only allocates the memory, and does not initialize the matrix. It can contain garbage values! Using gsl_matrix_calloc instead, the matrix is initialized with all entries to zero. For vectors, the corresponding functions are gsl_vector_alloc (n) and gsl_vector_calloc (n).

Note also that the size of the matrix is constant. You cannot add or remove rows or lines without creating a new matrix. It can seem inconvenient, but it actually helps making efficient code, by giving you explicit control of the memory. In languages such as Matlab, you may think you can easily change the size of a matrix, but it actually ends up in very costly hidden operations.

Of course, when you are done using a matrix, you have to free the allocated memory.

gsl_matrix_free (mat);

You can get the dimensions of the matrix mat with mat->size1 and mat->size2. You can set a value to an entry with gsl_matrix_set, and read it with gsl_matrix_get.

for (uint i = 0; i < mat->size1; ++i)
  for (uint j = 0; j < mat->size2; ++j)
    gsl_matrix_set (mat, i, j, (double)rand()/RAND_MAX);

When passing a GSL matrix as argument of a function, you only pass the pointer to this matrix. It makes code efficient, but be careful to not modify a matrix by mistake! The following function prints a matrix in a simple and readable way.

void
print_gslmat (const gsl_matrix * m)
{
  for (uint i = 0; i < m->size1; ++i)
    {
      for (uint j = 0; j < m->size2; ++j)
	    printf ("%g ", gsl_matrix_get (m, i, j));
      printf ("\n");
    }
}

You can create views of vectors and matrices. They allow you to work on only a part of a vector or matrix, as if they were independent objects.

/* Create a view of the j-th column of mat */
gsl_vector_view col = gsl_matrix_column (mat, j);

Using BLAS operations makes your code easier to write and read. It also makes it faster, because it uses low-level optimized code. The functions’ naming can seem a bit strange, but it actually follows the BLAS standard.

/* Compute the l1-norm of this column (sum of entries) */
double d = gsl_blas_dasum (&col.vector);
/* Divide all entries of this column by its norm */
gsl_blas_dscal (1 / d, &col.vector);

Before coding an operation by hand, always check if a built-in function exists!

Remarks

As opposed to Matlab for example, matrices in GSL are stored in a row-major order. This means that entries of a given row are contiguous in memory, so it is significantly faster to traverse a matrix row after row. Think about it when coding!

The nice thing about BLAS and its GSL wrapper, is that it is a standard interface independent of the implementation, so you can code using high-level functions no matter what is going on under the hood. There are many libraries implementing BLAS, and I personally use ATLAS. It’s quite faster than the default library CBLAS. To use it in your projects, just install the library and compile with the -latlas option.

To see all of the possibilities offered by GSL, have a look at the official documentation.