[Haskell-cafe] Numerical Analysis

[Gregory Crosswhite]

On May 16, 2010, at 4:51 AM, Roman Leshchinskiy wrote:

> You are quite right that vector only supports nested arrays but not multidimensional ones. This is by design, however - the library's only goal is to provide efficient one-dimensional, Int-indexed arrays. I'm thinking about how to implement multidimensional arrays on top of vector but it's not that easy. While repa is a step in that direction, I also need to support mutable arrays and interoperability with C which complicates things immensely.

Indeed;  I eventually decided to hand-roll my own package for dealing with N-dimensional arrays using StorableArray under the hood since I needed to be able to pass them to and from Fortran routines.  It has some nice features such as arbitrary type-safe cuts and slicing (i.e., you can take index 1 of the first dimension and every 3rd index between 2 and 9 of the third dimension and it will returns a new array with a type that has one fewer dimension), folding, converting to/from lists, mutable and immutable versions, and direct pointer access --- which is just the amount of functionality I needed for my purposes.  The problem (and part of the reason I haven't gotten around to uploading it to Hackage) is that the package is not geared around performing efficient computations in Haskell because I have been writing numeric routines themselves in Fortran;  its only real purpose was to let me to examine the arrays in Haskell code.

As an aside, while there are advantages to writing numerical analysis routines in Haskell, it might be better strategy to instead link in something like LAPACK and provide nice wrappers to it in Haskell, since this way you can harness the work of the experts who have spent a lot of time perfecting their code rather than re-inventing the wheel.  One downside of this, though, is that the LAPACK routines only achieve parallelism by making extensive use of Level 3 BLAS routines whenever possible and assuming that they are heavily optimized and parallelized (which they are), so there *might* be cases were a pure Haskell implementation might be able to out-perform them by exploiting even more opportunities parallelism within the algorithm then that provided by the BLAS calls.    Another downside is that using LAPACK requires the arrays to be in pinned memory --- though only for the duration of the LAPACK call.  Finally, I have been experiencing problems linking Fortran and Haskell code on OSX for reasons that I don't understand, since I have no problem on my Linux machine, I have made sure that all the code is compiled for the 32-bit architecture on my OSX machine, and linking C and Fortran programs on the OSX machine does not result in any problems.

Cheers,
Greg