26 April 2007

Linear Algebra in PLT Scheme

Last night, I posted a library which uses BLAS and LAPACK from within PLT Scheme to perform simple numerical linear algebra.

In doing so, I've learned a lot about the FFI for PLT Scheme. The conclusion: it's great! In particular, the "custom _fun types", discussed here are tremendously helpful! Particularly when interfacing to Fortran libraries, where every argument to a function must be a pointer. For example, here's how I make a call to dgesv from LAPACK to solve linear systems:

  (define matrix-solve-many
    (get-ffi-obj 'dgesv_ *lapack*
                 (_fun (m b) ::
                       ((_ptr i _int) = (matrix-rows m))
                       ((_ptr i _int) = (matrix-cols b))
                       (_matrix = (matrix-copy m))
                       ((_ptr i _int) = (matrix-rows m))
                       (_u32vector o (matrix-rows m))
                       (x : _matrix = (matrix-copy b))
                       ((_ptr i _int) = (matrix-rows b))
                       (_ptr o _int) ->
                       _void ->
                       x)))
all the (_ptr o ...) calls allocate the appropriate memory blocks and pass pointers to them for dgesv, and then handle the translation back to the scheme world automatically. I didn't have to write any C code to construct the library.

Up next: C code and PLT FFI for fast manipulation of SO(3,1) group elements (i.e. Lorentz transformations). Stay tuned.