Re: FFTW 3.0 beta available

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

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