[Beowulf] A start in Parallel Programming?

[Robert G. Brown]

On Mon, 19 Mar 2007, Greg Lindahl wrote:

> On Mon, Mar 19, 2007 at 08:02:40PM -0400, Joe Landman wrote:
>
>> The issue (last I looked last year) with GSL was parallelism.  There was
>> an OpenMP "port" attempt at one time.  No MPI-er-ization (is there a
>> reasonable word to describe the process of converting to MPI?)
>
> Then you should be sticking with BLAS, ATLAS, Scalapack, etc.
>
> I can't believe the GSL people invented a new interface for the one
> numerical interface which is now universal. Not to mention that they
> ignored lots of faster, free libraries like ATLAS and FFTW... what's
> the point of reinventing the wheel badly?

I think that this does them a disservice.  From the BLAS support
section:

   The interface for the gsl_cblas layer is specified in the file
   gsl_cblas.h. This interface corresponds to the blas Technical Forum's
   draft standard for the C interface to legacy blas implementations.
   Users who have access to other conforming cblas implementations can
   use these in place of the version provided by the library.

That is, you can either use their built in CBLAS, which should "satisfy
the needs of most users" (to quote from the same document) or replace it
with e.g. ATLAS-optimized BLAS.  See:

2.2.2 Linking with an alternative BLAS library

The following command line shows how you would link the same application
with an alternative cblas library called libcblas,

      $ gcc example.o -lgsl -lcblas -lm

For the best performance an optimized platform-specific cblas library
should be used for -lcblas. The library must conform to the cblas
standard. The atlas package provides a portable high-performance blas
library with a cblas interface. It is free software and should be
installed for any work requiring fast vector and matrix operations. The
following command line will link with the atlas library and its cblas
interface,

      $ gcc example.o -lgsl -lcblas -latlas -lm



Similarly, they are actually using reimplementations of e.g. netlib code
for fft -- fftpack.  I don't think that they let you link in an
arbitrary external fft library though.  The GSL "philosophy" isn't to
reinvent wheels, but on the other hand it is to reimplement them inside
a single consistent call interface and library.  There are good things
and bad things about doing this (apropos to the extended discussion).  A
single interface means that if you use their data objects for e.g.
matrices and vectors and the like you don't have to screw around with
the details of those objects when passing them through GSL routines the
way you would if you had to pass them between two or three "independent"
external libraries, each with its own expectations for how a vector or
matrix is supposed to be represented in a struct.  However, it also
means that you're stuck with the level of efficiency implicit to those
data objects, and so optimizing a core loop may no longer be possible.
In some cases they actually provide an opaque interface to the objects
that "prevents" one from writing code that directly accesses the data
objects, except that this is a C library so of course it doesn't.

So you can use gsl_vector_get(v,i) to get the ith element of the
gsl_vector v with all sorts of range and bounds checking, or you can
look at:

      typedef struct
      {
        size_t size;
        size_t stride;
        double * data;
        gsl_block * block;
        int owner;
      } gsl_vector;

and do a wee bit of pointer magic to access it directly -- no
"protection" in C (thank god:-).  This is something I actually dislike
in the GSL, especially because the GSL doesn't support tensors of
arbitrary rank (which are important in physics, at least).  Not all
matrices are 2x2.  My proposed "tensor" type would have required that
they redefine vector and matrix, though, and required a cast to make
output typing transparent and I think that this was too daunting.

As always, the coder has full choice in any piece of code as to how they
wish to do something.  The blas is fairly fundamental to the efficiency
of many things that use it by inheritance from within the library
itself, so it is replaceable.  FFTs on the other hand are maybe more
likely to be done just for the purpose at hand and not as part of
another library routine that uses the integrated FFT.  So one can just
not use the GSL library FFT -- use whatever FFT library you like AND the
GSL, just not the GSL FFT.

Overall, I think that the GSL is remarkably high quality code, and it is
certainly aggressively maintained and tested.  It blows away Numerical
Recipes (not difficult, I understand:-) and I'm guessing that it is
rapidly approaching or surpassing the quality of well-known commercial
libraries, e.g. IMSL or NAG, at least in specific areas.  And as always,
it's real advantage in the long run is that it is full GPL (viral).  So
if you, or anyone else, feels that its FFT sucks, well, write a
replacement using the source of YOUR favorite FFT library, integrate the
result with GSL data objects (hopefully without getting TOO irritated by
them) and propose it as an addition or replacement.  GSL contains lots
of additions -- many call routines are shells that can front any of a
whole list of background algorithms.

The entire RNG library is of this form, which is why I use it both for
simulations and for dieharder -- I can test any of 60 well-known RNGs by
just changing an index into this list, and can extend the list to
include another five or ten RNGs of my own or that are otherwise
contributed without changing a line of code that uses the rngs
(including higher level GSL routines that use the specified/initialized
rng to generate e.g. gaussian deviates or exponential deviates).  This
is overwhelmingly useful to simulation folks -- the standard recurring
nightmare is that you are FORCED to pick and use some presumably "good"
RNG (tested if you are lucky on strings of perhaps 10^8 rngs) to use in
your simulation.  10^23 or so rng calls later, you publish a result.

However, your simulation is ALWAYS an RNG tester -- and I mean this
quite literally.  The null hypothesis is that the rng is good, and that
your results are within a sigma or two of where they should be, and that
you know sigma.  However, you DON'T KNOW WHAT THE EXACT ANSWER IS -- if
you did you could compute the p-value of the answer you got and use it
to assess whether or not the rng you used was in fact a good one!

It is therefore very desirable to be able to do a whole computation with
(say) 10^16 rng calls or even more with each of >>several<< "presumed
good" but algorithmically very different rngs and apply a collective
null hypothesis to their answers -- if they all agree within their
mutual sigmas (or so) then it seems more likely that any of these
answers is in fact reliable than it would if one got three or four
answers that were separated by 10\sigma from one another.  The GSL makes
doing this quite easy, and dieharder lets you test candidate generators
much more stringently than diehard used to (and still does for many
people).

So yeah, perhaps it is a shame that the GSL isn't parallelized, but as
full GPL open source, all that stops it from being parallelized is
somebody with a vision or a need to take on parallelization of all or
part of it as a project.  If necessary you can even fork the whole damn
library and call an entirely new version of it gsl_mpi and never look
back.  However, I don't think Brian would be averse to implementing a
transparently parallel version of the library, and there are ways within
the project for libgsl_mpi to be written on top of the GSL (replacing
only those parts of it that need to be replaced) and simply linked in
along with the GSL.

This is unfortunately not that easy to consistently DO -- as the list
has discussed, MPI doesn't really have an ABI yet.  libgsl_mpi would
therefore very likely need to be written and tested with just one "MPI",
and if one was very lucky it would work when relinked with another but
one probably wouldn't be that lucky...

    rgb

>
> -- greg
>
>
>
> _______________________________________________
> Beowulf mailing list, Beowulf at beowulf.org
> To change your subscription (digest mode or unsubscribe) visit http://www.beowulf.org/mailman/listinfo/beowulf
>

-- 
Robert G. Brown	                       http://www.phy.duke.edu/~rgb/
Duke University Dept. of Physics, Box 90305
Durham, N.C. 27708-0305
Phone: 1-919-660-2567  Fax: 919-660-2525     email:rgb at phy.duke.edu