[Haskell-cafe] Large, numerical calculations in Haskell

[Carter Schonwald]

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 <emilhedevang at gmail.com>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 <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
>>
>
>
> _______________________________________________
> 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/05957d83/attachment-0001.htm>