GSL API design flaws

["E. Robert Tisdale" <edwin at netwood dot net>]

There appear to be several fundamental design flaws
in the GNU Scientific Library (GSL)
Application Programmers' Interface (API).

The GSL API designers need to answer a few basic questions:

  1. Who is supposed to use the GSL?
  2. What are the priorities?
  3. What Abstract Data Type (ADT) does the GSL support?
  4. How should the API conceal the details
     of the representation and implementation?

1. Who uses the GSL?
   Numerical application programmers?
   Or numerical library developers?
   Numerical application programmers write function main()
   and other application specific functions.
   The Basic Linear Algebra Subroutine (BLAS) library,
   on the other hand, was designed to help implement
   numerical libraries like the Linear Algebra PACKage (LAPACK)
   and not to help implement numerical application programs.

2. What are the priorities?
   0. performance,
   1. safety and
   2. convenience.

   Error handling implies error checking
   which can have a very negative impact on performance.
   If your first concern is performance,
   it is best to leave error checking and handling
   up to the application programmer.

   The GSL API should avoid specifying any function
   which might generate an error or exception.
   It is better, for example, to provide matrix decompositions
   and solvers instead of a matrix inverse
   so that the numerical application programmer
   can compute a condition number from the decomposition
   before applying the solver.

3. What ADT does the GSL support?
   The GSL supports something a little more specific
   than abstract vector and matrix objects.
   It appears that the GSL supports

     a one dimensional array of numbers referenced through
     strided real and complex multidimensional views

   except that gsl_matrix views are restricted
   to unit stride between the elements of any row.
   This restriction will certainly be violated frequently
   by application programs.
   It would be better to support arbitrary strides
   between the elements in the rows of a matrix view
   and create a condensed matrix if it is necessary
   to pass the matrix to a BLAS function.
   Both gsl_vector and gsl_matrix views
   appear to be restricted to positive strides
   but it would be better if the GSL API
   would support zero and negative strides as well.

4. How should the API conceal the details
   of the representation and the implementation?
   Exposing too many details of the representation
   and implementation may unnecessarily preclude
   some important implementation options
   such as reference counting for example.
   The GSL API should hide the data representation
   and provide methods to initialize and retrieve
   gsl_block, gsl_vector and gsl_matrix attributes.
   But the gsl_block, gsl_vector and gsl_matrix
   struct definitions must be public
   if the methods are to be defined
   by inline functions or C preprocessor macros.

        inline
        int gsl_vector_initialize(gsl_vector* v,
          gsl_block* b, size_t o, size_t n, ptrdiff_t s) {
          v->block = b;
          v->data = b->data + o;
          v->size = n;
          v->stride = s;
          return 0;
          }

        inline
        gsl_block* gsl_vector_block(const gsl_vector* v) {
          return v->block; }
        inline
        size_t gsl_vector_offset(const gsl_vector* v) {
          return v->data - v->block->data; }
        inline
        size_t gsl_vector_size(const gsl_vector* v) {
          return v->size; }
        inline
        ptrdiff_t gsl_vector_stride(const gsl_vector* v) {
          return v->stride; }

        /* alternately */
        /*
        #define gsl_vector_block(v) (v)->block
        #define gsl_vector_offset(v) ((v)->data - (v)->block->data)
        #define gsl_vector_size(v) (v)->size
        #define gsl_vector_stride(v) (v)->stride
        */

   Numerical application programmers should be advised
   that GSL library developers are at liberty
   to change their representation and implementation
   at any time without warning so any application program
   which circumvents these methods
   to access data members directly may break.

   Since the struct definitions are public,
   objects of type gsl_block, gsl_vector and gsl_matrix
   can be allocated on the stack instead of the heap.

        int f(size_t n) {
          gsl_vector vector;
          gsl_vector* v = &vector;
          gsl_vector_initialize(v,
            gsl_block_alloc(n), (size_t)0, n, (size_t)1);
          /* Do something with gsl_vector vector. */
          gsl_block_free(gsl_vector_block(v));
          return 0;
          }

   This can be an important optimization
   for application programmers who are concerned about performance.

Take a look at
the Vector, Signal and Image Processing Library API standard

        http://www.vsipl.org/

and
The C++ Scalar, Vector, Matrix and Tensor class library

        http://www.netwood.net/~edwin/svmt/

These are proposals for standard numerical library APIs.
Each includes a portable reference library
but avoids specifying any details that might preclude
and important implementation.