[Help-gsl] 2d FFTs, DCTs, etc.

James Bergstra james.bergstra@umontreal.ca
Fri May 12 17:01:00 GMT 2006

Watch out, it's a long one.  Perhaps un-necessarily so... but nonetheless:

> *  Wrapping all of the functionality of FFTW in GSL (all interfaces, 
> wisdom saving, etc.) would be, I feel, too complex for a general library 
> like the GSL, but then again that's not my call, and it's not what 
> you're proposing either (is it?).
Admittedly, what I had in mind was in fact a third interface, that adds
considerable code, and no new functionality.  I don't believe that it should
replace the gsl fft api.  I like the extension system, and I think this is the
sort of thing that would make a good extension.

My goal was an API that is as easy to call as gsl's, and that admits an FFTW
implementation, but doesn't require it ( gsl has a priority of being


I believe that even an interface that is technically superfluous, can have
practical, and pedagogical value.  An FFT interface that is *very* suboptimal,
such as:

    my_fft(data_in, data_out);

can be accompanied by a comment along the lines of "this interface is
suboptimal, refer to the source code for example use of the underlying API, and
tips for enhancing performance."  

New User can get his program up and running quickly with this interface (because
it is practical), and then when optimization becomes an issue, he or she can
cut-and paste the source of my_fft(), and slowly do his own loop-hoisting,
memory allocation, etc.  (Because the interface source is in a sense, a lesson
on how to circumvent the interface itself!) This removes the barrier-to-entry
erected by complicated interfaces such as the advanced FFTW interfaces (and
LAPACK, but that was another long-running debate a few months ago!).


> *  Your interface is simpler, but then it's entirely incompatible with 
> the current GSL FFT routines and I don't know whether Brian Gough and 
Right, this is a pretty fat whole in my plan.  If someone asks for a gsl
implementation of a 2d-fft, my api will have to throw an error.  As it stands, 
FFTW is the way to go for multi-d parallel transforms.  The question is how easy
can we make the transition, for people who have invested in the gsl data types.

> the other GSL developers would be willing to change the API.  You would 
> probably still have to make a few choices though (so you'll have to keep 
> implementation details in mind when all you want is a DFT) and I'm not 
> sure whether we would be able to get away with not using saved wisdom 
> files as calculating plans takes some time (from what I've read).  I 
> fear things will have to get complicated.
I don't think it's as bad as you think.  You can use FFTW to implement the gsl
API for example.  The loading of wisdom files is optional, and the search can be
limited by choosing the FFTW_ESTIMATE flag.  By estimating, and using no wisdom,
FFTW uses a default algo that can depend on your alignment and the input size.
Even in this default case it's still fast.  I just wanted to set up the api so
that it's easy to optionally "optimize" later by changing a flag from
FFTW_ESTIMATE to FFTW_MEASURE and adding a load-wisdom statement somewhere.

> one post about FFTW.  This leads me to theorize that someone in need of 
> a smoking fast FFT routine probably doesn't need a generic math library 
> and is happy with a specialized package.  On the other hand if most GSL 
I dispute this, even someone who is happy with a specialized package may balk
at debugging code that calls a function like this:

fftw_plan fftw_plan_many_dft_r2c(
    int rank,
    const int *n,
    int howmany,
    double *in,
    const int *inembed,
    int istride,
    int idist,
    fftw_complex *out,
    const int *onembed,
    int ostride,
    int odist,
    unsigned flags);

It smooths/speeds things appreciably to start with a working example or
demonstration of how to do it.  An API would provide the mapping of how to
convert gsl_matrix and gsl_matrix_complex arguments into this stuff.

> So I was thinking of patching up GSL's FFT routines (by adding 
> gsl_vector functions) if the developers of GSL don't have a good reason 
> for the asymmetry in the FFT interface (the wavelet interface also lacks 
> float functions) and add a DCT. 
I'm sorry that I distracted the thread from this (original) very good point. 

> As for libmsl, if what I have in mind 
> still interests you I'd be happy to host it there.  I haven't released 
> code until now so I'm not familiar with the use of autotools and the 
> general procedure of polishing something for release.  By the way, home 
> come libmsl isn't mentioned in the GSL site?  Is it supposed to be a 
> separate library of extensions or extensions to be incorporated into GSL 
> at some point?
It's got nothing to do with the GSL, it's just a little collection of relatively
stable code that I wrote in pursuit of my masters.  It was on the GSL extension
list at one time, but it wasn't always the *clean and well-organized* collection
that it is now.  I asked BJG about removing the link while I cleaned it up,
and I guess it hasn't made it's way back on.  Maybe some day it will be back :P

( Actually, I'd rather that GSL get the cvs server I was plugging earlier today.
I think that's a better idea than libmsl. )

James Bergstra

More information about the Gsl-discuss mailing list