[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