Re: simplicity vs. efficiency of algorithms in GSL

[Gerard Jungman <jungman at lanl dot gov>]

On Mon, 2008-09-22 at 16:59 -0400, Robert G. Brown wrote:
> 
> For that reason (and for completeness) I'd be happy to see LAPACK
> "swallowed" by the GSL -- either imported and made its own or fronted by
> thinly disguising macro-level GSL calls (that can always be replaced by
> actual local library calls if/when appropriate) that are DOCUMENTED in
> the GSL online manual, with examples and so on. 

I wanted to do this ten years ago. It's kind of a messed-up situation.
Here's the way I see it.

  1) GSL was a militant single-language design.
  2) Ten years ago (and still today, I think) there was no
     reasonable C interface to LAPACK. By reasonable I mean
     something that deals rationally with the row-major
     vs. column-major nonsense, as well as all the
     specialized packing schemes. Also, f2c probably
     does not count as reasonable (i.e. CLAPACK).
  3) License issues must be considered.
  4) We would all like to see ATLAS become a complete solution.
     But currently ATLAS is only complete for BLAS; it has only
     a small portion of LAPACK functionality.
  5) The C++ world has gone hog-wild with half-finished
     linear algebra projects. Normally such things would
     be irrelevant (the world is full of living-dead software),
     but some of these projects are really promising, and
     I worry/wonder how they deal with interface issues.


CBLAS was introduced around that time (~1998), as an actual
documented interface standard (in draft form). That's why it
is possible to link any CBLAS implementation with GSL, and
why the native gslcblas implementation is a separate library.
I insisted on doing it that way so you could drop in any
conforming implementation at link time.

Doing the same with LAPACK was not possible, due to the
lack of any rationally designed and standardized C interface.
Dongarra et. al claim to be working on the newest LAPACK-thing,
whatever it is, which is supposed to have a rationally designed
multi-language interface. I have no idea what state that
project is in, but I haven't seen anything yet, not even
a draft interface standard.

Ten years ago, we kind of gave up, and the "native GSL"
faction won the argument. Personally, I just wanted to
take the fortran LAPACK code, stick it in a subdirectory,
design a rational C interface, and call it done. But
maybe that was a bad idea after all...

At least we have standards for BLAS, and ATLAS is BLAS-complete.
To the extent that LAPACK functionality benefits from BLAS
tuning, it's not a total disaster. The native GSL linear algebra
crap benefits too, since it's supposed to be using BLAS whenever
possible. But somebody still needs to sit down and think through
all the usage cases. What BLAS should you use, which LAPACK, or native
GSL, under what conditions, etc.

In my view, there will be no gsl-2.0 without cleaning up this mess.

As an aside, a similar problem exists with FFTs. It's kinda pathetic
when the GSL designers urge people to "just use FFTW". It's true,
you should just use FFTW. But exactly how should you use it, how do
you make GSL play nice with it, and why then does GSL have FFTs at all?
Clearly GSL should have interfaces to FFTs, but why does it have
a native implementation?

In fairness to the developers (including myself, I share the guilt),
many things were different in 1998. And there was a general fear of
building in any dependencies on anybody else; that can be a risky
way to live.

But I still remember when we implemented GSL native CBLAS.
Oh my god. I never want to do anything like that again. Not ever.


> Another good reason for doing this is that linear algebra is important
> in a lot of other scientific numerical work and could easily be useful
> during the development process of new GSL routines or functions.  Having
> internal, efficient, linear algebra routines would permit the GSL to use
> those routines in other functions -- perhaps in curve fitting, conjugate
> gradient etc -- in a self-contained way with no external dependencies.

Absolutely. This was my number one argument for having a best-of-class
native implementation (by dropping fortran-LAPACK into a subdirectory!).

One of the ingredients in the original GSL kool-aid was the idea that
the sub-libraries should be independent of each other, each devoted
to some task. I pointed out very early on that this was naive, since
you would need to use some functionality to build up other stuff. I
proposed an explicitly designed "level" scheme, making clear which
sub-libraries were foundational, and which had dependencies. People
didn't get it. In the end, we adopted a build procedure where all
the libraries got mashed together into one giant libgsl.a anyway;
so much for independence. [Except for libgslcblas.a of course...]

I do still think that it is fundamentally wrong to have a native GSL
implementation for linear algebra, even given the desire for the
library to be self-contained. What we need is a good standard
interface, which requires people outside GSL to agree on something,
and at least one good implementation, which requires experts to
implement (and none of the GSL developers are linear algebra experts).

We currently have these things for BLAS functions (ATLAS
with CBLAS interfaces). We have LAPACK, in stable fortran
implementations, etc.; now we just need a good C interface.
I hope Dongarra et. al do something. Somebody please do something...


> With a GPL code base, there is no reason for it not to eventually get
> there.  I'd like to see it become a true swiss-army-knife universal
> toolbox -- the only library one ever needs to link to accomplish heavy
> numerical lifting.  But at the moment, it seems to avoid adding things
> like tensors, multidimensional quadrature, full-service linear algebra
> that are in principle (but not always in practice, at least "easy"
> practice) available in other packages.

To achieve this, we have to think bigger than GSL itself. We must
open the door to other implementations and we must conform to
established interfaces, when they exist. Saves us a lot of work too!
What is needed is an integration effort; getting everything we need
under one roof, and usable through one well-designed consistent
interface. As a by-product, I think we will get multi-language
interfaces for free, and the whole world of PyGSL and the rest
will benefit.

Hi-ho.

--
G. Jungman