mpi: fix local_size_many_transposed when yblock == n[1] (#174)

[meng004]

The transposed output distribution was gated by a strict yblock < n[1], so passing yblock == n[1] fell through to the non-transposed fallback and left dims[1].ob = n[1]. That made rank 1 non-idle and report local_ny = n1 instead of 0, disagreeing with fftw_mpi_plan_many_transpose which distributes n1 by yblock unconditionally. Use <= so the boundary case takes the transposed distribution the planner already produces.

Problem

fftwf_mpi_local_size_many_transposed returns an inconsistent partition of the
transposed (y) dimension when the explicitly-specified yblock equals n[1]. As
reported in #174, for grid n = {6, 5, 4} with block0 = block1 = 5 on 2 MPI
ranks, the transposed y-domain local size is reported as 5 on both ranks
(sum = 10 ≠ n1 = 5). Rank 1's y-size should be 0. The reporter also observes
that swapping the dimension order (so n1 lands on the input side instead of
the output side) produces the correct result — pinpointing the output-dimension
sizing path.

Root cause

mpi/api.c, XM(local_size_many_transposed) (local_size_many_transposed):

/* default 1d block distribution, with transposed output
   if yblock < n[1] */
dims[0].ib = xblock;
if (rnk > 1) {
     if (yblock < n[1])          /* strict inequality */
          dims[1].ob = yblock;   /* transposed: distribute dim 1 by yblock */
     else
          dims[0].ob = xblock;   /* dims[1].ob stays = n[1] from simple_dims() */
}

simple_dims() initializes dims[i].ob = n[i]. When the caller passes
yblock == n[1] (5), the guard 5 < 5 is false, so the else branch runs,
dims[1].ob is never set to yblock and remains n[1] — i.e. "dimension 1 not
distributed, full extent on every non-idle rank". The subsequent
local_size_guru then keeps dim 0 distributed on output while treating dim 1 as
full, so rank 1 is non-idle and reports local_ny = 5 instead of 0.

This is a genuine defect, not a documented boundary behavior, because the sizing
function disagrees with the actual transpose planner for the same arguments:
fftwf_mpi_plan_many_transpose(n0, n1, howmany, block0, block1, …) distributes
n1 by yblock unconditionally (num_blocks(5, 5) = 1 block ⇒ rank 0 gets all
5 y-rows, rank 1 gets 0). A user who sizes allocations from local_size and
then plans gets a mismatch. The reporter's dims-swap observation maps exactly to
this guard: the input-side assignment dims[0].ib = xblock is unconditional and
yields the correct size-0 idle rank when block == n, while the output-side
dim-1 assignment is gated by the strict < — so n1 is handled correctly on
the input side and incorrectly on the output side.

Fix

Change the strict < to <= so yblock == n[1] also takes the transposed
distribution (dims[1].ob = yblock), matching what the planner does:

/* default 1d block distribution, with transposed output
   if yblock <= n[1] */
dims[0].ib = xblock;
if (rnk > 1) {
     if (yblock <= n[1])
          dims[1].ob = yblock;
     else
          dims[0].ob = xblock;
}

Minimal and safe by construction:

• Default path unchanged. For FFTW_MPI_DEFAULT_BLOCK (yblock = 0), both
  0 < n1 and 0 <= n1 are true — same branch, default_sz() fills the block
  as before.
• yblock < n1 path unchanged — both old and new take the if branch.
• Only the single value yblock == n1 changes: from a silent non-transposed
  fallback (wrong sizing) to the correct 1-block transposed distribution that
  the planner already uses.

Test evidence

Validated on FFTW 3.3.10 (single precision + MPI, --enable-mpi --enable-float,
MPICH, mpiexec -n 2) with two arms: pristine vs. patched.

check	before (pristine)	after (patched)
oracle y-partition, n1=5: sum(local_ny) disjoint	r0 ny=5 / r1 ny=5 → sum 10 ≠ 5	r0 ny=5 / r1 ny=0 → sum 5, valid partition
control: non-transposed x-domain local_size_many (n0=6)	r0 nx=5 / r1 nx=1 (sum 6)	byte-identical
control: square {5,5,4}, default blocks	r0 (3,3) / r1 (2,2)	byte-identical
round-trip: alloc via local_size_many_transposed, plan many_transpose, execute, compare to analytic transpose	FAIL — rank 1 believes it owns 5 output rows the plan never wrote, reads uninitialized data (maxerr ≈ 7e18)	PASS — rank 0 checks 5 rows maxerr=0, rank 1 owns 0 rows, clean

Both controls are byte-identical across arms (verified by diff), so the change
raises no false alarm on the non-transposed path or on a symmetric grid. The
round-trip on the fixed arm proves the sizes are usable end-to-end; the
round-trip on the buggy arm proves the defect is harmful (a correct consumer
reads data the planner never produced), not a benign over-estimate.

How to reproduce

#include <fftw3-mpi.h>
#include <stdio.h>
int main(int argc, char **argv) {
    ptrdiff_t n[3] = {6, 5, 4};       /* {zn, yn, xn} */
    MPI_Init(&argc, &argv);
    fftwf_mpi_init();
    int rank; MPI_Comm_rank(MPI_COMM_WORLD, &rank);   /* run with -n 2 */
    ptrdiff_t lnx, lxs, lny, lys;
    fftwf_mpi_local_size_many_transposed(
        3, n, /*howmany=*/1, /*block0=*/5, /*block1=*/5, MPI_COMM_WORLD,
        &lnx, &lxs, &lny, &lys);
    printf("rank %d: local_ny=%td y_start=%td\n", rank, lny, lys);
    /* before the fix both ranks print local_ny=5 (sum 10, should be 5) */
    fftwf_mpi_cleanup(); MPI_Finalize(); return 0;
}

Build with mpicc repro.c -lfftw3f_mpi -lfftw3f and run mpiexec -n 2 ./a.out.
Before the fix both ranks report local_ny=5; after, rank 0 reports 5 and rank 1
reports 0.

---