[MSLUG] Floating point speed of Gambit-C
Bradley Lucier
lucier at math.purdue.edu
Sat Dec 8 13:49:53 EST 2007
Perhaps I should say what my design goals were when designing the
numerical PDE library. This is not really an argument against using
a mixed-language programming model, but I will say that I like to use
Scheme as a single language for the entire system for the following
reason. In the scripting language/low-level language model that is
popular these days in scientific computing, the "interface" between
these languages is fairly fixed, as it is often difficult to achieve
high performance in the scripting language or high flexibility in the
low-level language. In Scheme, I can move that boundary just by
choosing different implementation strategies for parts of the code.
For my numerical PDE class I wanted something that allowed graduate
students to be able to program new (for them) parts of algorithms
that were fairly sophisticated (e.g., multigrid) and be able to do it
in the context of a one-semester project. So nearly all the value in
the system is in the high-level parts, being able to take an
algorithm from a textbook or paper and translate it into code nearly
verbatim (after you struggle to really understand the half-page
algorithm ;-).
The reason that the system is (nearly) as fast as one programmed in C
or C++ is that almost all the floating-point operations in a
multigrid method for solving a finite-element method for elliptic or
parabolic PDEs, say, are in sparse-matrix--vector multiplication, and
that operation is limited by memory bandwidth in either language. So
the fact that the final assembly code for floating-point vector
accesses in Gambit-C--generated code is about 1/2 the speed of that
in C doesn't matter, we're always waiting for memory in either case.
I'll teach the course again next semester, and I plan to have a
preliminary project where students write a complete code for a two-
point boundary-value problem; this will involve rational arithmetic,
symbolic manipulation of polynomials, etc., but also finding the
roots of those polynomials in order to calculate Gaussian quadrature
rules dynamically. The best way to do this is by finding the
eigenvalues of specially-constructed matrices, and a very good way of
doing that is to use the LAPACK code, and, no, I don't want to
rewrite the LAPACK code for that, so I'm going to provide some sort
of FFI interface to that code for the students. I'm not into the
business of rewriting code.
I didn't mind rewriting the level-1 BLAS code in Scheme, however, as
it didn't seem worth the trouble to get a small speedup in the entire
system just to use a dot-product or saxpy written in C and called
from an FFI. I have enough difficulty in getting students to admit
to themselves that, yes, this system is fast (just about as fast as
any expert could have written it, and probably much faster than your
average graduate student could have written it) and it's flexible
(it's only at the end of the course that some students reluctantly
admit that they could not have finished their semester project in
their favorite language, whether C or C++), and that the speed
doesn't arise from level-1 BLAS written in C (because they're written
in Scheme). One point that most students seem to take away from the
class is that multigrid is one hell of a lot faster than conjugate
gradient, and many of them have been using conjugate gradient in
their own projects simply because it's too hard to program multigrid
in C or C++ (at least the first time you try it). So one gets a lot
of speedup simply by being able to program more sophisticated
algorithms.
Another example might be codes for large integer arithmetic. There
is a very simple radix-4 recursive floating-point FFT in the Gambit
runtime; it's half as fast as the same code in FFTW (which won SIAM's
award for a numerical package), and a lot faster than many, many
other FFT codes written in Fortran or C. It's also ~200 lines of
fairly straightforward code translated from a textbook. Using Scheme
in this one-person project allowed me to try a new algorithm for FFT-
based integer multiplication that gets up to about 1/2 the speed of
the FFT multiplication in gmp. I tried to get the gmp folks
interested, but they already knew about the principles the algorithm
is based on (I didn't invent the algorithm) and they aren't
interested; perhaps written in C it would beat their current code,
and these guys are *all* about speed. The next GCD code in gmp
*will* be a variant of what is now in Gambit's runtime library, put
into an already-existing half-gcd framework by Niels Möller of Sweden
after some correspondence between him and myself; I used this
algorithm by Schönhage after some correspondence with him because it
was especially well-suited to Scheme's computational model, but now
it seems that it is faster than any other generally-known GCD
algorithm. (Schönhage might know a faster one, but he doesn't talk
too much about his current work.)
Anyway, I'm going on too long. I use Scheme for numerical stuff
because it allows me to try new algorithms, or more sophisticated
algorithms, that would just be too difficult to play with in other
languages. And I get good enough speed (and the new or more
sophisticated algorithms often more than make up for the constant
factor I might gain by using C instead of Scheme).
Brad
More information about the MSLUG
mailing list