[Haskell-cafe] Numerical Analysis

Pierre-Etienne Meunier pierreetienne.meunier at gmail.com
Sun May 16 12:52:57 EDT 2010


> 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.

I understand. What complicates it even more (at least in what I imagine) is that C uses the same syntax for multidimensional and nested arrays, and I do not believe that for instance GHC's FFI allows for array types such as "int x[19][3]".

Data.Ix allows for indexing arrays with (Int,Int), for instance, but with costly index conversions, constructing lists in the memory for instance. Also, sometimes the user of these libraries may find a better bijection between N and N^2, better in the sense that for his particular problem, a slightly more complicated arithmetic sequence of operations for computing the index would optimize cache misses better.

Would it be for instance possible however to generate a default bijection using preprocessing (template haskell ?) ?

> That said, if all you need is a matrix it's quite easy to implement the necessary index calculations yourself. Also, since you are working with numerics I highly recommend that you use either Data.Vector.Unboxed or Data.Vector.Storable instead of Data.Vector as boxing tends to be prohibitively expensive in this domain.

I'm actually thinking about rewriting parts of my code now ! It seems indeed that given how the algorithms are specified for instance in "numerical recipes" or the papers in its bibliography, unboxing is the only option, since there is no such thing as "boxing" in cobol and fortran ;-)

I was also wondering about how to do linear algebra : an infinite number of types would be needed to express all the constraints on matrix multiplication : we need types such as "array of size m * n". Is there a way to generate these automatically with for instance template haskell (again ! But I know nothing of template haskell, neither, sorry !)

Cheers,
PE



More information about the Haskell-Cafe mailing list