Clarification on intended support for negative strides in FFTW

[haiwenheaven]

Hello,

I am currently exploring and comparing the capability boundaries of Intel oneMKL and FFTW for strided transforms, especially transforms using negative input or output strides.

During this investigation, I found that negative strides appear to be only partially covered and not always handled consistently in FFTW. In particular:

• Some negative-stride cases appear to cause the FFTW planner to continue indefinitely while consuming increasing amounts of memory, eventually exhausting the available memory;
• The existing correctness tests do not appear to systematically cover negative strides;
• Some negative-stride cases expose bugs in the test infrastructure.

I understand that not every negative-stride layout, particularly an in-place layout, is necessarily valid or expected to be supported. Before preparing patches, I would like to clarify the intended support boundary:

1. Are negative strides intended to be supported by the advanced, guru, and guru64 interfaces when the supplied memory layout is otherwise valid?
2. Do you consider negative-stride support, or bugs exposed by negative strides, important enough to address upstream?
3. Would regression tests and focused fixes for these cases be welcome?
4. Are there specific negative-stride combinations that should be treated as unsupported?

If this is considered relevant, I plan to prepare focused patches that add deterministic negative-stride regression tests, fix test-infrastructure problems separately, and address confirmed FFTW implementation bugs in individual changes.

Thank you.

[matteo-frigo]

Negative strides are supposed to work and in fact they are used internally by the "half-complex" real-transform format. Do you have a case where they don't work?

[haiwenheaven]

1. Case

#include <fftw3.h>
#include <stdio.h>

int main()
{
    fftw_iodim transform_dimensions[2] = {
        {3, -12, -4},
        {4,   3,  1},
    };
    fftw_iodim batch_dimension = {2, 32, 22};

    fftw_complex *input_storage = fftw_alloc_complex(128);
    fftw_complex *output_storage = fftw_alloc_complex(128);

    fftw_plan plan = fftw_plan_guru_dft(
        2,
        transform_dimensions,
        1,
        &batch_dimension,
        input_storage + 32,
        output_storage + 32,
        FFTW_FORWARD,
        FFTW_ESTIMATE);

    fftw_destroy_plan(plan);
    fftw_free(input_storage);
    fftw_free(output_storage);
    return 0;
}

2. Analysis

2.1. rank-geq2

The case above entered the rank-geq2(reg_id=0) solver. And split

transform tensor:
    {n=3, is=-24, os=-8},
    {n=4, is=  6, os= 2}
vector tensor:
    V0 = {n=2, is=64, os=44}
placement:
    out-of-place, input buffer -> output buffer

into cld1(Length-4 FFTs Along D1)

transform:
    {n=4, is=6, os=2}
vectors:
    {n=2, is=64,  os=44}
    {n=3, is=-24, os=-8}
placement:
    out-of-place, input buffer -> output buffer

and cld2(Length-3 In-Place FFTs Along D0)

transform:
    {n=3, is=-8, os=-8}
vectors:
    {n=2, is=44, os=44}
    {n=4, is= 2, os= 2}
placement:
    in-place on the output buffer

2.2. indirect-transpose

cld2 entered the indirect-transpose solver, and partitioned into:

transform:
    {n=3, is=-8, os=-8}
vectors:
    {n=2, is=44, os=44}
    {n=1, is= 2, os= 2}
placement:
    in-place on the output buffer

and 3x3 tile:

transform:
    {n=3, is=2, os=-8}
vectors:
    {n=2, is=44, os=44}
    {n=3, is=-8, os=2}
placement:
    in-place on the output buffer

The complete 3 x 3 tile is transposed, while the remainder column is handled separately.

2.3. indirect-before/indirect-after

For the 3x3 tile problem, the planner tries the indirect-after and indirect-before solvers in that order. In principle, this problem is better suited to indirect-after; however, it is rejected by indirect-after and instead accepted by indirect-before.

After entering indirect-before, the solver transforms the 3x3 tile into the following problem:

transform:
    {n=3, is=-8, os=-8}
vectors:
    {n=2, is=44, os=44}
    {n=4, is= 2, os= 2}
placement:
    in-place on the output buffer

This problem is exactly cld2. From that point on, the planner enters an infinite recursion cycle between indirect-transpose and indirect-before, eventually exhausting all available memory.

2.4. tensor_strides_decrease

The incorrect check in tensor_strides_decrease causes the 3x3 tile problem to enter indirect-before instead of indirect-after.

static int tensor_strides_decrease(const tensor *sz, inplace_kind k)
{
     if (FINITE_RNK(sz->rnk)) {
          int i;
          for (i = 0; i < sz->rnk; ++i)
               if ((sz->dims[i].os - sz->dims[i].is)
                   * (k == INPLACE_OS ? (INT)1 : (INT)-1) < 0)
                    return 1;
     }
     return 0;
}

tensor_strides_decrease() currently determines whether a stride decreases by comparing the signed stride values through (os - is). This works for positive strides, but it is not a valid ordering when negative strides are allowed, because a numerically smaller signed stride may have a larger magnitude. For example, with is = 2 and os = -8, the current code considers replacing 2 with -8 a decrease because -8 < 2, even though the stride magnitude actually increases from 2 to 8.

3. Proposed Fix

As a possible fix, I propose changing tensor_strides_decrease() to compare strides by magnitude rather than by their signed values. The current implementation effectively compares os < is for INPLACE_OS and is < os for INPLACE_IS, which does not provide the intended ordering when negative strides are involved. The revised comparison first compares X(iabs)(a) and X(iabs)(b), and uses the signed values only as a deterministic tie-breaker when the magnitudes are equal. tensor_strides_decrease() can then use stride_lt(os, is) for INPLACE_OS and stride_lt(is, os) for INPLACE_IS. This preserves the existing invariant that unequal strides have exactly one decreasing direction while making the ordering consistent with stride magnitude. For the problematic 3x3 tile, the new comparison correctly recognizes that |2| < |-8|, causing indirect-after to be accepted and indirect-before to be rejected. As a result, the indirect-before -> indirect-transpose recursion cycle is avoided.

static int stride_lt(INT a, INT b)
{
     INT aa = X(iabs)(a);
     INT ab = X(iabs)(b);

     return aa < ab || (aa == ab && a < b);
}

I believe this issue is worth addressing because certain valid layouts can trigger unbounded recursion in the planner, eventually exhausting all available memory. After the fix, the 3x3 tile problem is correctly accepted by indirect-after, which then successfully produces a valid plan. So I have submitted a pull request.#415

@matteo-frigo