[Haskell-cafe] Large, numerical calculations in Haskell

[Emil Hedevang]

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 <dominic at steinitz.org> wrote:

> Emil Hedevang <emilhedevang <at> gmail.com> 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 at haskell.org
> http://www.haskell.org/mailman/listinfo/haskell-cafe
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.haskell.org/pipermail/haskell-cafe/attachments/20121212/65295ea6/attachment.htm>