Emil, there is an actively maintained set of FFTW bindings on hackage, would that help you? cheers -Carter On Wed, Dec 12, 2012 at 5:03 AM, Emil Hedevang wrote:

Hi A.M. and Dominic,

I have been following the development of Repa for some time, mostly reading blogs and papers about it, not actually writing any code. It's some quite cool stuff the developers are doing.

A discrete convolution is a stencil computation. There are (to my knowledge) essentially two approaches to calculate a discrete convolutions Y = K * X: (1) from the definition, and (2) using discrete Fourier transforms

Approach no. 1 runs in O(nK nX) where nK and nX are the number of elements in the arrays K and X, respectively, whereas approach no. 2 runs in O(nX log nX). If nK is very, very small, then approach no. 1 is fastest, otherwise approach no. 2 is fastest. If nK is somewhat smaller than nX, then some tricks can be applied to obtain something a bit faster than approach no. 2.

To my knowledge, Repa applies approach no. 1 in its stencil computations and discrete convolutions. Therefore it cannot be applied in my case, since I have nK = 10^6 and nX = 10^9 (or something similar), and so I must use approach no. 2 with discrete Fourier transforms (DFT). In Repa, the DFT for arbitrary array sizes is implemented using the slow O(n^2) algorithm, i.e., implemented using the definition of the DFT. If the size of the array is a power of two, then the DFT is implemented using the fast Fourier transform (FFT). Unfortunately, the Repa documentation states that the Repa FFT is approximately 50 times slower than the FFT from FFTW. (And the FFT from FFTW is even capable of calculating fast DFTs of arbitrary size). This essentially forces my to use FFTW for the actual number crunching.

I suppose it would be a very interesting research project to implement FFTW entirely in Haskell. Unfortunately, I currently do not have the time to embark on such a quest.

Regarding the generation of random numbers, Repa implements what they call a minimal standard Lehmer generator. I am not sure if that random number generator is of high enough qualify for my needs.

Cheers, Emil

On 12 December 2012 10:15, Dominic Steinitz wrote:

...

Emil Hedevang writes:

...

Hi Haskell Cafe,

I need to perform very large numerical computations requiring tens of

GB of memory. The computations consist essentially of generation of random numbers and discrete convolutions of large arrays of random numbers with somewhat smaller kernel arrays. Should I use Haskell and call some C libraries when necessary? Should I just write everything in C? Or should I use e.g. Chapel (chapel.cray.com)?

...

Hi Emil,

The first place I would look would be repa http://repa.ouroborus.net/. IIRC it supports discrete convolutions and repa folks seem quite responsive.

Dominic.

PS I am planning to use repa to solve PDEs and address, at least partially, "the curse of dimensionality" with it.

_______________________________________________ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe

_______________________________________________ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe