Saturday, June 13, 2009

GSL: GNU Scientific Library

I was looking into the GSL for a project. I am puzzled about the design of the LU solver. To solve a system of linear equations Ax=b, one would call gsl_linalg_LU_decomp followed by gsl_linalg_LU_solve. The strange thing is: non of these functions has a return code for “sorry: the matrix was singular”. When passing a singular matrix both routines will return GSL_SUCCESS.

I thought it would be a good idea to give the GSL people some feedback on this:

Looking at the linear algebra code, I have to say the implementation looks rather unsophisticated compared to say LINPACK or LAPACK.

Update: the GSL people answered here. Better than nothing, but not optimal. You really want this flag in the decomposition as it is detected there, and because it makes it easier to write:

call decomp(A);
if not singular then
   for i := 1 to NumRHS do
      call solve(rhs[i]);

Never thought I would say this, but it looks like LINPACK/LAPACK is more user-friendly than GSL in this respect.