Re: Guile numerical work and uniform arrays

["Harvey J. Stein" <hjstein@bfr.co.il>]

I think uniform arrays is a side issue.  The two things I need for
numerical work are:

  1. to easily (and with little overhead) call numerical code written
     in Fortran & C, and
  2. to easily (and with little overhead) manipulate the data going
     to/from the Fortran & C code.

What this really means is you want a good foreign function interface &
the ability to manipulate foreign data.

The Fortran case is much easier than the C case because the types are
much more limited.  You basically just pass in pointers to blocks of
memory which can be viewed as arrays of the native types - real,
double, integer, etc.  The only tricky part is the calling convention
for strings, where a length argument must also be passed in &
different compilers may have different conventions as to how to do
this.

With the Fortran code you can basically get away with things like
dynamically linking in the raw Fortran library & then calling it
directly from scheme.  All you need really are uniform arrays of the
native Fortran types.  xlispstat has nice ways of doing this, as in
(make-array ... :type c-double), for making a uniform array &
functions for manipulating them.

What's nice about xlispstat is that all the numeric operations are
vectorized.  You can add/subtract/multiply/... vectors, matrices,
& lists.  All operations are elementwise except for operating on, for
example, a vector & a scalar, in which case the obvious appropriate
action is taken - (* '(1 2 3) 2) is (2 4 6).

Unfortunately, this is also a bit of a performance loss.

Two substantial problems with xlispstat's implementation are:

   1. doing operations in place, and
   2. preserving uniformity.

One needs to provide all the numerical primitives in a set! form to
avoid massive garbage creation.  Linking in BLAS gets you some of what
you need, but not necessarily all of it.  Also, when one *does* want
to create new objects, the system must be careful to create them in
the appropriate types - multiplying an integer scalar by a uniform
integer array should yield a uniform integer array.  It shouldn't give
a normal array, or contain bignums, etc.

xlispstat is weak in giving higher level fcns for doing operations in
place, which often makes it much slower than it needs to be.  And,
none of its generic functions preserve uniformity, so you can't even
use these conveniences except by doing the computations & then
coercing the result back to a uniform array of the appropriate type.

However, I think that Luke Tierney might have fixed some of this in
the latest versions.

With the C stuff you get into all the problems of structures &
pointers that don't exist with the Fortran code.

If you want to support numerics with guile then I think the thing to
do is to take inspiration from xlispstat.  You could look at lots of
other packages too (R, Splus, matlab, octave, ...), but xlispstat
might be closest to what's needed in guile, and you could do alot
worse than just doing what xlispstat does.

-- 
Harvey J. Stein
BFM Financial Research
hjstein@bfr.co.il