GSL, VSIPL, and the ultimate numerical library

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

Robert W. Brewer wrote:

> I just read a good chunk of the year 2000 list archives.
> I've been following the VSIPL forum work on and off
> for the past couple of years, and I was pretty excited
> when I saw the announcement of GSL 0.7 on freshmeat.net,
> because it sounded like the GSL might have already
> beaten the VSIPL to the punch so to speak.
> It appears that a lot of good work has gone into GSL so far.
> 
> I've been thinking about what I would want
> in a grand numerical library.
> The things that I myself would probably use
> a good numerical library for are
> some real-time speech processing under Linux
> and maybe as a basis for a neural networks library.
> Performance is important for both of those, otherwise,
> for offline use, I would just use Octave or Matlab.  
> If the Intel Signal Processing library ran under Linux,
> that would be great but it doesn't
> and it is also a closed-source library
> which I'm generally against. 
> 
> The TASP VSIPL code is a pretty good choice
> but it looks hard to maintain
> since there is so much code replication
> to support the different data types.
> I recall an early version the Hughes VSIPL code
> which used a preprocessor template scheme
> that seemed pretty slick.
> It looks like GSL has a similar problem.
> Many functions are written 3 and 4 times
> just to support the different float and int types
> the user might want to instantiate.
> And after all that, it still can't support neat things
> like arbitrary precision numbers,
> while a templating scheme might be able to.
> 
> To boil it down, the library that I would like to use
> would be (most important first):
>    high performance
>    generic to allow optimizing under the hood 
>       and advanced algorithms
>    templated for data type choices
>       and easier maintenance open source
>    easily tested for correctness with an automated suite
> 
> I think the VSIPL API has a lot of good ideas
> for the  functions that it supports
> but the reference implementation would be easier
> to deal with if it were templatized.
> Maybe it's my personality,
> but making the same change in more than one place
> just wears me out.
> 
> The GSL team has implemented a large breadth of functions,
> many of which have no equivalent in the VSIPL specification
> but which could probably benefit from some of its principles.
> Since both the TASP reference implementation and GSL
> are open source,
> it might be possible to gradually combine them,
> as Randall Judd seems interested in doing.
> Many GSL functions not present in VSIPL might be wrapped
> to present an interface similar to the VSIPL API.
> Then gradually the underlying GSL code could be moved
> into the wrapper to increase performance,
> better conform to VSIPL, etc.
> 
> It might also be nice to go back to the technique
> in the Hughes code of implementing most of the functions
> in terms of other VSIPL/GSL functions
> or lower-level helper functions whenever possible.
> This could have the advantage of allowing
> a mostly optimized version for a specific hardware platform
> by just optimizing the lowest-level vector routines
> and helpers for that platform and allowing
> the more complicated routines to ride on top and still benefit.
> The hope would be that, say, 80% of the Altivec
> (or Intel's SIMD) performance gain
> could come from  writing Altivec assembly
> in only, say, 20% of the functions.  
> I mention this because the neat thing, to me,
> about having a nice numerical library is being able to do some
> impressive real-time signal processing tricks on a desktop PC.
> Intel's MMX and Streaming SIMD Extensions,
> AMD's 3D Now instructions and PowerPC's Altivec technology
> all allow for amazingly fast computation
> but since standard compilers can't take advantage of them,
> writing assembly code with a high-level API is the way to go.
> But the idea of writing Rob's Gee Whiz Altivec Library
> is not very appealing to me, first because VSIPL and GSL
> have put a lot of thought into how to do it right
> and second because I'd rather not start from scratch
> if I switch to a non-Altivec machine.
> 
> I guess in the end I can understand
> why Tisdale is so interested in the vision and goals of GSL,
> so it is easier to decide if it is the same direction 
> that he would like to go.  I am in a similar quandary,
> there are several interesting-looking projects
> with various tradeoffs, and looking into the future a little
> would help me bind to one or decide to start my own.

Hi Rob,

The Vector, Signal and Image Processing Library (VSIPL)
has a limited scope.  But Randy Judd and I are both Physicists
so we are interested in general purpose numerical computing
as well as digital signal and image processing.
My hope was that the GNU Scientific Library (GSL)
would not be restricted to linear algebra
like the Basic Linear Algebra Subroutine (BLAS) library
and the Linear Algebra PACKage (LAPACK)
and that it would subsume the VSIPL API standard.

I implemented the VSIPL Portable Reference Library
that you got from Hughes Research Laboratories, Inc.
(now HRL Laboratories, LLC).
The Portable Reference Library was designed to be
easy to read, understand and maintain
because there wasn't any money to pay a programmer
to maintain the Portable Reference Library.
We actually began work on the Portable Reference Library
at the same time that we began work on the VSIPL
Application Programmer's Interface (API) specification
so that VSIPL Forum members could test their proposals
using the Portable Reference Library
before they presented them to the VSIPL Forum.
I suspect that the Portable Reference Library
has already served it's purpose and I suppose
that is the reason why nobody cares
whether it is ever released to the public or not.

I used a modified GNU genclass shell script
which also called the m4 preprocessor
to convert prototype files (containing what you call templates)
into ANSI C header and source files.
It was never my intention to release the prototype files
to the public but I did distribute them
to the members of the VSIPL Forum
along with the header and source files produced from them.

Vector, matrix and tensor classes are not container classes.
They are arrays of NUMBERS and NOT arbitrary types.
The decision to use templates, prototypes
or any other GENERIC method is an IMPLEMENTATION DETAIL
best left up to the LIBRARY DEVELOPER.
If you need a generic implementation,
you will need to negotiate that with the library developer
which, I guess, is what you are attempting to do here.

Randy Judd has implemented a high performance version
of the VSIPL API standard targeted
for a restricted class of computational platforms --
workstations (including personal computers).
He elected not to use any generic method
to implement his library and his implementation
of the VSIPL API standard happens to be similar
to the current implementation of the GSL.
If they were just a little more similar,
they could be INTEROPERABLE.  That is to say that
you could call GSL functions from Randy's VSIPL implementation
and you could call VSIPL functions from the current GSL.

If I could find someone to pay me to implement
a high performance version of the VSIPL API today,
I would probably write a perl script to replace
the old genclass shell script and accept templates
from a prototype file to emit ANSI C header and
source files for each supported type.
I would certainly try to leverage off of any
existing high performance numerical libraries
appropriate for the target platform.

What would I want in a grand numerical library?
Well just open up your copy
of Numerical Recipes in C
and read the table of contents.
I think that would be a very good place to start.
The problem with Numerical Recipes in C
is not that the implementation is so bad
but that the API makes it impossible
to re-implement the library.
And that's exactly what's wrong with the GSL as well.
No matter how clever and gifted the programmer,
it is practically impossible to implement
a high performance numerical library
that will port everywhere.
If you aren't very careful about how you design the API,
you will paint yourself into a corner
from which you will not be able to escape
without changing your API and breaking
existing application programs that depend upon it.

What we really need from the GSL designers
is some documentation of their rationale.
Something that we could use to determine
whether or not the GSL is appropriate for our purposes.
We need answers to very simple, basic questions.

1. What is the target platform?
2. Who are the target users?
3. What is the application domain?
4. What are the priorities?
5. What is the expected life of the library?
6. Who will provide support for the library?
7. What are the plans for the future?
   etc.

Otherwise, we will be obliged to wait and see
what sort of creature the GSL evolves into.