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