Fortran 2008 strided subarrays are passed with the wrong layout

[eschnett]

While reviewing mpif I found an error in MPICH's implementation of strided arrays in the Fortran 2008 bindings. MPICH creates a datatype to describe the strided array layout, and this code contains an error. I have asked AI to describe the error, and I attach this report below.

This is for MPICH's commit ab53493dad85ffee0fc95812b250e1c8dacf7982 which is at or near the head of the main branch.

Summary

With MPI_SUBARRAYS_SUPPORTED = .true., passing an array section whose inner
dimension is strided and which has at least one dimension above it silently
transfers the wrong data. On integer :: a(20,3),

call MPI_Sendrecv(a(1:20:2, :), 30, MPI_INTEGER, 0, 0, &
                  flat, 30, MPI_INTEGER, 0, 0, MPI_COMM_SELF, MPI_STATUS_IGNORE)

delivers 20 of the 30 elements wrong. No error is returned. Every routine whose
cdesc wrapper walks a descriptor is affected — point-to-point, RMA, file I/O,
MPI_Bcast, the fixed sides of the collectives.

Where

src/binding/fortran/use_mpi_f08/wrappers_c/cdesc.c, in
cdesc_create_datatype:

    int accum_sm = cdesc->elem_len;                      /* line 19 */
    ...
        if (extent <= INT_MAX) {
            if (cdesc->dim[i].sm == accum_sm) {          /* line 59 */
                mpi_errno = PMPI_Type_contiguous(extent, types[i], &types[i + 1]);
            } else {
                mpi_errno = PMPI_Type_create_hvector(extent, 1, cdesc->dim[i].sm,
                                                     types[i], &types[i + 1]);
            }
        } else {
            ...                                          /* line 66, same test */
        }
        ...
        accum_sm = cdesc->dim[i].sm * cdesc->dim[i].extent;   /* line 84 */

accum_sm is the span the levels below cover: stride times extent of the
outermost of them. But MPI_Type_contiguous replicates at multiples of the
inner type's extent, and the two are not the same number once any lower
level was built with MPI_Type_create_hvector: an hvector of n blocks at
stride sm over a type of extent e has extent

sm*(n-1) + e     not     sm*n

So dim[i].sm == accum_sm is satisfied — the dimension really is dense against
the span below it — while PMPI_Type_contiguous then places replica i at
i * (accum_sm - (sm - e)). Every replica but the first is short, by a
multiple of the deficit. Both the extent <= INT_MAX branch and the _c
branch have the same test.

For a(1:20:2, :): dimension 0 has sm = 8, accum_sm = elem_len = 4, so it
becomes hvector(10, 1, 8, MPI_INT), of extent 76. accum_sm then becomes
8 * 10 = 80. Dimension 1 has sm = 80, which equals accum_sm, so it becomes
contiguous(3, that) — and columns 2 and 3 land at bytes 76 and 152 instead of
80 and 160, i.e. 1 and 2 elements early.

Verifiable without any Fortran:

MPI_Datatype hv;
MPI_Type_create_hvector(10, 1, 8, MPI_INT, &hv);
MPI_Aint lb, ext;
MPI_Type_get_extent(hv, &lb, &ext);   /* 76; the span it covers is 80 */

Measured on the build below:

extent(hvector(10,1,8,MPI_INT)) = 76   span it covers = 80
extent(contiguous(3, that))     = 228   span needed  = 240

Reproducer

program mpich_cdesc_bug
  use mpi_f08
  implicit none
  integer, parameter :: n = 20, m = 3
  integer :: a(n, m), flat(n / 2 * m)
  integer :: i, j, p, nwrong

  call MPI_Init()
  print '(a,l1)', 'MPI_SUBARRAYS_SUPPORTED = ', MPI_SUBARRAYS_SUPPORTED

  do j = 1, m
     do i = 1, n
        a(i, j) = 100 * j + i
     end do
  end do
  flat = -1
  call MPI_Sendrecv(a(1:n:2, :), 30, MPI_INTEGER, 0, 0, &
                    flat, 30, MPI_INTEGER, 0, 0, &
                    MPI_COMM_SELF, MPI_STATUS_IGNORE)

  nwrong = 0
  p = 0
  do j = 1, m
     do i = 1, n, 2
        p = p + 1
        if (flat(p) /= a(i, j)) nwrong = nwrong + 1
     end do
  end do
  print '(a,i0,a)', 'wrong elements: ', nwrong, ' of 30'
  print '(a)', 'expected:'
  print '(10i5)', ((a(i, j), i = 1, n, 2), j = 1, m)
  print '(a)', 'received:'
  print '(10i5)', flat
  call MPI_Finalize()
end program mpich_cdesc_bug

Output, one rank:

MPI_SUBARRAYS_SUPPORTED = T
wrong elements: 20 of 30
expected:
  101  103  105  107  109  111  113  115  117  119
  201  203  205  207  209  211  213  215  217  219
  301  303  305  307  309  311  313  315  317  319
received:
  101  103  105  107  109  111  113  115  117  119
  120  202  204  206  208  210  212  214  216  218
  219  301  303  305  307  309  311  313  315  317

a(20,1) = 120 arrives where a(1,2) = 201 belongs — the 4-byte shift — and
column 3 is 8 bytes early.

The rule is general: any dimension above a strided one. Rank 3 with the stride
in the middle, integer :: t(2,4,3) and the section t(:, 1:4:2, :) (dense
innermost, strided, then dense) gives 8 of 12 wrong:

expected:
  111  112  131  132  211  212  231  232  311  312  331  332
received:
  111  112  131  132  141  142  221  222  231  232  311  312

A receiving section is corrupted the same way, and the wrong type also writes
outside the section the caller named. Receiving 30 elements into
b(1:20:2, :), with b pre-filled with -7:

section elements wrong: 20 of 30
elements clobbered outside the section: 10 of 30
b(20,1) = 1011 (should still be -7)

The displacements only shrink, so nothing lands past the end of b in this
case — but the even-numbered rows of b, which the section does not name, are
overwritten.

Environment

• MPICH main at ab53493dad85ffee0fc95812b250e1c8dacf7982, reporting
  MPICH Version: 5.1.0a1
• ./configure --disable-dependency-tracking --disable-doc --enable-cxx=no --enable-fortran --enable-mpi-abi --enable-shared=yes --enable-static=no --with-device=ch3 --with-hwloc
• macOS 26.6 (Darwin 25.6.0), arm64; gfortran 15 (MacPorts), clang for C
• MPI_SUBARRAYS_SUPPORTED is .true. in this build

The walker is independent of --enable-mpi-abi; nothing in the affected path
consults it.

Suggested fix

Either would do:

1. Take the contiguous branch only while the invariant
   extent(types[i]) == accum_sm still holds — that is, only while every level
   built so far has been dense. Once one level is an hvector, keep using
   hvector for the levels above it: an hvector places its blocks at exact byte
   offsets, so it is correct whatever the inner extent is, and it costs one
   datatype either way.
2. Or resize each hvector level to the span it covers
   (MPI_Type_create_resized(..., 0, accum_sm, ...)), restoring the invariant
   so that contiguous stays usable above it. This costs an extra datatype per
   strided level and gives the type lb/ub markers.

We took (1) in a downstream Fortran binding, seeding the flag by measuring the
base type's extent rather than assuming extent(oldtype) == elem_len (a
resized oldtype would break the innermost level the same way). The one-line
version of the change is to gate line 59 and line 66 on a dense flag and
clear that flag whenever the hvector branch is taken.

Separate observation, not measured

In the same function accum_sm and extent are int (lines 19 and 22) while
cdesc->dim[i].sm and cdesc->dim[i].extent are CFI_index_t. The
extent <= INT_MAX guard at line 58 covers the argument passed to the
constructors but not the accum_sm = cdesc->dim[i].sm * cdesc->dim[i].extent
product at line 84, which is computed in int. A section spanning more than
2 GiB would overflow it and then compare wrongly at line 59. I have not tried
to trigger this.