FFTW 3.0 beta available

FFTW 3.0 beta available

[Steven G. Johnson]

Dear GSL folks,

Since we're both providing GPL'ed numerical libraries, we thought you
might be interested to hear that we've released a beta of FFTW 3.0 on our
web site (www.fftw.org).

Version 3.0 of FFTW was rewritten from scratch to improve the speed,
flexibility, and generality of the library and API.  It also adds a number
of other enhancements from SIMD (SSE2/3DNow!/Altivec) support, to
optimized transforms of real even/odd data (discrete cosine/sine
transforms), to true in-place 1d transforms.  See the release notes
(http://www.fftw.org/release-notes.html) for a more complete list of
changes.

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.

Regardless, we'd be happy to hear any comments you might have.

Happy hacking,
Steven G. Johnson

Re: FFTW 3.0 beta available

[Brian Gough]

Steven G. Johnson writes:
 > Since we're both providing GPL'ed numerical libraries, we thought you
 > might be interested to hear that we've released a beta of FFTW 3.0 on our
 > web site (www.fftw.org).

Hi,

Thanks for the information.  

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 ;-)

I plan to keep the existing FFT routines in GSL for backwards
compatibility, and because GSL does not have any dependencies on
external libraries.  

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.

best regards,
-- 
Brian Gough

----------------------------------------------------------------------
Network Theory Ltd           Phone: +44 117 3179309 (UK: 0117 3179309)
15 Royal Park                  Fax: +44 117 9048108 (UK: 0117 9048108)
Bristol BS8 3AL                WWW: http://www.network-theory.co.uk/
United Kingdom               Email: bjg@network-theory.co.uk     
----------------------------------------------------------------------

Re: FFTW 3.0 beta available

[Steven G. Johnson]

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.

Re: FFTW 3.0 beta available

[Gerard Jungman]

I haven't seen Brian respond publicly to your second message,
so I thought I should say some things.


On Mon, 2003-03-17 at 22:12, Steven G. Johnson wrote:

>  What are your plans for the FFTs in GSL, then?  Deprecation?

Deprecation is too strong a term. Anyway, there are two issues. Having
a native GSL implementation is not a bad thing, although it has the
potential to create some confusion. Ideally I would like to see a
mechanism like the gslcblas/external-cblas one, where GSL provides
an interface but the user can choose whether to use the native gslcblas
library or an external cblas implementation, at link time.

So, what we really need is a set of interfaces which are general
enough to encompass FFTW and the native GSL implementation. This
probably means we can just use the FFTW interfaces as a base.
But there will almost certainly be some issues to sort out.
One the one hand I feel a little saddened that the FFTW interface
hides everything behind the (essentially) opaque fftw_plan object;
for instance, I would have hoped that the plan would be separate
from the data. Maybe I misunderstand the philosophy.
Anyway, it looks like any general interface will
have to have the same design. On the other hand, this sort of
encapsulation makes it easier to use and to wrap; it is not
a big deal.

Furthermore, I don't see any reason to hide or deprecate the current
GSL interfaces; they can be used by the more general wrapper. The
small amount of confusion which this will create can be offset by
good documentation, including examples of how to use an external
implementation (FFTW).


> 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.)

This is a good idea, as long as we have a reasonable plan for what to
include and what to leave out. The hardest part is deciding what to
leave out and sticking with that decision.


> 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.[*]

I absolutely agree that we need a "GNU Scientific Environment".
GSL itself is really intended to be one component of such an
environment (although a central one because of the need for consistent
interfaces), rather than the only entry point for it.
Of course, we can't claim that this was clear in our
minds from the start; GSL has evolved since it began, more than
seven years ago. So GSL still needs some real work to realize
this component-wise view.

At some level, we knew that we would have to implement
a bunch of bits and pieces, and we needed a loosely organized project
in which to drop them, so that those bits and pieces could start
being used and debugged, etc. That is certainly how I see GSL-1.x.
In 1997 it was not at all clear what a more general environment might
look like, and we needed to start somewhere. There was some idea at the
time (still expressed in the documentation and design document) that it
should be a very loosely coupled set of libraries with a strong
"numerical recipes" flavor. This bothered me from the start; as soon
as you realize that you need to use some library components in the
implementation of other library components, the over-simplification
of loose coupling evaporates. I think we reached some compromise
over this, but it is far from ideal.

If we are going to be serious about the "GNU Scientific Environment",
then we will need to have some concrete and reasonable notions of
how things do get coupled together. I hope that the GSL-1.x experience,
especially with BLAS and now with FFTW, will be enough to get us going
in the right direction. The other big external object is LAPACK...
see below.


> Like GNU/Linux, there is
> also a need for prepackaged distributions that include all relevant
> software, documentation, and build scripts.

This is a really big deal, and it will have to be a major force
in the design of any components and wrappers. These things have to
work well from top to bottom, from the numerical aspects all the
way up through the build and install process.


> 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).

I hope the above at least partly answers this. In short, GSL-1.x is a
mixture of the two, with a heavy lean toward "numerical recipes"
because of the history of its development. Parts of GSL will continue
that way in perpetuity, like the special functions, which are kind of
naturally self-contained. But other parts will have to grow up and
learn to talk with the rest of the world.

The goals for GSL-1.x were to fill a shopping cart with generic
numerical tools, solving as many local problems as we could. The
main goal for GSL-2.x must be to organize all these tools together
and make a real attempt to solve the global problems. Of course,
this is the step that kills a large fraction of software projects,
so we will need to be well-prepared to succeed.


> 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.

Again, it's a mixed bag. Some of GSL is straight re-implementation of
high quality code. Some has the look-and-feel of "numerical recipes",
but with improved implementations. And, of course, there are a few
embarrassing shady corners (for example, I have a strong urge to
remove the implementation of the confluent hypergeometric function;
objectively, that code  is a failed experiment).


> 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.

I know what you mean. I try not to make such statements. I think this
type of statement should be taken as shorthand for the following:
"GSL ffts are a well understood re-implementation of high quality
code, and you should feel good about using them. However, we need to
do a better job in interfacing the GSL world to high-performance
implementations such as FFTW. Until we get around to doing this,
and if you really need the extra performance and flexibility,
you have to solve the interface problem on your own. We don't
feel good about this either, but it will take time for us to
do it right."

Long-winded, but accurate.


> > 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.

Yup. On closer examination, these problem domains always turn out to be
more complex than we thought. That's why we will need better interfaces.


> 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).

Yes, it would be good to help people explicitly. Like much of the
documentation, I'll defer to Brian (I admit, I am a real problem-boy
when it comes to documentation).


> [*] For example, it would be great to have a raw clapack interface,
> analogous to the cblas interface.

I know Brian looked into this. Based on what he learned we decided
to postpone a clapack companion for the cblas interface. As I recall,
the main issue was how to deal with the workspace management issues,
which greatly complicate things. On the other hand, this may have
been a serious issue only for wrapping the fortran implementation;
I don't recall the details. That was about 4 years ago.
This needs to be looked at again, and it will have to be
solved as part of GSL-2.x.


Your original question:
> Is there any interest in calling FFTW from GSL?

In case Brian's response wasn't clear, the answer is definitely 'yes'.
We have to make it possible to do this in an easy way, preferably
allowing the user to defer the decision to link time, like with
cblas. We'll have to see what we can do.


Thanks for the message, and for FFTW.


-- 
Gerard Jungman <jungman@lanl.gov>
Los Alamos National Laboratory

Re: FFTW 3.0 beta available

[Brian Gough]

Yes, sorry for not replying yet.  I'm busy during the week and so
didn't have time to compose an adequate reply yet.  I will do so at
the weekend.

I try to CC the list for any replies so if you don't see anything it
means I haven't had time to reply.

Gerard Jungman writes:
 > 
 > I haven't seen Brian respond publicly to your second message,
 > so I thought I should say some things.
 > 

Re: FFTW 3.0 beta available

[Steven G. Johnson]

Gerard, thanks for your thoughtful response.

On 20 Mar 2003, Gerard Jungman wrote:
> One the one hand I feel a little saddened that the FFTW interface
> hides everything behind the (essentially) opaque fftw_plan object;
> for instance, I would have hoped that the plan would be separate
> from the data. Maybe I misunderstand the philosophy.

The main idea was that (a) in most high-performance applications (as far
as I can tell) you are usually transforming the same array over and over,
(b) if you need to transform another array of the same size, creating a
new plan once the first exists is a cheap operation, and (c) see below.

FFTW 3 does actually allow you use a given plan for a different
input/output array, but you have to satisfy certain constraints, like you
can't use an in-place plan for an out-of-place transform, the strides must
be the same, etcetera.  These constraints are important for performance,
it turns out, since a plan can use tricks very specific to a certain
stride, etcetera.  Also, if there is any chance that the program could be
linked to a SIMD version of FFTW (e.g. using SSE2), then the arrays have
to be equally aligned modulo 16 to those the plan was created
for...otherwise, the program could crash (SIMD is picky).  Because these
constraints require careful attention to the documention, we decided to
put the ability to apply a plan to a different array in our "guru"
interface, bedecked with warnings to scare away the casual user.  =)

(I need to add this to our FAQ.)

> Anyway, it looks like any general interface will have to have the same
> design. On the other hand, this sort of encapsulation makes it easier
> to use and to wrap; it is not a big deal.

One should also consider the possibility that one may not always want to
create wrappers in GSL for existing free codes, depending upon the problem
and upon the quality of the available software...one may just want to say
"we don't address this problem in GSL, because we recommend libfoo...but
here is how you use libfoo with GSL data types."

> I know Brian looked into this. Based on what he learned we decided
> to postpone a clapack companion for the cblas interface. As I recall,
> the main issue was how to deal with the workspace management issues,
> which greatly complicate things. On the other hand, this may have
> been a serious issue only for wrapping the fortran implementation;

I think you misunderstood me somewhat.  I don't mean creating a gsl_lapack
wrapper that uses underlying lapack code with gsl datatypes, at least not
right away.

I mean making a CLAPACK interface that just wraps the Fortran LAPACK in
the most direct possible way (ideally collaborating with the LAPACK
folks), like the CBLAS standard (not gsl_cblas) does for BLAS.  A "raw" C
lapack interface is needed to allow people maximum flexibility in how (and
whether) to use a higher-level interface and types; but right now, calling
LAPACK from C requires some knowledge of inter-language tricks.  
Probably, a raw CLAPACK (not to be confused with the f2c'ed LAPACK of
the same name), should be a logically separate project from GSL.

Then, on top of that, libraries like GSL could then make higher-level
interfaces that use the gsl types, automatically allocate workspaces of
the correct size, etcetera (there is no technical problem with this, I
have done it myself in my own code).

I think you'd be crazy to re-implement LAPACK in C instead of writing
wrappers; it's a lot more effort (repeated each time a new version of
LAPACK comes out) for no benefit (users don't care what code is doing the
underlying computation, only about the interface).  I hope you're not
considering that.

Steven

Re: FFTW 3.0 beta available

[Brian Gough]

GSL is intended as a free alternative to general purpose numerical
libraries such as NAG/IMSL/Numerical Recipes.  Historically this
accounts for the majority of users.  Consequently GSL's main design
goal is ease of use for the casual user, rather than performance.

I see high-performance libraries as being in a different part of the
ease-of-use/performance trade-off space from GSL.  So, for example, I
plan to keep the FFTs in GSL, since they fill a different need from
FFTW.

I'm not personally planning to expand the scope of GSL into an
umbrella project, since we're already resource limited - at the moment
I would just like to get more people involved in what we have.

Regarding LAPACK, I did once look at the possibility of making it more
accessible from C but basically decided that it was not worthwhile due
to (a) issues with row/column-major storage (b) licensing of LAPACK
itself and its non-free documentation.  I decided to leave the LAPACK
problem for someone else to solve.

-- 
Brian Gough

----------------------------------------------------------------------
Network Theory Ltd           Phone: +44 117 3179309 (UK: 0117 3179309)
15 Royal Park                  Fax: +44 117 9048108 (UK: 0117 9048108)
Bristol BS8 3AL                WWW: http://www.network-theory.co.uk/
United Kingdom               Email: bjg@network-theory.co.uk     
----------------------------------------------------------------------