Bound Larcform1 compile with folded flat DataLayout indexing [perf]

The unified DataLayouts back each layout with a 5-D parent array (and property
views with SubArrays of it). The finite-difference stencil getidx replicates a
point read 3^depth times; on the branch each read carried Base's full 5-D
index-decomposition (a plain Array's runtime col-major arithmetic, or a
property-view SubArray's linear-to-Cartesian reindex), so the deeply-inlined
stencil expression body exploded and Larcform1 compile-killed (OOM).

Add a constant-folded flat leaf-access path. `flat_mode(data)` is a pure type
function (Val(0/1/2), folds away) selecting:
  Val(1) IndexLinear parent  -> affine index into the parent (native arrayref).
  Val(2) shape-linear IndexCartesian property view (the column case) -> reduce
         to the root Array with a folded affine index (one hoistable
         first(indices) load), instantiating no Base reindex / div / rem.
  Val(0) everything else (dynamic Nh, box VIJFH view Nh>1) -> unchanged Cartesian
         get_struct fallback.
get_struct_linear/set_struct_linear! (+ a linear bitcast_struct_linear reusing
bitcast_struct_expr) address the Nf entries at base + (f-1)*fstride.

Measured (EDMF Larcform1 staged, interleaved): implicit_tendency! compile
32-33.5s vs 108-110s on the branch (~main's ~30-38s); peak RSS lowest of all.
Column FD-op runtime improves (e.g. GradientC2F 75.5->47.5 ns, div_interp
188->111 ns). Bit-identical values across scalar/multi-component, VIJFH/VIJHF,
Nh=1/Nh>1, property views, and setindex round-trip; device :kernel compile
unchanged (GPU-safe). Builds on the Slice-based property views (fast-linear
VIJHF) from the previous commit.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

Author: dennisYatunin

src/DataLayouts/bitcast_struct.jl

@@ -2,6 +2,14 @@
 field_expr(f, value_expr) = :(Core.getfield($value_expr, $f))
 field_expr(f, array_expr, index_expr) =
     :(@inbounds $array_expr[struct_index($f, $array_expr, $index_expr...)])
+# Flat/linear variant (see get_struct_linear): read field f of the value stored
+# at linear position `base` using a single linear load into the array's flat
+# storage. `base` and `fstride` are already fully constant-folded affine
+# expressions of the layout's static shape params, and `$(f - 1)` is a literal,
+# so `array[base + (f - 1) * fstride]` lowers to one arrayref with a precomputed
+# offset -- no Cartesian reindex, no div/rem, no per-dimension bounds math.
+field_expr(f, array_expr, base_expr, fstride_expr) =
+    :(@inbounds $array_expr[$base_expr + $(f - 1) * $fstride_expr])
 
 # Keep array element read instructions separate unless the full value is needed
 full_value_expr(@nospecialize(S), value_expr) = value_expr
@@ -155,3 +163,22 @@ For more information about `reinterpret` and padding, see the following:
     S = NTuple{num_indices, eltype(array)}
     return Expr(:block, :@inline, bitcast_struct_expr(T, S, :array, :index))
 end
+
+# Flat/linear analogue of the array method (see get_struct_linear): reads the
+# `num_indices` entries at linear positions `base, base + fstride, ...` using the
+# flat `field_expr(f, array, base, fstride)`. Reuses bitcast_struct_expr so the
+# type-reconstruction machinery is shared with the Cartesian path.
+@generated function bitcast_struct_linear(
+    ::Type{T},
+    array,
+    ::Val{num_indices},
+    base,
+    fstride,
+) where {T, num_indices}
+    S = NTuple{num_indices, eltype(array)}
+    return Expr(
+        :block,
+        :@inline,
+        bitcast_struct_expr(T, S, :array, :base, :fstride),
+    )
+end

src/DataLayouts/indexing.jl

@@ -70,17 +70,159 @@ is_invalid_linear(data, index) = index isa Integer && IndexStyle(data) isa Index
 # level) are identified by their length, since they keep their dimensions.
 is_trivial_point(data, index) = isone(length(data)) && index == CartesianIndex()
 
+# ---------------------------------------------------------------------------
+# Flat (constant-folded linear) leaf access.
+#
+# A point access reads `ncomponents(data)` entries out of the parent buffer. When
+# the parent has `IndexLinear` storage, `getindex(data, index)` used to route
+# through the Cartesian `get_struct`, which reads `parent[CartesianIndex(...)]`.
+# For a 5-D `Array` that generates the full 5-D linear-index arithmetic per tap,
+# and for a property-view `SubArray` it generates Base's linear-to-Cartesian
+# reindex -- the residual per-access bloat that inflates deeply-inlined stencil
+# expression bodies (and, for `SubArray`s, costs `div`/`rem` at runtime).
+#
+# `flat_mode(data)` picks one of three fully-static leaf-access strategies (the
+# `Val` is a pure function of the type, so the branch folds away and only the
+# selected `_{get,set}idx_flat` method is compiled per concrete layout):
+#
+#   Val(1)  parent is `IndexLinear` and every extent is inferable: index the
+#           parent directly with a constant-folded affine index. For a plain
+#           `Array` this is a native `arrayref`; for a fast property-view
+#           `SubArray` (e.g. `VIJHF`, F trailing) Base lowers `sub[i::Int]` to
+#           `parent(sub)[offset1 + i]` -- a single load into the root buffer.
+#   Val(2)  parent is an `IndexCartesian` property-view `SubArray` (e.g. a
+#           `VIJFH` field view, F not trailing) but the layout is shape-linear
+#           (all logical extents from the F axis onward are 1, the column case):
+#           reduce to the root `Array` directly with an affine index that is
+#           fully folded except for one hoistable `first(indices)` load, so no
+#           Base reindex / `div` / `rem` is instantiated.
+#   Val(0)  everything else (dynamic `Nh`, or a true box `VIJFH` field view with
+#           `Nh > 1`): the unchanged Cartesian path, which stays correct.
+@inline function flat_mode(data::DataLayout)
+    has_inferred_size(data) || return Val(0)
+    IndexStyle(data) isa IndexLinear && return Val(1)
+    F = f_dim(data)
+    F isa Integer || return Val(0)
+    # The root path omits every logical coordinate at or beyond the F axis, so it
+    # is only correct when those extents are all 1 (the column case). Unlike the
+    # layout IndexStyle, a `ncomponents <= 1` escape is NOT valid here: a scalar
+    # field view of a box (Nh > 1) still needs its h-term, which we drop.
+    all_ones(inferred_size(data)[F:end]...) &&
+        reducible_property_view(parent_type(data), Val(F)) && return Val(2)
+    return Val(0)
+end
+
+# A property-view SubArray (as built by struct_field_view) of an IndexLinear
+# array, whose indices are full `Slice`s except for a `UnitRange` at position F.
+# Reducing such a view to its root is a constant-stride affine map, so the flat
+# path can bypass Base's linear-to-Cartesian reindex.
+@inline reducible_property_view(::Type, ::Val) = false
+@inline reducible_property_view(
+    ::Type{<:SubArray{<:Any, N, PA, I}},
+    ::Val{F},
+) where {N, PA, I, F} =
+    F isa Integer && 1 <= F <= N && IndexStyle(PA) isa IndexLinear ?
+    _pv_indices_ok(I, Val(F), Val(N)) : false
+@generated function _pv_indices_ok(::Type{I}, ::Val{F}, ::Val{N}) where {I, F, N}
+    for k in 1:N
+        ft = fieldtype(I, k)
+        ok = k == F ? (ft <: AbstractUnitRange && !(ft <: Base.Slice)) : (ft <: Base.Slice)
+        ok || return :(false)
+    end
+    return :(true)
+end
+
+# Column-major linear offset (0-based) of coordinate tuple `idx` in an array of
+# static size `sz`. Tail recursion over tuples inlines to a folded Horner form;
+# `sz` entries are type-domain literals, only the coordinates are runtime values.
+@inline _col_major_offset(::Tuple{}, ::Tuple{}) = 0
+@inline _col_major_offset(sz::Tuple, idx::Tuple) =
+    (idx[1] - 1) + sz[1] * _col_major_offset(Base.tail(sz), Base.tail(idx))
+
+@inline _first_n(t::Tuple, ::Val{n}) where {n} = ntuple(k -> t[k], Val(n))
+
+# (base, fstride) for the Val(1) parent path. `base` is the linear position of
+# the first component of point `index`; `fstride` is the spacing between
+# successive components. For IndexLinear layouts the components form contiguous
+# length-`prod(inferred_size)` blocks, so `fstride == prod(inferred_size(data))`
+# and a linear point index maps to itself.
+@inline _flat_point_offset(data::DataLayout, index::Integer) =
+    (index, prod(inferred_size(data)))
+@inline _flat_point_offset(data::DataLayout, index::CartesianIndex) =
+    (_col_major_offset(inferred_size(data), Tuple(index)) + 1, prod(inferred_size(data)))
+
+# (root, base, fstride) for the Val(2) property-view root path. `fstart` (the
+# start of the field range) is the only runtime input; everything else folds
+# from the static shape. Shape-linearity guarantees every logical coordinate at
+# or beyond the F axis is 1, so only the `F - 1` leading dims contribute.
+@inline function _flat_root_offset(
+    data::DataLayout,
+    index::CartesianIndex,
+    fstart,
+    ::Val{F},
+) where {F}
+    before_sz = _first_n(inferred_size(data), Val(F - 1))
+    before_idx = _first_n(Tuple(index), Val(F - 1))
+    fstride = prod(before_sz)                       # root stride of the F axis
+    base = _col_major_offset(before_sz, before_idx) + (fstart - 1) * fstride + 1
+    return (base, fstride)
+end
+@inline _root_and_fstart(data::DataLayout, ::Val{F}) where {F} =
+    (parent(parent(data)), first(parentindices(parent(data))[F]))
+
+@propagate_inbounds _getidx_flat(::Val{0}, data, index) =
+    get_struct(parent(data), eltype(data), index, Val(f_dim(data)))
+@propagate_inbounds function _getidx_flat(::Val{1}, data, index)
+    (base, fstride) = _flat_point_offset(data, index)
+    return get_struct_linear(
+        parent(data),
+        eltype(data),
+        base,
+        fstride,
+        Val(ncomponents(data)),
+    )
+end
+@propagate_inbounds function _getidx_flat(::Val{2}, data, index)
+    (root, fstart) = _root_and_fstart(data, Val(f_dim(data)))
+    (base, fstride) = _flat_root_offset(data, index, fstart, Val(f_dim(data)))
+    return get_struct_linear(root, eltype(data), base, fstride, Val(ncomponents(data)))
+end
+
+@propagate_inbounds _setidx_flat!(::Val{0}, data, value, index) =
+    set_struct!(parent(data), convert(eltype(data), value), index, Val(f_dim(data)))
+@propagate_inbounds function _setidx_flat!(::Val{1}, data, value, index)
+    (base, fstride) = _flat_point_offset(data, index)
+    return set_struct_linear!(
+        parent(data),
+        convert(eltype(data), value),
+        base,
+        fstride,
+        Val(ncomponents(data)),
+    )
+end
+@propagate_inbounds function _setidx_flat!(::Val{2}, data, value, index)
+    (root, fstart) = _root_and_fstart(data, Val(f_dim(data)))
+    (base, fstride) = _flat_root_offset(data, index, fstart, Val(f_dim(data)))
+    return set_struct_linear!(
+        root,
+        convert(eltype(data), value),
+        base,
+        fstride,
+        Val(ncomponents(data)),
+    )
+end
+
 # Always convert to the element type of a DataLayout when modifying its values.
 # Represent every single-point DataLayout view using a zero-dimensional DataF.
 @propagate_inbounds Base.setindex!(data::DataLayout, value, index::PointIndex) =
     is_invalid_linear(data, index) ? setindex!(data, value, CartesianIndices(data)[index]) :
     is_trivial_point(data, index) ?
     set_struct!(parent(data), convert(eltype(data), value)) :
-    set_struct!(parent(data), convert(eltype(data), value), index, Val(f_dim(data)))
+    _setidx_flat!(flat_mode(data), data, value, index)
 @propagate_inbounds Base.getindex(data::DataLayout, index::PointIndex) =
     is_invalid_linear(data, index) ? getindex(data, CartesianIndices(data)[index]) :
     is_trivial_point(data, index) ? get_struct(parent(data), eltype(data)) :
-    get_struct(parent(data), eltype(data), index, Val(f_dim(data)))
+    _getidx_flat(flat_mode(data), data, index)
 @propagate_inbounds Base.view(data::DataLayout, index::PointIndex) =
     is_invalid_linear(data, index) ? view(data, CartesianIndices(data)[index]) :
     is_trivial_point(data, index) ? data :

src/DataLayouts/struct_storage.jl

@@ -227,6 +227,54 @@ julia> get_struct([0 2; 0 0; 0 1; 0 0], Tuple{Int32, Int32, Int128}, 5, Val(1))
     return bitcast_struct(T, array, Val(Nf), index...)
 end
 
+"""
+    get_struct_linear(array, T, base, fstride, Val(Nf))
+    set_struct_linear!(array, value, base, fstride, Val(Nf))
+
+Flat-storage analogues of [`get_struct`](@ref) and [`set_struct!`](@ref) that
+address the `Nf` entries of a single value with a fully affine linear index into
+`array`'s storage: entry `f` lives at `base + (f - 1) * fstride`. Both `base` and
+`fstride` are expected to be constant-folded functions of a [`DataLayout`](@ref)'s
+static shape parameters (see `_flat_point_offset` in `indexing.jl`).
+
+The key property is that `array[base + (f - 1) * fstride]` is a *linear* index:
+- for a plain `Array` it lowers to a native `arrayref` with a precomputed offset;
+- for a fast (`IndexLinear`) property-view `SubArray` Base lowers it to
+  `parent[offset1 + stride1 * i]`, i.e. a single load into the root buffer with
+  no Cartesian reindex, no `div`/`rem`, and no `SignedMultiplicativeInverse`.
+
+The caller must only use these when `array` has `IndexLinear` storage; otherwise
+the linear index would trigger Base's expensive linear-to-Cartesian conversion.
+"""
+@inline function get_struct_linear(
+    array,
+    ::Type{T},
+    base,
+    fstride,
+    ::Val{Nf},
+) where {T, Nf}
+    # Bounds-check only the two endpoints of the arithmetic progression, which
+    # bound every intermediate index; constructing a StepRange here would emit a
+    # checked-remainder + overflow chain that dwarfs the actual read.
+    @boundscheck (checkbounds(array, base); checkbounds(array, base + (Nf - 1) * fstride))
+    return bitcast_struct_linear(T, array, Val(Nf), base, fstride)
+end
+
+@inline function set_struct_linear!(
+    array,
+    value::T,
+    base,
+    fstride,
+    ::Val{Nf},
+) where {T, Nf}
+    @boundscheck (checkbounds(array, base); checkbounds(array, base + (Nf - 1) * fstride))
+    entries = bitcast_struct(NTuple{Nf, eltype(array)}, value)
+    unrolled_foreach(enumerate(entries)) do (i, entry)
+        @inbounds array[base + (i - 1) * fstride] = entry
+    end
+    return array
+end
+
 # Indices for a view of one value's entries, chosen so the view never has more than one
 # dimension: Cartesian indices split into scalar components plus a range along the F axis.
 # A Cartesian range would instead build a much costlier view with all singleton dimensions.