In-place ops and caller-provided destinations

benchmark_report.md

@@ -494,6 +494,76 @@ the rank-1 layout, GPU `add` at 16 M now runs at **236.4 GB/s against PyTorch's
 232.4** — 1.02× faster, where the previous report had it 1.47× slower. `relu`
 went from 1.9× slower to parity, and `x*s+t` from parity to 1.97× faster.
 
+## 9. Every op allocated its own result — in-place and `out=` forms added
+
+Until now every `Tensor` method was `const` and returned a freshly allocated
+tensor. In a training loop that is one allocation *and* one free per op per
+iteration, and a peak footprint equal to the whole expression graph rather than
+the working set. PyTorch exposes `add_`/`mul_` and `out=` for exactly this.
+
+The launchers already took a separate `dst`, so the work was on the public
+surface: `add_out(rhs, out)` is now the single implementation of each op, the
+allocating form builds `out` and calls it, and `add_(rhs)` calls it with
+`*this` as `out`. See CLAUDE.md for the surface and the aliasing rules.
+
+**One op, measured 2026-09-04, GB10, Release** (min of 7 batches, through the
+Python bindings so allocation and synchronization are inside the timing —
+`scripts/bench_inplace.py`):
+
+| device | n | op | allocating | in-place | `out=` | in-place speedup |
+|---|---:|---|---:|---:|---:|---:|
+| cpu | 1 K | add | 1.79 µs | 1.22 µs | 1.38 µs | 1.47× |
+| cpu | 1 K | relu | 1.31 µs | 0.80 µs | 0.99 µs | 1.65× |
+| cpu | 64 K | add | 6.34 µs | 5.70 µs | 5.81 µs | 1.11× |
+| cpu | 16 M | add | 2029 µs | 1661 µs | 2015 µs | 1.22× |
+| cpu | 16 M | mul (scalar) | 1413 µs | 972 µs | 1340 µs | 1.45× |
+| cuda | 1 K | add | 7.93 µs | 6.48 µs | 6.65 µs | 1.22× |
+| cuda | 64 K | add | 8.21 µs | 6.75 µs | 6.92 µs | 1.22× |
+| cuda | 1 M | add | 14.05 µs | 12.60 µs | 12.72 µs | 1.12× |
+| cuda | 16 M | add | 862 µs | 827 µs | 823 µs | 1.04× |
+
+**A 32-step chain** (`w += g`, the shape a training step actually has):
+
+| device | n | allocating | in-place | speedup |
+|---|---:|---:|---:|---:|
+| cpu | 1 K | 61.4 µs | 39.4 µs | 1.56× |
+| cpu | 64 K | 206 µs | 187 µs | 1.10× |
+| cpu | 16 M | 63.1 ms | 54.6 ms | 1.16× |
+| cuda | 1 K | 264 µs | 212 µs | 1.25× |
+| cuda | 64 K | 273 µs | 221 µs | 1.24× |
+| cuda | 1 M | 489 µs | 412 µs | 1.19× |
+| cuda | 16 M | 30.8 ms | 26.6 ms | 1.16× |
+
+Read the shape of it, not the individual rows. The saving is largest where
+allocation is a large fraction of the op (1.2–1.6× below ~64 K on both
+backends, and it is the same ~1.4 µs per op on CPU and ~1.5 µs on CUDA that §4
+already identified as fixed overhead), and it shrinks toward 1.04× at 16 M
+where the kernel itself dominates and `HostPool` / `cudaMallocAsync` are
+already recycling the block. That last figure is the honest ceiling for the
+*time* argument.
+
+The memory argument does not shrink: an in-place chain holds one buffer whatever
+its length, and the test suite asserts exactly that (`data_ptr` is unchanged
+across a 100-step loop, in `tests/test_inplace.cpp` and
+`python/tests/test_inplace.py`) — a correct-but-reallocating implementation
+would pass every value check and fail those.
+
+**A trap this script already sprang.** An earlier draft reused one destination
+tensor across the three ops at each size and reported CUDA `relu_` at 0.38× —
+a 2.6× *regression* that does not exist. Re-measured with per-case operands it
+is 1.05×. Every case now builds its own; and the CPU rows at 1 M were dropped
+from the table because they swing between 11 µs and 250 µs run to run on a
+loaded machine, which is contention, not the op.
+
+The change also removed `__restrict__` from the three kernels in
+`contiguous.cuh`: in-place calls them with `dst` equal to a source, which is
+exactly the aliasing the qualifier promises does not happen. It cost nothing —
+`add` at 16 M went from 230.2–233.4 GB/s to 228.6–230.7 and `relu` from
+232.9–234.6 to 233.4–235.5, i.e. inside run-to-run spread. The read-only path
+comes from the explicit `__ldg` in `device_load`, not from restrict-driven
+inference, and each pointer is touched once per thread, so there was nothing
+for the qualifier to buy.
+
 ## Full results
 
 ### CPU, default threading (20 OpenMP threads)
@@ -668,6 +738,10 @@ OPENMAT_LIB=build-release/OpenMat.so PYTHONPATH=python \
 OPENMAT_LIB=build-release/OpenMat.so PYTHONPATH=python \
   bench-env/bin/python scripts/bench_rank_sweep.py
 
+# §9 — in-place / out= against the allocating forms
+OPENMAT_LIB=build-release/OpenMat.so PYTHONPATH=python \
+  bench-env/bin/python scripts/bench_inplace.py
+
 # §7 — OpenMP scaling, one invocation per thread count
 for t in 1 20; do
   OMP_NUM_THREADS=$t OPENMAT_LIB=build-release/OpenMat.so PYTHONPATH=python \

compile_commands.json

@@ -185,6 +185,12 @@
   "file": "/home/apalese/OpenMat/tests/test_contiguous.cpp",
   "output": "tests/CMakeFiles/test_contiguous.dir/test_contiguous.cpp.o"
 },
+{
+  "directory": "/home/apalese/OpenMat/build/tests",
+  "command": "/usr/bin/c++  -I/home/apalese/OpenMat/headers -I/usr/local/cuda/targets/sbsa-linux/include -isystem /home/apalese/OpenMat/build/_deps/googletest-src/googletest/include -isystem /home/apalese/OpenMat/build/_deps/googletest-src/googletest -g -fopenmp -o CMakeFiles/test_inplace.dir/test_inplace.cpp.o -c /home/apalese/OpenMat/tests/test_inplace.cpp",
+  "file": "/home/apalese/OpenMat/tests/test_inplace.cpp",
+  "output": "tests/CMakeFiles/test_inplace.dir/test_inplace.cpp.o"
+},
 {
   "directory": "/home/apalese/OpenMat/build/_deps/googletest-build/googlemock",
   "command": "/usr/bin/c++  -I/home/apalese/OpenMat/build/_deps/googletest-src/googlemock/include -I/home/apalese/OpenMat/build/_deps/googletest-src/googlemock -isystem /home/apalese/OpenMat/build/_deps/googletest-src/googletest/include -isystem /home/apalese/OpenMat/build/_deps/googletest-src/googletest -g -std=c++17 -Wall -Wshadow -Wno-error=dangling-else -DGTEST_HAS_PTHREAD=1 -fexceptions -Wextra -Wno-unused-parameter -Wno-missing-field-initializers -o CMakeFiles/gmock.dir/src/gmock-all.cc.o -c /home/apalese/OpenMat/build/_deps/googletest-src/googlemock/src/gmock-all.cc",

docs/roadmap.md

@@ -4,7 +4,7 @@ Questo documento elenca le funzionalità del framework, il loro stato attuale e
 
 **Legenda:** ✅ implementato · ⚠️ implementato parzialmente o con divergenze rispetto al progetto originale · ❌ da fare
 
-Ultimo allineamento al codice: 4 settembre 2026 (benchmark rieseguiti su GB10, Release; vedi [`benchmark_report.md`](../benchmark_report.md)).
+Ultimo allineamento al codice: 4 settembre 2026 (in-place / `out=` aggiunti, §8; benchmark rieseguiti su GB10, Release; vedi [`benchmark_report.md`](../benchmark_report.md)).
 
 ---
 
@@ -205,11 +205,9 @@ Il bug latente è risolto. `apply` e `apply_binary` in [`headers/tensor.inl`](..
 
 Entrambi passano per `float` internamente (`static_cast<float>(x) > 0.0f`, `expf(-static_cast<float>(x))`): è ciò che li fa funzionare senza modifiche su `float16_t`.
 
-### 4.4 Overload sullo stream per `apply_binary` ❌
+### 4.4 Overload sullo stream per `apply_binary` ✅
 
-Asimmetria residua: `apply` ha l'overload `(Op, const Stream&)`, `apply_binary` no. Di conseguenza `scale_shift`, `shift_scale` e i quattro `fused_*_*` sono chiamate sincrone sullo stream di default, mentre `relu` e `sigmoid` no.
-
-Il lavoro è meccanico e segue `apply`: accettare `const Stream& s`, costruire l'output con il costruttore privato `Tensor(shape, device, Stream)` e passare `s.get()` al launcher — `launch_apply_binary_op` accetta già un `cudaStream_t`.
+Risolto come effetto collaterale di §8: `apply_binary` è stato riscritto sopra `apply_binary_out(rhs, op, out, s)`, che uno stream lo prende per costruzione, e le forme allocante e in-place ci delegano. `scale_shift`, `shift_scale` e i quattro `fused_*_*` restano sincroni sullo stream di default perché non hanno un overload proprio — chiamarli su uno stream significa passare per `apply`/`apply_binary` con lo stesso functor.
 
 ---
 
@@ -347,6 +345,34 @@ Suite di test: `python/tests/test_tensor.py`, `test_tensor_api.py`, `test_dtypes
 
 ---
 
+## 8. Operazioni in-place e destinazione fornita dal chiamante ✅
+
+Ogni metodo di `Tensor` era `const` e allocava il risultato. In un ciclo di training questo significa un'allocazione **e** una deallocazione per operazione per iterazione, e un picco di memoria pari a tutto il grafo dell'espressione invece che al working set. PyTorch espone `add_`, `mul_` e `out=` esattamente per questo.
+
+I launcher accettavano già una `dst` separata, quindi il lavoro era tutto sulla superficie pubblica. Ogni operazione ha ora tre forme, di cui **una sola è un'implementazione vera**:
+
+```cpp
+auto c = a.add(b, s);        // alloca il risultato, poi chiama add_out
+a.add_out(b, out, s);        // ← il corpo: scrive in una destinazione che esiste già
+a.add_(b, s);                // == a.add_out(b, a, s)
+```
+
+`add_out` restituisce `Tensor&` (la destinazione) e `add_` restituisce `*this`, quindi le chiamate si concatenano; le forme senza stream sono delegate di una riga con `Stream::default_stream()`. Ci sono anche `operator+=`/`-=`/`*=`/`/=`.
+
+**Famiglie coperte:** `add`/`sub`/`mul`/`div` (tensore e scalare), `apply`, `apply_binary`, `relu`, `sigmoid`, `fill_`. `matmul`, `transpose` e `permute` hanno solo `_out`: i loro kernel leggono elementi che non scrivono, quindi una destinazione che condivide il buffer con un operando leggerebbe valori già sovrascritti — `_check_alias_none` lo rifiuta al call site invece di restituire un risultato plausibile e sbagliato.
+
+**Perché la famiglia elementwise *può* fare aliasing:** il loop CPU, il fast path contiguo GPU e i kernel per rank leggono l'indice i e scrivono l'indice i, quindi `dst == lhs` è corretto esattamente quanto un buffer separato. Questo dipende dal fatto che ogni buffer sia una singola sequenza piatta: `_check_alias_elementwise` ricontrolla `is_contiguous()`, così il giorno in cui una view con stride potrà puntare a una regione *diversa* della stessa allocazione (§2.1b) l'operazione solleva un'eccezione invece di calcolare la cosa sbagliata.
+
+Conseguenza collaterale: i tre kernel di [`contiguous.cuh`](../headers/ops/kernels/contiguous.cuh) hanno perso `__restrict__`, che è precisamente la promessa che l'in-place viola. Misurato, non è costato nulla (`add` a 16M da 230.2-233.4 a 228.6-230.7 GB/s, `relu` da 232.9-234.6 a 233.4-235.5): il percorso read-only viene dall'`__ldg` esplicito in `device_load`, non da inferenza sul restrict.
+
+Il `_out` è anche il punto in cui vive ora la validazione degli operandi (shape e device controllati prima del dispatch), quindi un `add` GPU fra shape diverse solleva un'eccezione invece di leggere oltre la fine di un buffer come faceva prima.
+
+**Numeri** ([`benchmark_report.md` §9](../benchmark_report.md#9-every-op-allocated-its-own-result--in-place-and-out-forms-added), harness `scripts/bench_inplace.py`): 1.2-1.6× sotto i ~64 K elementi su entrambi i backend, dove l'allocazione è una frazione grande dell'operazione, che si riduce a 1.04× a 16 M dove domina il kernel e `HostPool`/`cudaMallocAsync` stanno già riciclando il blocco. L'argomento sulla memoria invece non si riduce: `tests/test_inplace.cpp` e `python/tests/test_inplace.py` verificano che `data_ptr` non cambi lungo un ciclo di 100 passi — la proprietà che un'implementazione corretta ma ri-allocante fallirebbe superando ogni controllo sui valori.
+
+**Non fatto:** `matmul` batch/in-place resta fuori (§5.1), e i quattro `fused_*_*` non hanno una forma in-place dedicata — `apply_binary_(rhs, BinaryCompose<…>{…})` con lo stesso functor è il modo di ottenerla.
+
+---
+
 ## Priorità suggerite
 
 | Done | Priorità | Item |
@@ -367,7 +393,8 @@ Suite di test: `python/tests/test_tensor.py`, `test_tensor_api.py`, `test_dtypes
 | ✅ | — | 6.6 — `matmul_cpu` `ikj` + tiling L2 + OpenMP: 1.81 → 123 GFLOP/s a 1024³ (§5) |
 | ✅ | — | 6.7 — `reduce_sum_cpu` a 8 lane: 7.7 → 36.9 GB/s, ora davanti a NumPy (§1.2) |
 | ✅ | — | 6.8 — `HostPool`/`PinnedHostPool`: round-trip 100 × 64 MB da 4.3 a 51.1 GB/s |
-| | **Alta** | 4.4 — Overload stream per `apply_binary` (chiude l'asimmetria di 6.1) |
+| ✅ | — | 8 — Operazioni in-place (`add_`, `mul_`, `relu_`, `fill_`) e overload `_out` con destinazione fornita dal chiamante |
+| ✅ | — | 4.4 — Overload stream per `apply_binary` (arrivato con §8: `apply_binary_out` prende lo stream per costruzione) |
 | | **Alta** | 1.1 — Unarie `abs`, `sqrt`, `exp`, `log` via `apply` (costo basso, alto valore) |
 | | **Alta** | 6.4b — OpenMP su `apply`/`apply_binary`: stesso `if(_total > 65536)` di 6.4, ~2× atteso a 16M su tutta la famiglia fused |
 | | Media | 6.4c — `min`/`max` con `parallel for reduction(min:)/(max:)` fra thread (2.7× dietro NumPy); **non** la clausola `simd`, già misurata e rimossa |

headers/ops/kernels/contiguous.cuh

@@ -47,6 +47,16 @@
 // A grid-stride loop was also tried and is not used: capping the grid to a few
 // waves per SM costs 8-12 %, and at the exact block count the loop body runs
 // once and buys nothing over a plain bounds check.
+//
+// None of the three kernels marks its pointers __restrict__, and that is
+// deliberate: the in-place family (Tensor::add_, relu_, …) calls them with dst
+// equal to one of the sources, which is exactly the aliasing __restrict__
+// promises does not happen. The promise bought nothing here — each pointer is
+// read or written once per thread, and the read-only path comes from the
+// explicit __ldg in device_load, not from restrict-driven inference. Measured
+// on the reference GB10 at 16 M elements, removing it moved `add` from
+// 230.2-233.4 GB/s to 228.6-230.7 and `relu` from 232.9-234.6 to 233.4-235.5,
+// i.e. nothing outside run-to-run spread. Re-measure before putting it back.
 // ─────────────────────────────────────────────────────────────────────────────
 
 namespace om {
@@ -83,8 +93,8 @@ namespace detail {
     // ill-formed and a member-wise copy is free to lower to two 2-byte accesses,
     // which is the whole thing this path exists to avoid.
     template <typename T, int V, typename Op>
-    __global__ void contig_binary_kernel(const T* __restrict__ a, const T* __restrict__ b,
-                                         T* __restrict__ dst, size_t n_pack, size_t n, Op op)
+    __global__ void contig_binary_kernel(const T* a, const T* b,
+                                         T* dst, size_t n_pack, size_t n, Op op)
     {
         const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
 
@@ -114,7 +124,7 @@ namespace detail {
     }
 
     template <typename T, int V, typename Op>
-    __global__ void contig_unary_kernel(const T* __restrict__ src, T* __restrict__ dst,
+    __global__ void contig_unary_kernel(const T* src, T* dst,
                                         size_t n_pack, size_t n, Op op)
     {
         const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
@@ -141,7 +151,7 @@ namespace detail {
     }
 
     template <typename T, int V>
-    __global__ void contig_fill_kernel(T* __restrict__ dst, T value, size_t n_pack, size_t n)
+    __global__ void contig_fill_kernel(T* dst, T value, size_t n_pack, size_t n)
     {
         const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;