Huang Jianshi wrote:

Hi,

I'm new to Gambit. I read the benchmark and Gambit performs on the top of the list. I want to ask if there is any tutorial/manual/ examples on showing how to write efficient code on number crunching task, especially matrix operations and floating point operations.

Brad Lucier is the main number crunching specialist here. He has some code for matrix operations on his web page (though I don't know how current).

If you want fast handling of single floats (as opposed to vectors or matrices), garbage collection is an issue since floats are normally boxed; you can get around this by using (declare (flonum) (not safe)) in the relevant sections, and/or ##flonum.{*,+,-,/,..} operators. (Floats are still boxed across real function calls, though.)

If you have sparse matrices or want to transparently use multiple cpu's or a cluster without decomposing the work into smaller tasks yourself, the PETSc library may be for you. I did start implementing bindings some time ago but have suspendet that work for now. (I did also plan an infrastructure and did some coding experiments for working with big time series (in disk files) efficiently, but that project has then been cancelled by the customer.)

Christian.

On Apr 19, 2007, at 1:36 PM, Huang Jianshi wrote:

tutorial/manual/ examples on showing how to write efficient code on number crunching task, especially matrix operations and floating point operations.

Christian gave good advice. You'll almost need to write C-style code and declarations (removing the flexibility that is otherwise built in to Scheme).

Start with

(declare (standard-bindings) (extended-bindings) (block) (not safe))

These mean

(standard-bindings): All R4RS/R5RS procedure names point to the standard R4RS/R5RS procedures (i.e., you haven't redefined + to something different) (extended-bindings): The same thing for the Gambit-specific procedures. (block): If you haven't set! a variable in the file in which it is defined, its value is never changed via set! in any other file. (This is why you'll see things like

(define foo ...) (set! foo foo)

in some of the scheme files in lib; this tells gambit that even though the file has (declare (block)), this particular variable can be redefined in other files.) (not safe): You don't want Gambit to insert checking code to save yourself from, e.g., (car 1) crashing gambit.

Then you should use fixnum and/or flonum-specific numeric operations as appropriate. You can put a (declare ...) form anywhere that you can put a (define ...) form, so

(let () (declare (flonum)) ...)

means that all numeric operations in that "let" have flonum arguments and results, and tells gambit to just use the usual floating-point operations (so (sqrt -1.) is +nan.0 instead of +1.i). Similarly for (declare (fixnum)). To make it easier to use these declarations in simple expressions, I usually define the macros:

(define-macro (FLOAT . rest) `(let () (declare (flonum)) ,@rest))

(define-macro (FIX . rest) `(let () (declare (fixnum)) ,@rest))

so one can simply say (FIX (+ i 1)), for example. There's less need for this given the fx+, etc., operations, but it does allow you a bit more flexibility for testing, etc. (one can simply redefine (FIX . body) as (begin ,@body).

The gambit runtime system catches all interrupts and handles them itself, and it checks for interrupts at the same time it checks for stack overflows. (See POLL macros in the C code generated by gsc.) Sometimes it pays to disable those checks in tight loops with

(declare (not interrupts-enabled))

If you do this in a loop that allocates stuff on the stack, then gambit will crash, since there will be a stack overflow that won't be caught. Also, you cannot interrupt a loop that has interrupts disabled, and this can be annoying. In 4.0b22, Marc uses the __builtin_expect built-in function in gcc-3.0 and later to tell gcc that these POLLs are unlikely to be taken, so maybe POLLs will slow down loops less than in previous versions of gambit.

Finally, the Gambit-C compiler generates code that boxes flonums across any function call, (if ...) boundary, etc. And loops are function calls in Scheme, so it boxes flonums across loops. Sometimes allocating the 12-16 bytes (depending on whether Gambit's running on a 32-bit or 64-bit machine) for a boxed flonum takes longer that to compute the value of that flonum.

If this is a real problem, I tend to lesson it for a given loop by unrolling the loop a number of times by hand or, if you're accumulating a sum into a single flonum result, I may allocate a f64vector of length one, read its value at the beginning of the loop, and then store back into it at the end of the loop.

All that being said, I tend to write numerical code by identifying which loops are computationally intensive, optimize them, and then use the late-binding flexibility of normal Scheme semantics to make the rest of the code easier to write and maintain. It's not so easy to find out where your code is spending it's time, since 'gcc -pg' is pretty useless (each file is compiled to one C function). If necessary, I end up configuring gambit with

env CC='gcc -ftest-coverage -fprofile-arcs' ./configure --enable- single-host --enable-gcc-opts

and then run gcov on the output. The problem with doing this is that you need to understand the macros used in the C code generate by Gambit, the so-called Gambit Virtual Machine (GVM), and how to associate that code with lines in your code. (It's not hard, GVM is a simple machine with a few registers, a stack, and the usual abstract operations.) Or, you can compile the scheme file with

gsc -track-scheme file.scm

and then gcov will report the number of times each line of C code is executed. (I just did this for an image-processing program I'm writing, and found out that all source files should be in the same directory or gcov gets confused.)

You can find code for numerical solution of PDEs at

http://www.math.purdue.edu/~lucier/615/software/index.html

that uses all these techniques in one place or another. Search for (FLOAT ...), (FIX ...), and (declare ...). Understanding this code is more complicated because it uses the Meroon object system, but the basic numerical tricks are visible.

Brad