Weird results from FFTW -- banding and large roundoff

Yes Steve, I did read the manual and I don't "guess". I have a PhD from a
top university and 20 years in the business. The point is that I made a
hueristic inference and it turned out wrong. And speaking of hueristics,
why even use the "x" in "nx", or similarly for the other axes? Are you not
thinking of Cartesian x, y, and z axes? If so, then which direction is
conventional for the x axis? And then which direction would "nx" reference
the count of? I, for one, don't think it is unambiguous. You describe the
array as "n1 x n2 x n3... in row major order". But the prototypes use
nx,ny,nz, and not n1,n2,n3. Where's the association of nx to n1, or n3 to
nx, or whatever? Is it right alongside?

Here's my suggestion: How about scrapping "nx,ny,nz" altogether for the
function prototypes and using instead "n1,n2,n3" ? Then state simply that
n1 refers to the most slowly varying subscript.

Michael

>On Feb 26, 12:37 pm, "Michael48" <m...@mirametrics.com> wrote:
>> I really do appreciate your work of FFTW, and I understand how array
>> ordering works in Fortran vs C. But that wasn't my point. My point was
>> that your doc explicltly states the parameter passing order as "foo(
nx,
>> ny, nz, ..."  So my mistake was to associate (nx) with the number of
>> columns and (ny) with number of rows. If you are looking at a 2-D image
on
>> a screen, then don't you think that is somewhat hueristic? If I have a
C
>> program, with its multi-D arrays in a native row-major order, and my
code
>> calls your C-library and the docs show the parameter order as nx, ny,
nz,
>> then can you see that there's room for confusion at the implementation
>> end?
>
>There is room for confusion whether we list them as (nx,ny,nz) or
>(nz,ny,nx): that's why you have to read the manual instead of
>guessing.
>
>In English, it is conventional to name the dimensions "x", "y", "z"
>listed from left to right.  It is also conventional in discussing row-
>major order to read the dimensions from left to right, with the right-
>most being the contiguous dimension.
>
>Conversely, it's not reasonable to expect FFTW's data formats to be
>determined by how you happen to display your 2d data on a screen.
>
>The fact is, when you use any library that processes multidimensional
>data, you have to find out how the data must be stored---it is foolish
>to guess.   And we say *over and over* in the manual that the data is
>in row-major order.  In the tutorial.  In the reference section.  In
>the section on multi-dimensional array formats.  As far as I can tell,
>in every place the dimensions are discussed.
>
>I'm not sure what you expect us to do.  If we label the dimensions
>(nz,ny,nx), then at *least* as many people will be confused in the
>opposite direction.  And for people who just read the function
>prototype and don't read the accompanying explanations in the manual,
>nothing we write will make any difference...
>
>Regards,
>Steven G. Johnson
>
>

On Feb 26, 3:31 pm, "Michael48" <m...@mirametrics.com> wrote:
> Yes Steve, I did read the manual and I don't "guess". I have a PhD from a
> top university and 20 years in the business. The point is that I made a
> hueristic inference and it turned out wrong.

When you write "heuristic inference," I read "guess".  The question
is, why did you make a "heuristic inference" from the prototype,
rather than following the explicit documentation?

You read the manual and you missed something that is clearly spelled
out several times.  This happens to the best of us.  But I'm not sure
the blame lies entirely with the manual.

> And speaking of hueristics,
> why even use the "x" in "nx", or similarly for the other axes?

Because 2d and 3d FFTs most commonly correspond to spatial dimensions,
and that is how I'm used to thinking of them.

> Are you not
> thinking of Cartesian x, y, and z axes? If so, then which direction is
> conventional for the x axis?

In row-major order, the first (discontiguous) dimension is
conventional for the x axis, since the dimensions are usually listed
from left to right alphabetically.   How the data is stored in memory
has nothing to do with how the data is displayed on the screen.

> And then which direction would "nx" reference
> the count of? I, for one, don't think it is unambiguous.

I never claimed it was unambiguous.   I don't think that there *is* an
unambiguous naming convention.  I think you have no choice but to read
the manual and see what it defines the data layout to be.

> You describe the
> array as "n1 x n2 x n3... in row major order". But the prototypes use
> nx,ny,nz, and not n1,n2,n3. Where's the association of nx to n1, or n3 to
> nx, or whatever? Is it right alongside?

Yes, it's right alongside.  If you look in the reference manual, it
states that the array is "nx x ny x nz" with "row-major order".  In
the tutorial section, it is also right alongside, stating, "All of
these transforms operate on contiguous arrays in the C-standard row-
major order, so that the last dimension has the fastest-varying index
in the array."

> Here's my suggestion: How about scrapping "nx,ny,nz" altogether for the
> function prototypes and using instead "n1,n2,n3" ? Then state simply that
> n1 refers to the most slowly varying subscript.

That's a reasonable, constructive suggestion.  Thank you.  I don't
think there's anything wrong with calling them "nx, ny, nz", but "n0,
n1, n2" has the advantage of using a nomenclature consistent with the
arbitrary-dimensional case.

Regards,
Steven G. Johnson