Re: simplicity vs. efficiency of algorithms in GSL

["Robert G. Brown" <rgb at phy dot duke dot edu>]

On Mon, 22 Sep 2008, Patrick Alken wrote:

> > About this point: I have been wondering where the balance of GSL should
> > be on simplicity vs. efficiency of the code. Of course, even a simple
> > algorithm is better than none, but I think it would be good to offer
> > efficient routines in GSL. If somebody publishes a well written,
> > efficient eigenvalue algorithm, would it get included in GSL? Or would
> > it be better off as an extension library?
>
> I think it depends on how simple and easy to understand it is. The
> algorithm in Golub and Van Loan is relatively simple and someone
> who is interested could read that book and then look at the GSL
> code and understand what is going on fairly easily.
>
> A year or so ago I tried to implement the current LAPACK algorithm
> for nonsymmetric eigenvalues, called the "Small Bulge Agressive
> Early Deflation" algorithm. This algorithm is easily 50-70% faster
> than the current GSL algorithm however it is extremely complex
> and it would take a non-expert several weeks to understand what
> the code is doing.
>
> Since that algorithm is available in LAPACK, a free library, I
> really don't see the need to put faster more complicated algorithms
> in GSL since LAPACK already has the fastest algorithms currently
> available.

My only comment on this is that by in large, from a user point of view I
don't really care what the code is doing as long as it is written by
experts and tested and delivered to me ready to run (this is an
exagerration, sure, but still, one point of a library is to not HAVE to
understand the algorithm to use the -- reliable -- tool).  It isn't,
therefore, the algorithm that matters so much as the front end and in
particular its documentation.

I recently had call to do some linear algebra in code (something I
fortunately don't have to do too often) and discovered that LAPACK is
not, actually, terribly well documented.  It took me many hours and a
lot of googling to find useful and relevant examples I could adapt to my
own needs.

OTOH, the GSL is superbly documented.  From the beginning I've been able
to go into the online manual and get a very accurate idea of what the
API looks like for all its many functions, with at least adequate
examples ranging from code fragments to working programs to help smooth
the way.  I use it so much I went ahead and bought the paper copy of the
manual not so much because I need it (paper, alas, becomes obsolete) but
to help support the project.

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.  To me the hardest
thing about library linear algebra calls is not the algorithm at all --
it is figuring out just where one is supposed to put the matrix being
worked on, the vector(s) being worked on, whether or not one has to call
this routine before (to, say, put the matrix in tridiagonal form) or
that routine after in order to take a half-cooked answer and turn it
into what you are looking for.

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.

The same thing is true in other areas.  There exists at least one fairly
advanced C++ library that does multidimensional (high dimensional)
quadrature, but the GSL does not have it.  I think it should (especially
since C and C++ don't exactly have the same program interface), and
offered to port at least a few routines from Hint with no takers.  The
GSL doesn't manage high dimensional tensors -- it pretty much stops with
2x2 -- but a lot of science is currently done in N>2 dimensional spaces
where having e.g. 4th rank tensors would be very useful.  I
> > personally<< would like the GSL to become the one-stop shop for
scientific and mathematical subroutines, well documented (with plenty of
examples) and exhaustive.

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.

    rgb

> However if you want to implement these as extensions some would
> certainly find that useful.
>
> Patrick Alken

--
Robert G. Brown                            Phone(cell): 1-919-280-8443
Duke University Physics Dept, Box 90305
Durham, N.C. 27708-0305
Web: http://www.phy.duke.edu/~rgb
Book of Lilith Website: http://www.phy.duke.edu/~rgb/Lilith/Lilith.php
Lulu Bookstore: http://stores.lulu.com/store.php?fAcctID=877977