Re: FFTW 3.0 beta available

["Steven G. Johnson" <stevenj at ab-initio dot mit dot edu>]

On Mon, 17 Mar 2003, Brian Gough wrote:
> Steven G. Johnson writes:
>  > Is there any interest in calling FFTW from GSL?  GSL's FFT performance is
>  > respectable, but our latest benchmarks (not yet posted) show that FFTW 3
>  > is significantly faster on most machines, as well as being more general.
> 
> In the FFT Chapter of the GSL manual we recommend that people use
> FFTW.  There is no competitor to FFTW and it is GPL'ed, so what more
> can I say ;-)

Thanks!  Sorry; I missed the FFTW reference, since it was at the end of
the section.  What are your plans for the FFTs in GSL, then?  Deprecation?

On your TODO list I noticed the item:

* Sine and Cosine Transforms from FFTPACK

Note that FFTW 3 has DCTs and DSTs (types I-IV, whereas FFTPACK has types
I-III).  (Like FFTPACK, we currently do this by pre/post-processing of
real DFTs, although we may add more specialized algorithms in the future,
or at least some hard-coded transforms of small sizes.  Note that some of
the array passes in the FFTPACK DCT/DST code can be combined.)

One thing that might be quite useful would be if you guys would implement
(or collect from somewhere) the various windowing functions and things
that people use to pre/post-process FFTs for spectral estimation,
etcetera, in DSP.  (We aren't DSP experts ourselves and prefer to stick
with the pure transforms.)

I think the GSL project is ideally suited to fulfill a role much like that
of the GNU project in creating a free Unix replacement.  GNU started by
looking at the existing free tools (such as LaTeX and X11), and
systematically identified gaps in what was available and filled them
(including writing glue programs like Texinfo).  In the same way, I would
suggest that GSL should start with the best existing free numerical code,
such as LAPACK and ATLAS/BLAS for linear algebra, and then systematically
fill in the gaps and create glue if necessary.[*] Like GNU/Linux, there is
also a need for prepackaged distributions that include all relevant
software, documentation, and build scripts.

Reading the GSL documentation, it's not clear to me whether this is your
goal, or whether you are going for a more "Numerical Recipes" approach of
one-stop shopping and re-implementation (instead of achieving a similar
goal by packaging and documenting various code and libraries in one
place).  The NR approach leaves you with a lot of work re-implementing
existing free code, since I'm sure you don't want to emulate NR's record
of poor-to-mediocre implementations.  I get a mixed message when you
present your own code and then say "but go to FFTW if you are
serious," because I know that you guys are serious too.

> Where performance is an issue I'd recommend that people write their
> code to call FFTW directly.  But if someone writes a wrapper to make
> GSL FFT calls use FFTW I'd be happy to make it available somewhere.

Probably, that's not very practical to do efficiently, at least not with
the current GSL interface and with FFTW 3, since in FFTW 3 you need to
create a different plan if you e.g. change the stride.

What might be nice would be to simply document directly the FFTW calls
analogous to the GSL calls (and in general, provide one-stop shopping for
introductory documentation to best-of-breed free libraries).  For example,
calls analogous to:

#include <gsl.h>
...
gsl_fft_complex_wavetable *wt = gsl_fft_complex_wavetable_alloc(n);
gsl_fft_complex_workspace *wk = gsl_fft_complex_workspace_alloc(n);
...
gsl_fft_complex_forward(data, 1, n, wt, wk); /* repeat as needed */
...
gsl_fft_complex_workspace_free(wk);
gsl_fft_complex_wavetable_free(wt);

are (in FFTW 3) (assuming the same double *data as in GSL):

#include <fftw3.h>
...
fftw_plan p = fftw_plan_dft_1d(n, (fftw_complex *) data, 
                                  (fftw_complex *) data,
			       FFTW_FORWARD, FFTW_ESTIMATE);
...
fftw_execute(p); /* repeat as needed */
...
fftw_destroy_plan(p);

Strided transforms are also available in FFTW, of course, but you can
refer to the FFTW manual for the more complicated and less common cases.

Cordially,
Steven G. Johnson

[*] For example, it would be great to have a raw clapack interface,
analogous to the cblas interface.  Otherwise, it is too hard to call
Fortran from C without using e.g. autoconf to figure out the linking
conventions. On top of the raw clapack, one can then build higher-level
abstractions if desired, of course, but it's important to have the raw
calls available for flexibility and because of the pre-eminence of the
LAPACK-style interface.