diff --git a/include/fastfields/core/autocast.h b/include/fastfields/core/autocast.h
index 6886480..cc1238a 100644
--- a/include/fastfields/core/autocast.h
+++ b/include/fastfields/core/autocast.h
@@ -1,6 +1,12 @@
 #pragma once
 #include <cstddef>
 #include <cstdint>
+// `hostNew` / `hostDelete` throw std::runtime_error on the pinned-host
+// (cudaMallocHost) path. This header never included <stdexcept> for them: it
+// compiled only because every translation unit that includes it happens to
+// include <stdexcept> first. A standalone .cu that includes just this header
+// fails to compile without this line.
+#include <stdexcept>
 #include <fastfields/core/dlpack.h>
 #include <fastfields/core/cuda_switch.h>
 
@@ -195,5 +201,114 @@ inline void free_if_needed(OutPointer ptr)
     _copy_if_needed<OutPointer, InpPointer>::free(ptr);
 }
 
+/***********************************************************************
+ *                    THE SAME THING, BUT SCOPED                       *
+ ***********************************************************************/
+
+/**
+ * `IndexArray<offset_t>` is the RAII form of the `copy_if_needed` /
+ * `free_if_needed` pair above, for the one job that pair is actually used
+ * for: handing a `const int64_t *` shape or stride array to a kernel
+ * templated on `offset_t`.
+ *
+ * Two things it fixes.
+ *
+ * **The manual pair leaks whenever anything between the two throws**, and in
+ * every dispatch wrapper in this tree something between them can. The shape
+ * is always
+ *
+ *     const offset_t * _size = copy_if_needed<offset_t*>(size, n);   // allocates
+ *     ... as_weights(), new reduce_t[], the impl call ...            // can throw
+ *     free_if_needed<int64_t*>(_size);                               // skipped
+ *
+ * and on the CUDA side the impl call throws by design -- every
+ * `FF_CUDA_LAUNCH` does, and so does every `copyToDevice`. Measured with
+ * ASan/LSan on the `reg_field` wrapper shape: 32 bytes in 2 allocations
+ * escape per throwing call on the narrow arm. It is a per-call leak on a
+ * path user code is expected to hit (an invalid argument), not a
+ * once-per-process one.
+ *
+ * **It does not allocate at all for the sizes that actually occur.** These
+ * arrays are `nbatch + ndim + 1` long -- single digits in every caller -- so
+ * the elements live in the object itself and the narrow arm costs no
+ * allocator traffic. Only a genuinely large rank falls back to `hostNew`.
+ * That matters most under `FF_AUTOCAST_PINNED_HOST`, where the fallback is
+ * `cudaMallocHost`: a page-locking syscall, on the per-call path, three to
+ * five times per launch.
+ *
+ * The wide arm (`offset_t == int64_t`) is specialised below to borrow the
+ * caller's array verbatim -- no copy, no storage -- exactly as
+ * `_copy_if_needed<T*, const T*>` already does for that case. So with
+ * `FF_INDEX32=0`, where both arms name `int64_t`, this whole class is a
+ * pointer copy.
+ *
+ * Deliberately not convertible to a *non*-const pointer: nothing in the
+ * dispatch layer writes through a shape or stride array, and the wide arm
+ * aliases the caller's memory.
+ */
+
+// How many elements live inside the object before it reaches for the heap.
+// `nbatch + ndim + 1` with ndim <= 3; 8 covers every shape this tree builds
+// and keeps the object at 32 bytes on the narrow arm.
+#ifndef FF_INDEX_ARRAY_INLINE
+#  define FF_INDEX_ARRAY_INLINE 8
+#endif
+
+template <class offset_t>
+class IndexArray
+{
+public:
+    // A null `src` yields a null array and copies nothing. An absent
+    // optional operand is passed as a null-data descriptor whose stride
+    // array is never read, and the call sites that have one spelled it
+    // `wgt ? copy_if_needed<offset_t*>(stride_wgt, n) : nullptr`.
+    IndexArray(const int64_t * src, size_t numel)
+        : _ptr(nullptr), _heap(nullptr)
+    {
+        if (!src) return;
+        offset_t * dst;
+        if (numel <= FF_INDEX_ARRAY_INLINE) {
+            dst = _inln;
+        } else {
+            dst = _heap = hostNew<offset_t>(numel);
+        }
+        for (size_t i = 0; i < numel; ++i)
+            dst[i] = static_cast<offset_t>(src[i]);
+        _ptr = dst;
+    }
+
+    ~IndexArray() { if (_heap) hostDelete<offset_t>(_heap); }
+
+    operator const offset_t * () const { return _ptr; }
+    const offset_t * get()       const { return _ptr; }
+
+private:
+    IndexArray(const IndexArray &);
+    IndexArray & operator=(const IndexArray &);
+
+    const offset_t * _ptr;
+    offset_t       * _heap;
+    offset_t         _inln[FF_INDEX_ARRAY_INLINE];
+};
+
+// Wide arm: there is nothing to narrow to, so borrow the caller's array.
+// Same no-op the `_copy_if_needed<T*, const T*>` specialisation provides, and
+// the reason `FF_INDEX32=0` costs nothing here either.
+template <>
+class IndexArray<int64_t>
+{
+public:
+    IndexArray(const int64_t * src, size_t /* numel */) : _ptr(src) {}
+
+    operator const int64_t * () const { return _ptr; }
+    const int64_t * get()       const { return _ptr; }
+
+private:
+    IndexArray(const IndexArray &);
+    IndexArray & operator=(const IndexArray &);
+
+    const int64_t * _ptr;
+};
+
 FF_NAMESPACE_END(FF_DEVICE)
 FF_NAMESPACE_END(FF_NS)
diff --git a/src/lib-cpu/distance.cpp b/src/lib-cpu/distance.cpp
index 45ca949..1f7ba20 100644
--- a/src/lib-cpu/distance.cpp
+++ b/src/lib-cpu/distance.cpp
@@ -51,14 +51,12 @@ inline void _dt_euclidean(
     const int64_t * size      ,     // [ndim] data shape   == (*batch, n)
     const int64_t * stride    )     // [ndim] data strides
 {
-    const offset_t * _size   = copy_if_needed<offset_t *>(size,   ndim);
-    const offset_t * _stride = copy_if_needed<offset_t *>(stride, ndim);
+    const IndexArray<offset_t> _size   (size, ndim);
+    const IndexArray<offset_t> _stride (stride, ndim);
           scalar_t * _f      = static_cast<scalar_t *>(f);
     const offset_t   _ndim   = static_cast<offset_t  >(ndim);
     const scalar_t   _w      = static_cast<scalar_t  >(w);
-    distance_e::dt(_ndim, _f, _w, _size, _stride);
-    free_if_needed<int64_t *>(_size);
-    free_if_needed<int64_t *>(_stride);
+    distance_e::dt(_ndim, _f, _w, _size.get(), _stride.get());
 }
 }
 
@@ -92,14 +90,12 @@ inline void _dt_l1(
     const int64_t * size      ,     // [ndim] data shape   == (*batch, n)
     const int64_t * stride    )     // [ndim] data strides
 {
-    const offset_t * _size   = copy_if_needed<offset_t *>(size,   ndim);
-    const offset_t * _stride = copy_if_needed<offset_t *>(stride, ndim);
+    const IndexArray<offset_t> _size   (size, ndim);
+    const IndexArray<offset_t> _stride (stride, ndim);
           scalar_t * _f      = static_cast<scalar_t *>(f);
     const offset_t   _ndim   = static_cast<offset_t  >(ndim);
     const scalar_t   _w      = static_cast<scalar_t  >(w);
-    distance_l1::dt(_ndim, _f, _w, _size, _stride);
-    free_if_needed<int64_t *>(_size);
-    free_if_needed<int64_t *>(_stride);
+    distance_l1::dt(_ndim, _f, _w, _size.get(), _stride.get());
 }
 }
 
@@ -190,12 +186,12 @@ inline void _dt_spline_table(
     //   stride_loc    -> loc   (*batch, ndim)                 == nbatch+1
     //   stride_coeff  -> coeff (*batch, npoints, ndim)        == nbatch+2
     //   stride_times  -> times (*batch, ntimes)               == nbatch+1
-    const offset_t * _size        = copy_if_needed<offset_t *>(size,         nbatch+2);
-    const offset_t * _int64_time = copy_if_needed<offset_t *>(int64_time,  nbatch);
-    const offset_t * _stride_dist = copy_if_needed<offset_t *>(stride_dist,  nbatch);
-    const offset_t * _stride_loc  = copy_if_needed<offset_t *>(stride_loc,   nbatch+1);
-    const offset_t * _stride_coeff= copy_if_needed<offset_t *>(stride_coeff, nbatch+2);
-    const offset_t * _int64_times= copy_if_needed<offset_t *>(int64_times, nbatch+1);
+    const IndexArray<offset_t> _size        (size, nbatch+2);
+    const IndexArray<offset_t> _int64_time (int64_time, nbatch);
+    const IndexArray<offset_t> _stride_dist (stride_dist, nbatch);
+    const IndexArray<offset_t> _stride_loc  (stride_loc, nbatch+1);
+    const IndexArray<offset_t> _stride_coeff(stride_coeff, nbatch+2);
+    const IndexArray<offset_t> _int64_times(int64_times, nbatch+1);
           scalar_t * _time        = static_cast<      scalar_t *>(time);
           scalar_t * _dist        = static_cast<      scalar_t *>(dist);
     const scalar_t * _loc         = static_cast<const scalar_t *>(loc);
@@ -212,12 +208,6 @@ inline void _dt_spline_table(
         _size, _int64_time, _stride_dist, _stride_loc, _stride_coeff, _int64_times,
         _spline, _bound
     );
-    free_if_needed<int64_t *>(_size);
-    free_if_needed<int64_t *>(_int64_time);
-    free_if_needed<int64_t *>(_stride_dist);
-    free_if_needed<int64_t *>(_stride_loc);
-    free_if_needed<int64_t *>(_stride_coeff);
-    free_if_needed<int64_t *>(_int64_times);
 }
 }
 
@@ -312,11 +302,11 @@ inline void _dt_spline_brent(
     //   stride_dist  -> dist  (*batch)                      == nbatch
     //   stride_loc   -> loc   (*batch, ndim)                == nbatch+1
     //   stride_coeff -> coeff (*batch, npoints, ndim)       == nbatch+2
-    const offset_t * _size        = copy_if_needed<offset_t *>(size,         nbatch+2);
-    const offset_t * _int64_time = copy_if_needed<offset_t *>(int64_time,  nbatch);
-    const offset_t * _stride_dist = copy_if_needed<offset_t *>(stride_dist,  nbatch);
-    const offset_t * _stride_loc  = copy_if_needed<offset_t *>(stride_loc,   nbatch+1);
-    const offset_t * _stride_coeff= copy_if_needed<offset_t *>(stride_coeff, nbatch+2);
+    const IndexArray<offset_t> _size        (size, nbatch+2);
+    const IndexArray<offset_t> _int64_time (int64_time, nbatch);
+    const IndexArray<offset_t> _stride_dist (stride_dist, nbatch);
+    const IndexArray<offset_t> _stride_loc  (stride_loc, nbatch+1);
+    const IndexArray<offset_t> _stride_coeff(stride_coeff, nbatch+2);
           scalar_t * _time        = static_cast<      scalar_t *>(time);
           scalar_t * _dist        = static_cast<      scalar_t *>(dist);
     const scalar_t * _loc         = static_cast<const scalar_t *>(loc);
@@ -334,11 +324,6 @@ inline void _dt_spline_brent(
         _size, _int64_time, _stride_dist, _stride_loc, _stride_coeff,
         _max_iter, _tol, _step, _spline, _bound
     );
-    free_if_needed<int64_t *>(_size);
-    free_if_needed<int64_t *>(_int64_time);
-    free_if_needed<int64_t *>(_stride_dist);
-    free_if_needed<int64_t *>(_stride_loc);
-    free_if_needed<int64_t *>(_stride_coeff);
 }
 }
 
@@ -428,11 +413,11 @@ inline void _dt_spline_gaussnewton(
     //   stride_dist  -> dist  (*batch)                      == nbatch
     //   stride_loc   -> loc   (*batch, ndim)                == nbatch+1
     //   stride_coeff -> coeff (*batch, npoints, ndim)       == nbatch+2
-    const offset_t * _size        = copy_if_needed<offset_t *>(size,         nbatch+2);
-    const offset_t * _int64_time = copy_if_needed<offset_t *>(int64_time,  nbatch);
-    const offset_t * _stride_dist = copy_if_needed<offset_t *>(stride_dist,  nbatch);
-    const offset_t * _stride_loc  = copy_if_needed<offset_t *>(stride_loc,   nbatch+1);
-    const offset_t * _stride_coeff= copy_if_needed<offset_t *>(stride_coeff, nbatch+2);
+    const IndexArray<offset_t> _size        (size, nbatch+2);
+    const IndexArray<offset_t> _int64_time (int64_time, nbatch);
+    const IndexArray<offset_t> _stride_dist (stride_dist, nbatch);
+    const IndexArray<offset_t> _stride_loc  (stride_loc, nbatch+1);
+    const IndexArray<offset_t> _stride_coeff(stride_coeff, nbatch+2);
           scalar_t * _time        = static_cast<      scalar_t *>(time);
           scalar_t * _dist        = static_cast<      scalar_t *>(dist);
     const scalar_t * _loc         = static_cast<const scalar_t *>(loc);
@@ -449,11 +434,6 @@ inline void _dt_spline_gaussnewton(
         _size, _int64_time, _stride_dist, _stride_loc, _stride_coeff,
         _max_iter, _tol, _spline, _bound
     );
-    free_if_needed<int64_t *>(_size);
-    free_if_needed<int64_t *>(_int64_time);
-    free_if_needed<int64_t *>(_stride_dist);
-    free_if_needed<int64_t *>(_stride_loc);
-    free_if_needed<int64_t *>(_stride_coeff);
 }
 }
 
@@ -602,11 +582,11 @@ _dt_mesh(
 {
     // vertices (N, D) and faces (M, D) are always 2D (the impl only reads
     // stride[0] and stride[1]); their length is independent of loc's batch rank.
-    const offset_t * _size           = copy_if_needed<offset_t *>(size,           nbatch);
-    const offset_t * _stride_dist    = copy_if_needed<offset_t *>(stride_dist,    nbatch);
-    const offset_t * _stride_coord   = copy_if_needed<offset_t *>(stride_coord,   nbatch+1);
-    const offset_t * _stride_vertices= copy_if_needed<offset_t *>(stride_vertices,2);
-    const offset_t * _stride_faces   = copy_if_needed<offset_t *>(stride_faces,   2);
+    const IndexArray<offset_t> _size           (size, nbatch);
+    const IndexArray<offset_t> _stride_dist    (stride_dist, nbatch);
+    const IndexArray<offset_t> _stride_coord   (stride_coord, nbatch+1);
+    const IndexArray<offset_t> _stride_vertices(stride_vertices, 2);
+    const IndexArray<offset_t> _stride_faces   (stride_faces, 2);
           scalar_t * _dist           = static_cast<      scalar_t *>(dist);
           index_t  * _nearest_vertex = static_cast<      index_t  *>(nearest_vertex);
     const scalar_t * _coord          = static_cast<const scalar_t *>(coord);
@@ -616,11 +596,7 @@ _dt_mesh(
     const offset_t   _nb_faces       = static_cast<      offset_t  >(nb_faces);
     const offset_t   _nb_vertices    = static_cast<      offset_t  >(nb_vertices);
 
-    const offset_t * _stride_nearest = (
-        stride_nearest
-        ? copy_if_needed<offset_t *>(stride_nearest, nbatch)
-        : nullptr
-    );
+    const IndexArray<offset_t> _stride_nearest(stride_nearest, nbatch);
 
     distance_mesh::dt<
         ndim, scalar_t, index_t, offset_t
@@ -630,13 +606,6 @@ _dt_mesh(
         _stride_dist, _stride_nearest, _stride_coord, _stride_vertices, _stride_faces,
         _signed, naive
     );
-
-    free_if_needed<int64_t *>(_size);
-    free_if_needed<int64_t *>(_stride_dist);
-    free_if_needed<int64_t *>(_stride_coord);
-    free_if_needed<int64_t *>(_stride_vertices);
-    free_if_needed<int64_t *>(_stride_faces);
-    if (_stride_nearest) free_if_needed<int64_t *>(_stride_nearest);
 }
 
 void dt_mesh(
diff --git a/src/lib-cpu/posdef.cpp b/src/lib-cpu/posdef.cpp
index 8ac7b36..7d8b47a 100644
--- a/src/lib-cpu/posdef.cpp
+++ b/src/lib-cpu/posdef.cpp
@@ -98,20 +98,16 @@ inline void _sym_matvec(
     const int64_t * stride_hes ,
     const int64_t * stride_inp )
 {
-    const offset_t * _size       = copy_if_needed<offset_t *>(size,       nbatch+1);
-    const offset_t * _stride_out = copy_if_needed<offset_t *>(stride_out, nbatch+1);
-    const offset_t * _stride_hes = copy_if_needed<offset_t *>(stride_hes, nbatch+1);
-    const offset_t * _stride_inp = copy_if_needed<offset_t *>(stride_inp, nbatch+1);
+    const IndexArray<offset_t> _size       (size, nbatch+1);
+    const IndexArray<offset_t> _stride_out (stride_out, nbatch+1);
+    const IndexArray<offset_t> _stride_hes (stride_hes, nbatch+1);
+    const IndexArray<offset_t> _stride_inp (stride_inp, nbatch+1);
           scalar_t * _out = static_cast<      scalar_t *>(out);
     const scalar_t * _hes = static_cast<const scalar_t *>(hes);
     const scalar_t * _inp = static_cast<const scalar_t *>(inp);
     posdef::sym_matvec<C, reduce_t, scalar_t, offset_t>(
         static_cast<offset_t>(nbatch), static_cast<offset_t>(nchannel),
         _out, _hes, _inp, _size, _stride_out, _stride_hes, _stride_inp);
-    free_if_needed<int64_t *>(_size);
-    free_if_needed<int64_t *>(_stride_out);
-    free_if_needed<int64_t *>(_stride_hes);
-    free_if_needed<int64_t *>(_stride_inp);
 }
 }
 
@@ -157,20 +153,16 @@ inline void _sym_addmatvec_(
     const int64_t * size, const int64_t * stride_out,
     const int64_t * stride_hes, const int64_t * stride_inp)
 {
-    const offset_t * _size       = copy_if_needed<offset_t *>(size,       nbatch+1);
-    const offset_t * _stride_out = copy_if_needed<offset_t *>(stride_out, nbatch+1);
-    const offset_t * _stride_hes = copy_if_needed<offset_t *>(stride_hes, nbatch+1);
-    const offset_t * _stride_inp = copy_if_needed<offset_t *>(stride_inp, nbatch+1);
+    const IndexArray<offset_t> _size       (size, nbatch+1);
+    const IndexArray<offset_t> _stride_out (stride_out, nbatch+1);
+    const IndexArray<offset_t> _stride_hes (stride_hes, nbatch+1);
+    const IndexArray<offset_t> _stride_inp (stride_inp, nbatch+1);
           scalar_t * _out = static_cast<      scalar_t *>(out);
     const scalar_t * _hes = static_cast<const scalar_t *>(hes);
     const scalar_t * _inp = static_cast<const scalar_t *>(inp);
     posdef::sym_addmatvec_<C, reduce_t, scalar_t, offset_t>(
         static_cast<offset_t>(nbatch), static_cast<offset_t>(nchannel),
         _out, _hes, _inp, _size, _stride_out, _stride_hes, _stride_inp);
-    free_if_needed<int64_t *>(_size);
-    free_if_needed<int64_t *>(_stride_out);
-    free_if_needed<int64_t *>(_stride_hes);
-    free_if_needed<int64_t *>(_stride_inp);
 }
 
 template <int C, typename scalar_t, typename offset_t>
@@ -180,20 +172,16 @@ inline void _sym_submatvec_(
     const int64_t * size, const int64_t * stride_out,
     const int64_t * stride_hes, const int64_t * stride_inp)
 {
-    const offset_t * _size       = copy_if_needed<offset_t *>(size,       nbatch+1);
-    const offset_t * _stride_out = copy_if_needed<offset_t *>(stride_out, nbatch+1);
-    const offset_t * _stride_hes = copy_if_needed<offset_t *>(stride_hes, nbatch+1);
-    const offset_t * _stride_inp = copy_if_needed<offset_t *>(stride_inp, nbatch+1);
+    const IndexArray<offset_t> _size       (size, nbatch+1);
+    const IndexArray<offset_t> _stride_out (stride_out, nbatch+1);
+    const IndexArray<offset_t> _stride_hes (stride_hes, nbatch+1);
+    const IndexArray<offset_t> _stride_inp (stride_inp, nbatch+1);
           scalar_t * _out = static_cast<      scalar_t *>(out);
     const scalar_t * _hes = static_cast<const scalar_t *>(hes);
     const scalar_t * _inp = static_cast<const scalar_t *>(inp);
     posdef::sym_submatvec_<C, reduce_t, scalar_t, offset_t>(
         static_cast<offset_t>(nbatch), static_cast<offset_t>(nchannel),
         _out, _hes, _inp, _size, _stride_out, _stride_hes, _stride_inp);
-    free_if_needed<int64_t *>(_size);
-    free_if_needed<int64_t *>(_stride_out);
-    free_if_needed<int64_t *>(_stride_hes);
-    free_if_needed<int64_t *>(_stride_inp);
 }
 }
 
@@ -277,20 +265,16 @@ inline void _sym_matvec_backward(
     const int64_t * size, const int64_t * stride_out,
     const int64_t * stride_grd, const int64_t * stride_inp)
 {
-    const offset_t * _size       = copy_if_needed<offset_t *>(size,       nbatch+1);
-    const offset_t * _stride_out = copy_if_needed<offset_t *>(stride_out, nbatch+1);
-    const offset_t * _stride_grd = copy_if_needed<offset_t *>(stride_grd, nbatch+1);
-    const offset_t * _stride_inp = copy_if_needed<offset_t *>(stride_inp, nbatch+1);
+    const IndexArray<offset_t> _size       (size, nbatch+1);
+    const IndexArray<offset_t> _stride_out (stride_out, nbatch+1);
+    const IndexArray<offset_t> _stride_grd (stride_grd, nbatch+1);
+    const IndexArray<offset_t> _stride_inp (stride_inp, nbatch+1);
           scalar_t * _out = static_cast<      scalar_t *>(out);
     const scalar_t * _grd = static_cast<const scalar_t *>(grd);
     const scalar_t * _inp = static_cast<const scalar_t *>(inp);
     posdef::sym_matvec_backward<C, reduce_t, scalar_t, offset_t>(
         static_cast<offset_t>(nbatch), static_cast<offset_t>(nchannel),
         _out, _grd, _inp, _size, _stride_out, _stride_grd, _stride_inp);
-    free_if_needed<int64_t *>(_size);
-    free_if_needed<int64_t *>(_stride_out);
-    free_if_needed<int64_t *>(_stride_grd);
-    free_if_needed<int64_t *>(_stride_inp);
 }
 }
 
@@ -341,11 +325,11 @@ inline void _sym_solve(
     const int64_t * stride_inp, const int64_t * stride_hes,
     const int64_t * stride_wgt)
 {
-    const offset_t * _size       = copy_if_needed<offset_t *>(size,       nbatch+1);
-    const offset_t * _stride_out = copy_if_needed<offset_t *>(stride_out, nbatch+1);
-    const offset_t * _stride_inp = copy_if_needed<offset_t *>(stride_inp, nbatch+1);
-    const offset_t * _stride_hes = copy_if_needed<offset_t *>(stride_hes, nbatch+1);
-    const offset_t * _stride_wgt = wgt ? copy_if_needed<offset_t *>(stride_wgt, nbatch+1) : nullptr;
+    const IndexArray<offset_t> _size       (size, nbatch+1);
+    const IndexArray<offset_t> _stride_out (stride_out, nbatch+1);
+    const IndexArray<offset_t> _stride_inp (stride_inp, nbatch+1);
+    const IndexArray<offset_t> _stride_hes (stride_hes, nbatch+1);
+    const IndexArray<offset_t> _stride_wgt (wgt ? stride_wgt : nullptr, nbatch+1);
           scalar_t * _out = static_cast<      scalar_t *>(out);
     const scalar_t * _inp = static_cast<const scalar_t *>(inp);
     const scalar_t * _hes = static_cast<const scalar_t *>(hes);
@@ -353,11 +337,6 @@ inline void _sym_solve(
     posdef::sym_solve<reduce_t, scalar_t, offset_t>(
         static_cast<offset_t>(nbatch), static_cast<offset_t>(nchannel),
         _out, _inp, _hes, _wgt, _size, _stride_out, _stride_inp, _stride_hes, _stride_wgt);
-    free_if_needed<int64_t *>(_size);
-    free_if_needed<int64_t *>(_stride_out);
-    free_if_needed<int64_t *>(_stride_inp);
-    free_if_needed<int64_t *>(_stride_hes);
-    if (_stride_wgt) free_if_needed<int64_t *>(_stride_wgt);
 }
 }
 
@@ -412,20 +391,16 @@ inline void _sym_solve_(
     const int64_t * size, const int64_t * stride_out,
     const int64_t * stride_hes, const int64_t * stride_wgt)
 {
-    const offset_t * _size       = copy_if_needed<offset_t *>(size,       nbatch+1);
-    const offset_t * _stride_out = copy_if_needed<offset_t *>(stride_out, nbatch+1);
-    const offset_t * _stride_hes = copy_if_needed<offset_t *>(stride_hes, nbatch+1);
-    const offset_t * _stride_wgt = wgt ? copy_if_needed<offset_t *>(stride_wgt, nbatch+1) : nullptr;
+    const IndexArray<offset_t> _size       (size, nbatch+1);
+    const IndexArray<offset_t> _stride_out (stride_out, nbatch+1);
+    const IndexArray<offset_t> _stride_hes (stride_hes, nbatch+1);
+    const IndexArray<offset_t> _stride_wgt (wgt ? stride_wgt : nullptr, nbatch+1);
           scalar_t * _out = static_cast<      scalar_t *>(out);
     const scalar_t * _hes = static_cast<const scalar_t *>(hes);
     const scalar_t * _wgt = static_cast<const scalar_t *>(wgt);
     posdef::sym_solve_<reduce_t, scalar_t, offset_t>(
         static_cast<offset_t>(nbatch), static_cast<offset_t>(nchannel),
         _out, _hes, _wgt, _size, _stride_out, _stride_hes, _stride_wgt);
-    free_if_needed<int64_t *>(_size);
-    free_if_needed<int64_t *>(_stride_out);
-    free_if_needed<int64_t *>(_stride_hes);
-    if (_stride_wgt) free_if_needed<int64_t *>(_stride_wgt);
 }
 }
 
@@ -478,17 +453,14 @@ inline void _sym_invert(
           void * out, const void * hes,
     const int64_t * size, const int64_t * stride_out, const int64_t * stride_hes)
 {
-    const offset_t * _size       = copy_if_needed<offset_t *>(size,       nbatch+1);
-    const offset_t * _stride_out = copy_if_needed<offset_t *>(stride_out, nbatch+1);
-    const offset_t * _stride_hes = copy_if_needed<offset_t *>(stride_hes, nbatch+1);
+    const IndexArray<offset_t> _size       (size, nbatch+1);
+    const IndexArray<offset_t> _stride_out (stride_out, nbatch+1);
+    const IndexArray<offset_t> _stride_hes (stride_hes, nbatch+1);
           scalar_t * _out = static_cast<      scalar_t *>(out);
     const scalar_t * _hes = static_cast<const scalar_t *>(hes);
     posdef::sym_invert<reduce_t, scalar_t, offset_t>(
         static_cast<offset_t>(nbatch), static_cast<offset_t>(nchannel),
         _out, _hes, _size, _stride_out, _stride_hes);
-    free_if_needed<int64_t *>(_size);
-    free_if_needed<int64_t *>(_stride_out);
-    free_if_needed<int64_t *>(_stride_hes);
 }
 }
 
@@ -529,14 +501,12 @@ inline void _sym_invert_(
           void * hes,
     const int64_t * size, const int64_t * stride)
 {
-    const offset_t * _size   = copy_if_needed<offset_t *>(size,   nbatch+1);
-    const offset_t * _stride = copy_if_needed<offset_t *>(stride, nbatch+1);
+    const IndexArray<offset_t> _size   (size, nbatch+1);
+    const IndexArray<offset_t> _stride (stride, nbatch+1);
           scalar_t * _hes = static_cast<scalar_t *>(hes);
     posdef::sym_invert_<reduce_t, scalar_t, offset_t>(
         static_cast<offset_t>(nbatch), static_cast<offset_t>(nchannel),
         _hes, _size, _stride);
-    free_if_needed<int64_t *>(_size);
-    free_if_needed<int64_t *>(_stride);
 }
 }
 
diff --git a/src/lib-cpu/pushpull.cpp b/src/lib-cpu/pushpull.cpp
index d0f66f7..fca48ab 100644
--- a/src/lib-cpu/pushpull.cpp
+++ b/src/lib-cpu/pushpull.cpp
@@ -25,11 +25,11 @@ inline void _pull(
     const int64_t * size_grid,  const int64_t * size_splinc,
     const int64_t * stride_out, const int64_t * stride_inp, const int64_t * stride_grid)
 {
-    const offset_t * _sg  = copy_if_needed<offset_t *>(size_grid,   n1);
-    const offset_t * _ss  = copy_if_needed<offset_t *>(size_splinc, n1);
-    const offset_t * _so  = copy_if_needed<offset_t *>(stride_out,  n1);
-    const offset_t * _si  = copy_if_needed<offset_t *>(stride_inp,  n1);
-    const offset_t * _sgr = copy_if_needed<offset_t *>(stride_grid, n1);
+    const IndexArray<offset_t> _sg  (size_grid, n1);
+    const IndexArray<offset_t> _ss  (size_splinc, n1);
+    const IndexArray<offset_t> _so  (stride_out, n1);
+    const IndexArray<offset_t> _si  (stride_inp, n1);
+    const IndexArray<offset_t> _sgr (stride_grid, n1);
           scalar_t * _out  = static_cast<      scalar_t *>(out);
     const scalar_t * _inp  = static_cast<const scalar_t *>(inp);
     const scalar_t * _grid = static_cast<const scalar_t *>(grid);
@@ -37,10 +37,6 @@ inline void _pull(
     pushpull::pull<ndim, reduce_t, scalar_t, offset_t, I, B>(
         static_cast<offset_t>(nbatch), extrapolate, _out, _inp, _grid,
         _sg, _ss, _so, _si, _sgr, bvec, svec);
-
-    free_if_needed<int64_t *>(_sg);  free_if_needed<int64_t *>(_ss);
-    free_if_needed<int64_t *>(_so);  free_if_needed<int64_t *>(_si);
-    free_if_needed<int64_t *>(_sgr);
 }
 
 template <int ndim, spline::type I, bound::type B,
@@ -52,11 +48,11 @@ inline void _push(
     const int64_t * size_grid,  const int64_t * size_splinc,
     const int64_t * stride_out, const int64_t * stride_inp, const int64_t * stride_grid)
 {
-    const offset_t * _sg  = copy_if_needed<offset_t *>(size_grid,   n1);
-    const offset_t * _ss  = copy_if_needed<offset_t *>(size_splinc, n1);
-    const offset_t * _so  = copy_if_needed<offset_t *>(stride_out,  n1);
-    const offset_t * _si  = copy_if_needed<offset_t *>(stride_inp,  n1);
-    const offset_t * _sgr = copy_if_needed<offset_t *>(stride_grid, n1);
+    const IndexArray<offset_t> _sg  (size_grid, n1);
+    const IndexArray<offset_t> _ss  (size_splinc, n1);
+    const IndexArray<offset_t> _so  (stride_out, n1);
+    const IndexArray<offset_t> _si  (stride_inp, n1);
+    const IndexArray<offset_t> _sgr (stride_grid, n1);
           scalar_t * _out  = static_cast<      scalar_t *>(out);
     const scalar_t * _inp  = static_cast<const scalar_t *>(inp);
     const scalar_t * _grid = static_cast<const scalar_t *>(grid);
@@ -64,10 +60,6 @@ inline void _push(
     pushpull::push<ndim, reduce_t, scalar_t, offset_t, I, B>(
         static_cast<offset_t>(nbatch), extrapolate, _out, _inp, _grid,
         _sg, _ss, _so, _si, _sgr, bvec, svec);
-
-    free_if_needed<int64_t *>(_sg);  free_if_needed<int64_t *>(_ss);
-    free_if_needed<int64_t *>(_so);  free_if_needed<int64_t *>(_si);
-    free_if_needed<int64_t *>(_sgr);
 }
 
 template <int ndim, spline::type I, bound::type B,
@@ -79,19 +71,16 @@ inline void _count(
     const int64_t * size_grid,  const int64_t * size_splinc,
     const int64_t * stride_out, const int64_t * stride_grid)
 {
-    const offset_t * _sg  = copy_if_needed<offset_t *>(size_grid,   n1);
-    const offset_t * _ss  = copy_if_needed<offset_t *>(size_splinc, n1);
-    const offset_t * _so  = copy_if_needed<offset_t *>(stride_out,  n1);
-    const offset_t * _sgr = copy_if_needed<offset_t *>(stride_grid, n1);
+    const IndexArray<offset_t> _sg  (size_grid, n1);
+    const IndexArray<offset_t> _ss  (size_splinc, n1);
+    const IndexArray<offset_t> _so  (stride_out, n1);
+    const IndexArray<offset_t> _sgr (stride_grid, n1);
           scalar_t * _out  = static_cast<      scalar_t *>(out);
     const scalar_t * _grid = static_cast<const scalar_t *>(grid);
 
     pushpull::count<ndim, reduce_t, scalar_t, offset_t, I, B>(
         static_cast<offset_t>(nbatch), extrapolate, _out, _grid,
         _sg, _ss, _so, _sgr, bvec, svec);
-
-    free_if_needed<int64_t *>(_sg);  free_if_needed<int64_t *>(_ss);
-    free_if_needed<int64_t *>(_so);  free_if_needed<int64_t *>(_sgr);
 }
 
 // grad: out has an extra trailing (D) axis, so stride_out has length n2 = n1+1.
@@ -104,11 +93,11 @@ inline void _grad(
     const int64_t * size_grid,  const int64_t * size_splinc,
     const int64_t * stride_out, const int64_t * stride_inp, const int64_t * stride_grid)
 {
-    const offset_t * _sg  = copy_if_needed<offset_t *>(size_grid,   n1);
-    const offset_t * _ss  = copy_if_needed<offset_t *>(size_splinc, n1);
-    const offset_t * _so  = copy_if_needed<offset_t *>(stride_out,  n1 + 1);
-    const offset_t * _si  = copy_if_needed<offset_t *>(stride_inp,  n1);
-    const offset_t * _sgr = copy_if_needed<offset_t *>(stride_grid, n1);
+    const IndexArray<offset_t> _sg  (size_grid, n1);
+    const IndexArray<offset_t> _ss  (size_splinc, n1);
+    const IndexArray<offset_t> _so  (stride_out, n1 + 1);
+    const IndexArray<offset_t> _si  (stride_inp, n1);
+    const IndexArray<offset_t> _sgr (stride_grid, n1);
           scalar_t * _out  = static_cast<      scalar_t *>(out);
     const scalar_t * _inp  = static_cast<const scalar_t *>(inp);
     const scalar_t * _grid = static_cast<const scalar_t *>(grid);
@@ -121,12 +110,7 @@ inline void _grad(
         pushpull::grad<ndim, false, reduce_t, scalar_t, offset_t, I, B>(
             static_cast<offset_t>(nbatch), extrapolate, _out, _inp, _grid,
             _sg, _ss, _so, _si, _sgr, bvec, svec);
-
-    free_if_needed<int64_t *>(_sg);  free_if_needed<int64_t *>(_ss);
-    free_if_needed<int64_t *>(_so);  free_if_needed<int64_t *>(_si);
-    free_if_needed<int64_t *>(_sgr);
 }
-
 } // anonymous namespace
 
 /***********************************************************************
diff --git a/src/lib-cpu/pushpull_backward.cpp b/src/lib-cpu/pushpull_backward.cpp
index 53c8df9..505ed77 100644
--- a/src/lib-cpu/pushpull_backward.cpp
+++ b/src/lib-cpu/pushpull_backward.cpp
@@ -39,13 +39,13 @@ inline void _pull_backward(
     const int64_t * stride_inp,  const int64_t * stride_ginp,
     const int64_t * stride_grid)
 {
-    const offset_t * _sg   = copy_if_needed<offset_t *>(size_grid,   n1);
-    const offset_t * _ss   = copy_if_needed<offset_t *>(size_splinc, n1);
-    const offset_t * _so   = copy_if_needed<offset_t *>(stride_out,  n1);
-    const offset_t * _sgo  = copy_if_needed<offset_t *>(stride_gout, n1);
-    const offset_t * _si   = copy_if_needed<offset_t *>(stride_inp,  n1);
-    const offset_t * _sgi  = copy_if_needed<offset_t *>(stride_ginp, n1);
-    const offset_t * _sgr  = copy_if_needed<offset_t *>(stride_grid, n1);
+    const IndexArray<offset_t> _sg   (size_grid, n1);
+    const IndexArray<offset_t> _ss   (size_splinc, n1);
+    const IndexArray<offset_t> _so   (stride_out, n1);
+    const IndexArray<offset_t> _sgo  (stride_gout, n1);
+    const IndexArray<offset_t> _si   (stride_inp, n1);
+    const IndexArray<offset_t> _sgi  (stride_ginp, n1);
+    const IndexArray<offset_t> _sgr  (stride_grid, n1);
           scalar_t * _out  = static_cast<      scalar_t *>(out);
           scalar_t * _gout = static_cast<      scalar_t *>(gout);
     const scalar_t * _inp  = static_cast<const scalar_t *>(inp);
@@ -56,11 +56,6 @@ inline void _pull_backward(
         static_cast<offset_t>(nbatch), extrapolate, _out, _gout,
         _inp, _ginp, _grid,
         _sg, _ss, _so, _sgo, _si, _sgi, _sgr, bvec, svec);
-
-    free_if_needed<int64_t *>(_sg);   free_if_needed<int64_t *>(_ss);
-    free_if_needed<int64_t *>(_so);   free_if_needed<int64_t *>(_sgo);
-    free_if_needed<int64_t *>(_si);   free_if_needed<int64_t *>(_sgi);
-    free_if_needed<int64_t *>(_sgr);
 }
 
 template <int ndim, spline::type I, bound::type B,
@@ -75,13 +70,13 @@ inline void _push_backward(
     const int64_t * stride_inp,  const int64_t * stride_ginp,
     const int64_t * stride_grid)
 {
-    const offset_t * _sg   = copy_if_needed<offset_t *>(size_grid,   n1);
-    const offset_t * _ss   = copy_if_needed<offset_t *>(size_splinc, n1);
-    const offset_t * _so   = copy_if_needed<offset_t *>(stride_out,  n1);
-    const offset_t * _sgo  = copy_if_needed<offset_t *>(stride_gout, n1);
-    const offset_t * _si   = copy_if_needed<offset_t *>(stride_inp,  n1);
-    const offset_t * _sgi  = copy_if_needed<offset_t *>(stride_ginp, n1);
-    const offset_t * _sgr  = copy_if_needed<offset_t *>(stride_grid, n1);
+    const IndexArray<offset_t> _sg   (size_grid, n1);
+    const IndexArray<offset_t> _ss   (size_splinc, n1);
+    const IndexArray<offset_t> _so   (stride_out, n1);
+    const IndexArray<offset_t> _sgo  (stride_gout, n1);
+    const IndexArray<offset_t> _si   (stride_inp, n1);
+    const IndexArray<offset_t> _sgi  (stride_ginp, n1);
+    const IndexArray<offset_t> _sgr  (stride_grid, n1);
           scalar_t * _out  = static_cast<      scalar_t *>(out);
           scalar_t * _gout = static_cast<      scalar_t *>(gout);
     const scalar_t * _inp  = static_cast<const scalar_t *>(inp);
@@ -92,11 +87,6 @@ inline void _push_backward(
         static_cast<offset_t>(nbatch), extrapolate, _out, _gout,
         _inp, _ginp, _grid,
         _sg, _ss, _so, _sgo, _si, _sgi, _sgr, bvec, svec);
-
-    free_if_needed<int64_t *>(_sg);   free_if_needed<int64_t *>(_ss);
-    free_if_needed<int64_t *>(_so);   free_if_needed<int64_t *>(_sgo);
-    free_if_needed<int64_t *>(_si);   free_if_needed<int64_t *>(_sgi);
-    free_if_needed<int64_t *>(_sgr);
 }
 
 template <int ndim, spline::type I, bound::type B,
@@ -109,11 +99,11 @@ inline void _count_backward(
     const int64_t * stride_gout, const int64_t * stride_ginp,
     const int64_t * stride_grid)
 {
-    const offset_t * _sg   = copy_if_needed<offset_t *>(size_grid,   n1);
-    const offset_t * _ss   = copy_if_needed<offset_t *>(size_splinc, n1);
-    const offset_t * _sgo  = copy_if_needed<offset_t *>(stride_gout, n1);
-    const offset_t * _sgi  = copy_if_needed<offset_t *>(stride_ginp, n1);
-    const offset_t * _sgr  = copy_if_needed<offset_t *>(stride_grid, n1);
+    const IndexArray<offset_t> _sg   (size_grid, n1);
+    const IndexArray<offset_t> _ss   (size_splinc, n1);
+    const IndexArray<offset_t> _sgo  (stride_gout, n1);
+    const IndexArray<offset_t> _sgi  (stride_ginp, n1);
+    const IndexArray<offset_t> _sgr  (stride_grid, n1);
           scalar_t * _gout = static_cast<      scalar_t *>(gout);
     const scalar_t * _ginp = static_cast<const scalar_t *>(ginp);
     const scalar_t * _grid = static_cast<const scalar_t *>(grid);
@@ -121,10 +111,6 @@ inline void _count_backward(
     pushpull::count_backward<ndim, false, reduce_t, scalar_t, offset_t, I, B>(
         static_cast<offset_t>(nbatch), extrapolate, _gout, _ginp, _grid,
         _sg, _ss, _sgo, _sgi, _sgr, bvec, svec);
-
-    free_if_needed<int64_t *>(_sg);   free_if_needed<int64_t *>(_ss);
-    free_if_needed<int64_t *>(_sgo);  free_if_needed<int64_t *>(_sgi);
-    free_if_needed<int64_t *>(_sgr);
 }
 
 // grad_backward: `ginp` carries the extra trailing (D) axis of `grad`'s
@@ -141,13 +127,13 @@ inline void _grad_backward(
     const int64_t * stride_inp,  const int64_t * stride_ginp,
     const int64_t * stride_grid)
 {
-    const offset_t * _sg   = copy_if_needed<offset_t *>(size_grid,   n1);
-    const offset_t * _ss   = copy_if_needed<offset_t *>(size_splinc, n1);
-    const offset_t * _so   = copy_if_needed<offset_t *>(stride_out,  n1);
-    const offset_t * _sgo  = copy_if_needed<offset_t *>(stride_gout, n1);
-    const offset_t * _si   = copy_if_needed<offset_t *>(stride_inp,  n1);
-    const offset_t * _sgi  = copy_if_needed<offset_t *>(stride_ginp, n1 + 1);
-    const offset_t * _sgr  = copy_if_needed<offset_t *>(stride_grid, n1);
+    const IndexArray<offset_t> _sg   (size_grid, n1);
+    const IndexArray<offset_t> _ss   (size_splinc, n1);
+    const IndexArray<offset_t> _so   (stride_out, n1);
+    const IndexArray<offset_t> _sgo  (stride_gout, n1);
+    const IndexArray<offset_t> _si   (stride_inp, n1);
+    const IndexArray<offset_t> _sgi  (stride_ginp, n1 + 1);
+    const IndexArray<offset_t> _sgr  (stride_grid, n1);
           scalar_t * _out  = static_cast<      scalar_t *>(out);
           scalar_t * _gout = static_cast<      scalar_t *>(gout);
     const scalar_t * _inp  = static_cast<const scalar_t *>(inp);
@@ -164,13 +150,7 @@ inline void _grad_backward(
             static_cast<offset_t>(nbatch), extrapolate, _out, _gout,
             _inp, _ginp, _grid,
             _sg, _ss, _so, _sgo, _si, _sgi, _sgr, bvec, svec);
-
-    free_if_needed<int64_t *>(_sg);   free_if_needed<int64_t *>(_ss);
-    free_if_needed<int64_t *>(_so);   free_if_needed<int64_t *>(_sgo);
-    free_if_needed<int64_t *>(_si);   free_if_needed<int64_t *>(_sgi);
-    free_if_needed<int64_t *>(_sgr);
 }
-
 } // anonymous namespace
 
 /***********************************************************************
diff --git a/src/lib-cpu/resize.cpp b/src/lib-cpu/resize.cpp
index 9bc6428..ca239d9 100644
--- a/src/lib-cpu/resize.cpp
+++ b/src/lib-cpu/resize.cpp
@@ -29,10 +29,10 @@ inline void _resample(
     const int64_t * stride_inp )   // [nall] input strides
 {
     const int64_t nall = nbatch + ndim;   // == out.ndim == inp.ndim
-    const offset_t * _size_out   = copy_if_needed<offset_t *>(size_out,   nall);
-    const offset_t * _size_inp   = copy_if_needed<offset_t *>(size_inp,   nall);
-    const offset_t * _stride_out = copy_if_needed<offset_t *>(stride_out, nall);
-    const offset_t * _stride_inp = copy_if_needed<offset_t *>(stride_inp, nall);
+    const IndexArray<offset_t> _size_out   (size_out, nall);
+    const IndexArray<offset_t> _size_inp   (size_inp, nall);
+    const IndexArray<offset_t> _stride_out (stride_out, nall);
+    const IndexArray<offset_t> _stride_inp (stride_inp, nall);
           scalar_t * _out        = static_cast<      scalar_t *>(out);
     const scalar_t * _inp        = static_cast<const scalar_t *>(inp);
 
@@ -51,11 +51,6 @@ inline void _resample(
     resize::loop<ndim, scalar_t, offset_t, double, I, B>(
         static_cast<offset_t>(nbatch), _out, _inp, shift, _scale,
         _size_out, _size_inp, _stride_out, _stride_inp);
-
-    free_if_needed<int64_t *>(_size_out);
-    free_if_needed<int64_t *>(_size_inp);
-    free_if_needed<int64_t *>(_stride_out);
-    free_if_needed<int64_t *>(_stride_inp);
 }
 } // anonymous namespace
 
diff --git a/src/lib-cpu/restrict.cpp b/src/lib-cpu/restrict.cpp
index f9191db..1f45980 100644
--- a/src/lib-cpu/restrict.cpp
+++ b/src/lib-cpu/restrict.cpp
@@ -29,10 +29,10 @@ inline void _restriction(
     const int64_t * stride_inp )   // [nall] input strides
 {
     const int64_t nall = nbatch + ndim;   // == out.ndim == inp.ndim
-    const offset_t * _size_out   = copy_if_needed<offset_t *>(size_out,   nall);
-    const offset_t * _size_inp   = copy_if_needed<offset_t *>(size_inp,   nall);
-    const offset_t * _stride_out = copy_if_needed<offset_t *>(stride_out, nall);
-    const offset_t * _stride_inp = copy_if_needed<offset_t *>(stride_inp, nall);
+    const IndexArray<offset_t> _size_out   (size_out, nall);
+    const IndexArray<offset_t> _size_inp   (size_inp, nall);
+    const IndexArray<offset_t> _stride_out (stride_out, nall);
+    const IndexArray<offset_t> _stride_inp (stride_inp, nall);
           scalar_t * _out        = static_cast<      scalar_t *>(out);
     const scalar_t * _inp        = static_cast<const scalar_t *>(inp);
 
@@ -51,11 +51,6 @@ inline void _restriction(
     restrict::loop<ndim, scalar_t, offset_t, double, I, B>(
         static_cast<offset_t>(nbatch), _out, _inp, shift, _scale,
         _size_out, _size_inp, _stride_out, _stride_inp);
-
-    free_if_needed<int64_t *>(_size_out);
-    free_if_needed<int64_t *>(_size_inp);
-    free_if_needed<int64_t *>(_stride_out);
-    free_if_needed<int64_t *>(_stride_inp);
 }
 } // anonymous namespace
 
diff --git a/src/lib-cpu/solve_field.cpp b/src/lib-cpu/solve_field.cpp
index ba309c8..24b884c 100644
--- a/src/lib-cpu/solve_field.cpp
+++ b/src/lib-cpu/solve_field.cpp
@@ -76,9 +76,9 @@ inline reduce_t _dot(
     const int64_t * stride_y )
 {
     const int64_t nall1 = nall + 1;
-    const offset_t * _size     = copy_if_needed<offset_t *>(size,     nall1);
-    const offset_t * _stride_x = copy_if_needed<offset_t *>(stride_x, nall1);
-    const offset_t * _stride_y = copy_if_needed<offset_t *>(stride_y, nall1);
+    const IndexArray<offset_t> _size     (size, nall1);
+    const IndexArray<offset_t> _stride_x (stride_x, nall1);
+    const IndexArray<offset_t> _stride_y (stride_y, nall1);
 
     reduce_t out = solve_field::dot<reduce_t, scalar_t, offset_t>(
         static_cast<offset_t>(nall),
@@ -86,9 +86,6 @@ inline reduce_t _dot(
         static_cast<const scalar_t *>(y),
         _size, _stride_x, _stride_y);
 
-    free_if_needed<int64_t *>(_size);
-    free_if_needed<int64_t *>(_stride_x);
-    free_if_needed<int64_t *>(_stride_y);
     return out;
 }
 
@@ -104,19 +101,15 @@ inline void _axpby_(
     const int64_t * stride_x )
 {
     const int64_t nall1 = nall + 1;
-    const offset_t * _size     = copy_if_needed<offset_t *>(size,     nall1);
-    const offset_t * _stride_y = copy_if_needed<offset_t *>(stride_y, nall1);
-    const offset_t * _stride_x = copy_if_needed<offset_t *>(stride_x, nall1);
+    const IndexArray<offset_t> _size     (size, nall1);
+    const IndexArray<offset_t> _stride_y (stride_y, nall1);
+    const IndexArray<offset_t> _stride_x (stride_x, nall1);
 
     solve_field::axpby_<reduce_t, scalar_t, offset_t>(
         static_cast<offset_t>(nall),
         static_cast<      scalar_t *>(y),
         static_cast<const scalar_t *>(x),
         a, b, _size, _stride_y, _stride_x);
-
-    free_if_needed<int64_t *>(_size);
-    free_if_needed<int64_t *>(_stride_y);
-    free_if_needed<int64_t *>(_stride_x);
 }
 
 // dtype/offset dispatch for the two primitives. Both operands are always
@@ -157,7 +150,6 @@ static inline void axpby_(DLTensor & y, const DLTensor & x,
     SOLVE_DT_SWITCH(AXPBY_CALL)
 #undef AXPBY_CALL
 }
-
 } // anonymous namespace
 
 /***********************************************************************
diff --git a/src/lib-cpu/splinc.cpp b/src/lib-cpu/splinc.cpp
index acda3d5..a685729 100644
--- a/src/lib-cpu/splinc.cpp
+++ b/src/lib-cpu/splinc.cpp
@@ -71,13 +71,11 @@ inline void _splinc(
     const double  * poles   )   // [npoles] filter poles
 {
     const int64_t    ndim    = nbatch + 1;
-    const offset_t * _size   = copy_if_needed<offset_t *>(size,   ndim);
-    const offset_t * _stride = copy_if_needed<offset_t *>(stride, ndim);
+    const IndexArray<offset_t> _size   (size, ndim);
+    const IndexArray<offset_t> _stride (stride, ndim);
           scalar_t * _inp    = static_cast<scalar_t *>(inp);
     splinc::loop<npoles, B, scalar_t, offset_t, double>(
         static_cast<offset_t>(nbatch), _inp, _size, _stride, poles);
-    free_if_needed<int64_t *>(_size);
-    free_if_needed<int64_t *>(_stride);
 }
 } // anonymous namespace
 
diff --git a/src/lib-cuda/distance.cpp b/src/lib-cuda/distance.cpp
index 299da51..baf004e 100644
--- a/src/lib-cuda/distance.cpp
+++ b/src/lib-cuda/distance.cpp
@@ -52,14 +52,12 @@ inline void _dt_euclidean(
     const int64_t * stride    ,     // [ndim] data strides
           intptr_t  stream    )     // CUDA stream (0 == default stream)
 {
-    const offset_t * _size   = copy_if_needed<offset_t *>(size,   ndim);
-    const offset_t * _stride = copy_if_needed<offset_t *>(stride, ndim);
+    const IndexArray<offset_t> _size   (size, ndim);
+    const IndexArray<offset_t> _stride (stride, ndim);
           scalar_t * _f      = static_cast<scalar_t *>(f);
     const offset_t   _ndim   = static_cast<offset_t  >(ndim);
     const scalar_t   _w      = static_cast<scalar_t  >(w);
-    distance_e::dt(_ndim, _f, _w, _size, _stride, stream);
-    free_if_needed<int64_t *>(_size);
-    free_if_needed<int64_t *>(_stride);
+    distance_e::dt(_ndim, _f, _w, _size.get(), _stride.get(), stream);
 }
 }
 
@@ -95,14 +93,12 @@ inline void _dt_l1(
     const int64_t * stride    ,     // [ndim] data strides
           intptr_t  stream    )     // CUDA stream (0 == default stream)
 {
-    const offset_t * _size   = copy_if_needed<offset_t *>(size,   ndim);
-    const offset_t * _stride = copy_if_needed<offset_t *>(stride, ndim);
+    const IndexArray<offset_t> _size   (size, ndim);
+    const IndexArray<offset_t> _stride (stride, ndim);
           scalar_t * _f      = static_cast<scalar_t *>(f);
     const offset_t   _ndim   = static_cast<offset_t  >(ndim);
     const scalar_t   _w      = static_cast<scalar_t  >(w);
-    distance_l1::dt(_ndim, _f, _w, _size, _stride, stream);
-    free_if_needed<int64_t *>(_size);
-    free_if_needed<int64_t *>(_stride);
+    distance_l1::dt(_ndim, _f, _w, _size.get(), _stride.get(), stream);
 }
 }
 
@@ -194,12 +190,12 @@ inline void _dt_spline_table(
     //   stride_loc    -> loc   (*batch, ndim)                 == nbatch+1
     //   stride_coeff  -> coeff (*batch, npoints, ndim)        == nbatch+2
     //   stride_times  -> times (*batch, ntimes)               == nbatch+1
-    const offset_t * _size        = copy_if_needed<offset_t *>(size,         nbatch+2);
-    const offset_t * _int64_time = copy_if_needed<offset_t *>(int64_time,  nbatch);
-    const offset_t * _stride_dist = copy_if_needed<offset_t *>(stride_dist,  nbatch);
-    const offset_t * _stride_loc  = copy_if_needed<offset_t *>(stride_loc,   nbatch+1);
-    const offset_t * _stride_coeff= copy_if_needed<offset_t *>(stride_coeff, nbatch+2);
-    const offset_t * _int64_times= copy_if_needed<offset_t *>(int64_times, nbatch+1);
+    const IndexArray<offset_t> _size        (size, nbatch+2);
+    const IndexArray<offset_t> _int64_time (int64_time, nbatch);
+    const IndexArray<offset_t> _stride_dist (stride_dist, nbatch);
+    const IndexArray<offset_t> _stride_loc  (stride_loc, nbatch+1);
+    const IndexArray<offset_t> _stride_coeff(stride_coeff, nbatch+2);
+    const IndexArray<offset_t> _int64_times(int64_times, nbatch+1);
           scalar_t * _time        = static_cast<      scalar_t *>(time);
           scalar_t * _dist        = static_cast<      scalar_t *>(dist);
     const scalar_t * _loc         = static_cast<const scalar_t *>(loc);
@@ -216,12 +212,6 @@ inline void _dt_spline_table(
         _size, _int64_time, _stride_dist, _stride_loc, _stride_coeff, _int64_times,
         _spline, _bound
     );
-    free_if_needed<int64_t *>(_size);
-    free_if_needed<int64_t *>(_int64_time);
-    free_if_needed<int64_t *>(_stride_dist);
-    free_if_needed<int64_t *>(_stride_loc);
-    free_if_needed<int64_t *>(_stride_coeff);
-    free_if_needed<int64_t *>(_int64_times);
 }
 }
 
@@ -316,11 +306,11 @@ inline void _dt_spline_brent(
     //   stride_dist  -> dist  (*batch)                      == nbatch
     //   stride_loc   -> loc   (*batch, ndim)                == nbatch+1
     //   stride_coeff -> coeff (*batch, npoints, ndim)       == nbatch+2
-    const offset_t * _size        = copy_if_needed<offset_t *>(size,         nbatch+2);
-    const offset_t * _int64_time = copy_if_needed<offset_t *>(int64_time,  nbatch);
-    const offset_t * _stride_dist = copy_if_needed<offset_t *>(stride_dist,  nbatch);
-    const offset_t * _stride_loc  = copy_if_needed<offset_t *>(stride_loc,   nbatch+1);
-    const offset_t * _stride_coeff= copy_if_needed<offset_t *>(stride_coeff, nbatch+2);
+    const IndexArray<offset_t> _size        (size, nbatch+2);
+    const IndexArray<offset_t> _int64_time (int64_time, nbatch);
+    const IndexArray<offset_t> _stride_dist (stride_dist, nbatch);
+    const IndexArray<offset_t> _stride_loc  (stride_loc, nbatch+1);
+    const IndexArray<offset_t> _stride_coeff(stride_coeff, nbatch+2);
           scalar_t * _time        = static_cast<      scalar_t *>(time);
           scalar_t * _dist        = static_cast<      scalar_t *>(dist);
     const scalar_t * _loc         = static_cast<const scalar_t *>(loc);
@@ -338,11 +328,6 @@ inline void _dt_spline_brent(
         _size, _int64_time, _stride_dist, _stride_loc, _stride_coeff,
         _max_iter, _tol, _step, _spline, _bound
     );
-    free_if_needed<int64_t *>(_size);
-    free_if_needed<int64_t *>(_int64_time);
-    free_if_needed<int64_t *>(_stride_dist);
-    free_if_needed<int64_t *>(_stride_loc);
-    free_if_needed<int64_t *>(_stride_coeff);
 }
 }
 
@@ -432,11 +417,11 @@ inline void _dt_spline_gaussnewton(
     //   stride_dist  -> dist  (*batch)                      == nbatch
     //   stride_loc   -> loc   (*batch, ndim)                == nbatch+1
     //   stride_coeff -> coeff (*batch, npoints, ndim)       == nbatch+2
-    const offset_t * _size        = copy_if_needed<offset_t *>(size,         nbatch+2);
-    const offset_t * _int64_time = copy_if_needed<offset_t *>(int64_time,  nbatch);
-    const offset_t * _stride_dist = copy_if_needed<offset_t *>(stride_dist,  nbatch);
-    const offset_t * _stride_loc  = copy_if_needed<offset_t *>(stride_loc,   nbatch+1);
-    const offset_t * _stride_coeff= copy_if_needed<offset_t *>(stride_coeff, nbatch+2);
+    const IndexArray<offset_t> _size        (size, nbatch+2);
+    const IndexArray<offset_t> _int64_time (int64_time, nbatch);
+    const IndexArray<offset_t> _stride_dist (stride_dist, nbatch);
+    const IndexArray<offset_t> _stride_loc  (stride_loc, nbatch+1);
+    const IndexArray<offset_t> _stride_coeff(stride_coeff, nbatch+2);
           scalar_t * _time        = static_cast<      scalar_t *>(time);
           scalar_t * _dist        = static_cast<      scalar_t *>(dist);
     const scalar_t * _loc         = static_cast<const scalar_t *>(loc);
@@ -453,11 +438,6 @@ inline void _dt_spline_gaussnewton(
         _size, _int64_time, _stride_dist, _stride_loc, _stride_coeff,
         _max_iter, _tol, _spline, _bound
     );
-    free_if_needed<int64_t *>(_size);
-    free_if_needed<int64_t *>(_int64_time);
-    free_if_needed<int64_t *>(_stride_dist);
-    free_if_needed<int64_t *>(_stride_loc);
-    free_if_needed<int64_t *>(_stride_coeff);
 }
 }
 
@@ -607,11 +587,11 @@ _dt_mesh(
 {
     // vertices (N, D) and faces (M, D) are always 2D (the impl only reads
     // stride[0] and stride[1]); their length is independent of loc's batch rank.
-    const offset_t * _size           = copy_if_needed<offset_t *>(size,           nbatch);
-    const offset_t * _stride_dist    = copy_if_needed<offset_t *>(stride_dist,    nbatch);
-    const offset_t * _stride_coord   = copy_if_needed<offset_t *>(stride_coord,   nbatch+1);
-    const offset_t * _stride_vertices= copy_if_needed<offset_t *>(stride_vertices,2);
-    const offset_t * _stride_faces   = copy_if_needed<offset_t *>(stride_faces,   2);
+    const IndexArray<offset_t> _size           (size, nbatch);
+    const IndexArray<offset_t> _stride_dist    (stride_dist, nbatch);
+    const IndexArray<offset_t> _stride_coord   (stride_coord, nbatch+1);
+    const IndexArray<offset_t> _stride_vertices(stride_vertices, 2);
+    const IndexArray<offset_t> _stride_faces   (stride_faces, 2);
           scalar_t * _dist           = static_cast<      scalar_t *>(dist);
           index_t  * _nearest_vertex = static_cast<      index_t  *>(nearest_vertex);
     const scalar_t * _coord          = static_cast<const scalar_t *>(coord);
@@ -621,11 +601,7 @@ _dt_mesh(
     const offset_t   _nb_faces       = static_cast<      offset_t  >(nb_faces);
     const offset_t   _nb_vertices    = static_cast<      offset_t  >(nb_vertices);
 
-    const offset_t * _stride_nearest = (
-        stride_nearest
-        ? copy_if_needed<offset_t *>(stride_nearest, nbatch)
-        : nullptr
-    );
+    const IndexArray<offset_t> _stride_nearest(stride_nearest, nbatch);
 
     distance_mesh::dt<
         ndim, scalar_t, index_t, offset_t
@@ -635,13 +611,6 @@ _dt_mesh(
         _stride_dist, _stride_nearest, _stride_coord, _stride_vertices, _stride_faces,
         _signed, naive, stream
     );
-
-    free_if_needed<int64_t *>(_size);
-    free_if_needed<int64_t *>(_stride_dist);
-    free_if_needed<int64_t *>(_stride_coord);
-    free_if_needed<int64_t *>(_stride_vertices);
-    free_if_needed<int64_t *>(_stride_faces);
-    if (_stride_nearest) free_if_needed<int64_t *>(_stride_nearest);
 }
 
 void dt_mesh(
diff --git a/src/lib-cuda/posdef.cpp b/src/lib-cuda/posdef.cpp
index 4596c94..0bfa457 100644
--- a/src/lib-cuda/posdef.cpp
+++ b/src/lib-cuda/posdef.cpp
@@ -98,20 +98,16 @@ inline void _sym_matvec(
     const int64_t * stride_hes ,
     const int64_t * stride_inp )
 {
-    const offset_t * _size       = copy_if_needed<offset_t *>(size,       nbatch+1);
-    const offset_t * _stride_out = copy_if_needed<offset_t *>(stride_out, nbatch+1);
-    const offset_t * _stride_hes = copy_if_needed<offset_t *>(stride_hes, nbatch+1);
-    const offset_t * _stride_inp = copy_if_needed<offset_t *>(stride_inp, nbatch+1);
+    const IndexArray<offset_t> _size       (size, nbatch+1);
+    const IndexArray<offset_t> _stride_out (stride_out, nbatch+1);
+    const IndexArray<offset_t> _stride_hes (stride_hes, nbatch+1);
+    const IndexArray<offset_t> _stride_inp (stride_inp, nbatch+1);
           scalar_t * _out = static_cast<      scalar_t *>(out);
     const scalar_t * _hes = static_cast<const scalar_t *>(hes);
     const scalar_t * _inp = static_cast<const scalar_t *>(inp);
     posdef::sym_matvec<C, reduce_t, scalar_t, offset_t>(
         static_cast<offset_t>(nbatch), static_cast<offset_t>(nchannel),
         _out, _hes, _inp, _size, _stride_out, _stride_hes, _stride_inp);
-    free_if_needed<int64_t *>(_size);
-    free_if_needed<int64_t *>(_stride_out);
-    free_if_needed<int64_t *>(_stride_hes);
-    free_if_needed<int64_t *>(_stride_inp);
 }
 }
 
@@ -157,20 +153,16 @@ inline void _sym_addmatvec_(
     const int64_t * size, const int64_t * stride_out,
     const int64_t * stride_hes, const int64_t * stride_inp)
 {
-    const offset_t * _size       = copy_if_needed<offset_t *>(size,       nbatch+1);
-    const offset_t * _stride_out = copy_if_needed<offset_t *>(stride_out, nbatch+1);
-    const offset_t * _stride_hes = copy_if_needed<offset_t *>(stride_hes, nbatch+1);
-    const offset_t * _stride_inp = copy_if_needed<offset_t *>(stride_inp, nbatch+1);
+    const IndexArray<offset_t> _size       (size, nbatch+1);
+    const IndexArray<offset_t> _stride_out (stride_out, nbatch+1);
+    const IndexArray<offset_t> _stride_hes (stride_hes, nbatch+1);
+    const IndexArray<offset_t> _stride_inp (stride_inp, nbatch+1);
           scalar_t * _out = static_cast<      scalar_t *>(out);
     const scalar_t * _hes = static_cast<const scalar_t *>(hes);
     const scalar_t * _inp = static_cast<const scalar_t *>(inp);
     posdef::sym_addmatvec_<C, reduce_t, scalar_t, offset_t>(
         static_cast<offset_t>(nbatch), static_cast<offset_t>(nchannel),
         _out, _hes, _inp, _size, _stride_out, _stride_hes, _stride_inp);
-    free_if_needed<int64_t *>(_size);
-    free_if_needed<int64_t *>(_stride_out);
-    free_if_needed<int64_t *>(_stride_hes);
-    free_if_needed<int64_t *>(_stride_inp);
 }
 
 template <int C, typename scalar_t, typename offset_t>
@@ -180,20 +172,16 @@ inline void _sym_submatvec_(
     const int64_t * size, const int64_t * stride_out,
     const int64_t * stride_hes, const int64_t * stride_inp)
 {
-    const offset_t * _size       = copy_if_needed<offset_t *>(size,       nbatch+1);
-    const offset_t * _stride_out = copy_if_needed<offset_t *>(stride_out, nbatch+1);
-    const offset_t * _stride_hes = copy_if_needed<offset_t *>(stride_hes, nbatch+1);
-    const offset_t * _stride_inp = copy_if_needed<offset_t *>(stride_inp, nbatch+1);
+    const IndexArray<offset_t> _size       (size, nbatch+1);
+    const IndexArray<offset_t> _stride_out (stride_out, nbatch+1);
+    const IndexArray<offset_t> _stride_hes (stride_hes, nbatch+1);
+    const IndexArray<offset_t> _stride_inp (stride_inp, nbatch+1);
           scalar_t * _out = static_cast<      scalar_t *>(out);
     const scalar_t * _hes = static_cast<const scalar_t *>(hes);
     const scalar_t * _inp = static_cast<const scalar_t *>(inp);
     posdef::sym_submatvec_<C, reduce_t, scalar_t, offset_t>(
         static_cast<offset_t>(nbatch), static_cast<offset_t>(nchannel),
         _out, _hes, _inp, _size, _stride_out, _stride_hes, _stride_inp);
-    free_if_needed<int64_t *>(_size);
-    free_if_needed<int64_t *>(_stride_out);
-    free_if_needed<int64_t *>(_stride_hes);
-    free_if_needed<int64_t *>(_stride_inp);
 }
 }
 
@@ -277,20 +265,16 @@ inline void _sym_matvec_backward(
     const int64_t * size, const int64_t * stride_out,
     const int64_t * stride_grd, const int64_t * stride_inp)
 {
-    const offset_t * _size       = copy_if_needed<offset_t *>(size,       nbatch+1);
-    const offset_t * _stride_out = copy_if_needed<offset_t *>(stride_out, nbatch+1);
-    const offset_t * _stride_grd = copy_if_needed<offset_t *>(stride_grd, nbatch+1);
-    const offset_t * _stride_inp = copy_if_needed<offset_t *>(stride_inp, nbatch+1);
+    const IndexArray<offset_t> _size       (size, nbatch+1);
+    const IndexArray<offset_t> _stride_out (stride_out, nbatch+1);
+    const IndexArray<offset_t> _stride_grd (stride_grd, nbatch+1);
+    const IndexArray<offset_t> _stride_inp (stride_inp, nbatch+1);
           scalar_t * _out = static_cast<      scalar_t *>(out);
     const scalar_t * _grd = static_cast<const scalar_t *>(grd);
     const scalar_t * _inp = static_cast<const scalar_t *>(inp);
     posdef::sym_matvec_backward<C, reduce_t, scalar_t, offset_t>(
         static_cast<offset_t>(nbatch), static_cast<offset_t>(nchannel),
         _out, _grd, _inp, _size, _stride_out, _stride_grd, _stride_inp);
-    free_if_needed<int64_t *>(_size);
-    free_if_needed<int64_t *>(_stride_out);
-    free_if_needed<int64_t *>(_stride_grd);
-    free_if_needed<int64_t *>(_stride_inp);
 }
 }
 
@@ -341,11 +325,11 @@ inline void _sym_solve(
     const int64_t * stride_inp, const int64_t * stride_hes,
     const int64_t * stride_wgt)
 {
-    const offset_t * _size       = copy_if_needed<offset_t *>(size,       nbatch+1);
-    const offset_t * _stride_out = copy_if_needed<offset_t *>(stride_out, nbatch+1);
-    const offset_t * _stride_inp = copy_if_needed<offset_t *>(stride_inp, nbatch+1);
-    const offset_t * _stride_hes = copy_if_needed<offset_t *>(stride_hes, nbatch+1);
-    const offset_t * _stride_wgt = wgt ? copy_if_needed<offset_t *>(stride_wgt, nbatch+1) : nullptr;
+    const IndexArray<offset_t> _size       (size, nbatch+1);
+    const IndexArray<offset_t> _stride_out (stride_out, nbatch+1);
+    const IndexArray<offset_t> _stride_inp (stride_inp, nbatch+1);
+    const IndexArray<offset_t> _stride_hes (stride_hes, nbatch+1);
+    const IndexArray<offset_t> _stride_wgt (wgt ? stride_wgt : nullptr, nbatch+1);
           scalar_t * _out = static_cast<      scalar_t *>(out);
     const scalar_t * _inp = static_cast<const scalar_t *>(inp);
     const scalar_t * _hes = static_cast<const scalar_t *>(hes);
@@ -353,11 +337,6 @@ inline void _sym_solve(
     posdef::sym_solve<reduce_t, scalar_t, offset_t>(
         static_cast<offset_t>(nbatch), static_cast<offset_t>(nchannel),
         _out, _inp, _hes, _wgt, _size, _stride_out, _stride_inp, _stride_hes, _stride_wgt);
-    free_if_needed<int64_t *>(_size);
-    free_if_needed<int64_t *>(_stride_out);
-    free_if_needed<int64_t *>(_stride_inp);
-    free_if_needed<int64_t *>(_stride_hes);
-    if (_stride_wgt) free_if_needed<int64_t *>(_stride_wgt);
 }
 }
 
@@ -412,20 +391,16 @@ inline void _sym_solve_(
     const int64_t * size, const int64_t * stride_out,
     const int64_t * stride_hes, const int64_t * stride_wgt)
 {
-    const offset_t * _size       = copy_if_needed<offset_t *>(size,       nbatch+1);
-    const offset_t * _stride_out = copy_if_needed<offset_t *>(stride_out, nbatch+1);
-    const offset_t * _stride_hes = copy_if_needed<offset_t *>(stride_hes, nbatch+1);
-    const offset_t * _stride_wgt = wgt ? copy_if_needed<offset_t *>(stride_wgt, nbatch+1) : nullptr;
+    const IndexArray<offset_t> _size       (size, nbatch+1);
+    const IndexArray<offset_t> _stride_out (stride_out, nbatch+1);
+    const IndexArray<offset_t> _stride_hes (stride_hes, nbatch+1);
+    const IndexArray<offset_t> _stride_wgt (wgt ? stride_wgt : nullptr, nbatch+1);
           scalar_t * _out = static_cast<      scalar_t *>(out);
     const scalar_t * _hes = static_cast<const scalar_t *>(hes);
     const scalar_t * _wgt = static_cast<const scalar_t *>(wgt);
     posdef::sym_solve_<reduce_t, scalar_t, offset_t>(
         static_cast<offset_t>(nbatch), static_cast<offset_t>(nchannel),
         _out, _hes, _wgt, _size, _stride_out, _stride_hes, _stride_wgt);
-    free_if_needed<int64_t *>(_size);
-    free_if_needed<int64_t *>(_stride_out);
-    free_if_needed<int64_t *>(_stride_hes);
-    if (_stride_wgt) free_if_needed<int64_t *>(_stride_wgt);
 }
 }
 
@@ -478,17 +453,14 @@ inline void _sym_invert(
           void * out, const void * hes,
     const int64_t * size, const int64_t * stride_out, const int64_t * stride_hes)
 {
-    const offset_t * _size       = copy_if_needed<offset_t *>(size,       nbatch+1);
-    const offset_t * _stride_out = copy_if_needed<offset_t *>(stride_out, nbatch+1);
-    const offset_t * _stride_hes = copy_if_needed<offset_t *>(stride_hes, nbatch+1);
+    const IndexArray<offset_t> _size       (size, nbatch+1);
+    const IndexArray<offset_t> _stride_out (stride_out, nbatch+1);
+    const IndexArray<offset_t> _stride_hes (stride_hes, nbatch+1);
           scalar_t * _out = static_cast<      scalar_t *>(out);
     const scalar_t * _hes = static_cast<const scalar_t *>(hes);
     posdef::sym_invert<reduce_t, scalar_t, offset_t>(
         static_cast<offset_t>(nbatch), static_cast<offset_t>(nchannel),
         _out, _hes, _size, _stride_out, _stride_hes);
-    free_if_needed<int64_t *>(_size);
-    free_if_needed<int64_t *>(_stride_out);
-    free_if_needed<int64_t *>(_stride_hes);
 }
 }
 
@@ -529,14 +501,12 @@ inline void _sym_invert_(
           void * hes,
     const int64_t * size, const int64_t * stride)
 {
-    const offset_t * _size   = copy_if_needed<offset_t *>(size,   nbatch+1);
-    const offset_t * _stride = copy_if_needed<offset_t *>(stride, nbatch+1);
+    const IndexArray<offset_t> _size   (size, nbatch+1);
+    const IndexArray<offset_t> _stride (stride, nbatch+1);
           scalar_t * _hes = static_cast<scalar_t *>(hes);
     posdef::sym_invert_<reduce_t, scalar_t, offset_t>(
         static_cast<offset_t>(nbatch), static_cast<offset_t>(nchannel),
         _hes, _size, _stride);
-    free_if_needed<int64_t *>(_size);
-    free_if_needed<int64_t *>(_stride);
 }
 }
 
diff --git a/src/lib-cuda/pushpull.cpp b/src/lib-cuda/pushpull.cpp
index eaa27d2..19249b4 100644
--- a/src/lib-cuda/pushpull.cpp
+++ b/src/lib-cuda/pushpull.cpp
@@ -26,11 +26,11 @@ inline void _pull(
     const int64_t * stride_out, const int64_t * stride_inp, const int64_t * stride_grid,
           intptr_t  stream)
 {
-    const offset_t * _sg  = copy_if_needed<offset_t *>(size_grid,   n1);
-    const offset_t * _ss  = copy_if_needed<offset_t *>(size_splinc, n1);
-    const offset_t * _so  = copy_if_needed<offset_t *>(stride_out,  n1);
-    const offset_t * _si  = copy_if_needed<offset_t *>(stride_inp,  n1);
-    const offset_t * _sgr = copy_if_needed<offset_t *>(stride_grid, n1);
+    const IndexArray<offset_t> _sg  (size_grid, n1);
+    const IndexArray<offset_t> _ss  (size_splinc, n1);
+    const IndexArray<offset_t> _so  (stride_out, n1);
+    const IndexArray<offset_t> _si  (stride_inp, n1);
+    const IndexArray<offset_t> _sgr (stride_grid, n1);
           scalar_t * _out  = static_cast<      scalar_t *>(out);
     const scalar_t * _inp  = static_cast<const scalar_t *>(inp);
     const scalar_t * _grid = static_cast<const scalar_t *>(grid);
@@ -38,10 +38,6 @@ inline void _pull(
     pushpull::pull<ndim, reduce_t, scalar_t, offset_t, I, B>(
         static_cast<offset_t>(nbatch), extrapolate, _out, _inp, _grid,
         _sg, _ss, _so, _si, _sgr, bvec, svec, stream);
-
-    free_if_needed<int64_t *>(_sg);  free_if_needed<int64_t *>(_ss);
-    free_if_needed<int64_t *>(_so);  free_if_needed<int64_t *>(_si);
-    free_if_needed<int64_t *>(_sgr);
 }
 
 template <int ndim, spline::type I, bound::type B,
@@ -54,11 +50,11 @@ inline void _push(
     const int64_t * stride_out, const int64_t * stride_inp, const int64_t * stride_grid,
           intptr_t  stream)
 {
-    const offset_t * _sg  = copy_if_needed<offset_t *>(size_grid,   n1);
-    const offset_t * _ss  = copy_if_needed<offset_t *>(size_splinc, n1);
-    const offset_t * _so  = copy_if_needed<offset_t *>(stride_out,  n1);
-    const offset_t * _si  = copy_if_needed<offset_t *>(stride_inp,  n1);
-    const offset_t * _sgr = copy_if_needed<offset_t *>(stride_grid, n1);
+    const IndexArray<offset_t> _sg  (size_grid, n1);
+    const IndexArray<offset_t> _ss  (size_splinc, n1);
+    const IndexArray<offset_t> _so  (stride_out, n1);
+    const IndexArray<offset_t> _si  (stride_inp, n1);
+    const IndexArray<offset_t> _sgr (stride_grid, n1);
           scalar_t * _out  = static_cast<      scalar_t *>(out);
     const scalar_t * _inp  = static_cast<const scalar_t *>(inp);
     const scalar_t * _grid = static_cast<const scalar_t *>(grid);
@@ -66,10 +62,6 @@ inline void _push(
     pushpull::push<ndim, reduce_t, scalar_t, offset_t, I, B>(
         static_cast<offset_t>(nbatch), extrapolate, _out, _inp, _grid,
         _sg, _ss, _so, _si, _sgr, bvec, svec, stream);
-
-    free_if_needed<int64_t *>(_sg);  free_if_needed<int64_t *>(_ss);
-    free_if_needed<int64_t *>(_so);  free_if_needed<int64_t *>(_si);
-    free_if_needed<int64_t *>(_sgr);
 }
 
 template <int ndim, spline::type I, bound::type B,
@@ -82,19 +74,16 @@ inline void _count(
     const int64_t * stride_out, const int64_t * stride_grid,
           intptr_t  stream)
 {
-    const offset_t * _sg  = copy_if_needed<offset_t *>(size_grid,   n1);
-    const offset_t * _ss  = copy_if_needed<offset_t *>(size_splinc, n1);
-    const offset_t * _so  = copy_if_needed<offset_t *>(stride_out,  n1);
-    const offset_t * _sgr = copy_if_needed<offset_t *>(stride_grid, n1);
+    const IndexArray<offset_t> _sg  (size_grid, n1);
+    const IndexArray<offset_t> _ss  (size_splinc, n1);
+    const IndexArray<offset_t> _so  (stride_out, n1);
+    const IndexArray<offset_t> _sgr (stride_grid, n1);
           scalar_t * _out  = static_cast<      scalar_t *>(out);
     const scalar_t * _grid = static_cast<const scalar_t *>(grid);
 
     pushpull::count<ndim, reduce_t, scalar_t, offset_t, I, B>(
         static_cast<offset_t>(nbatch), extrapolate, _out, _grid,
         _sg, _ss, _so, _sgr, bvec, svec, stream);
-
-    free_if_needed<int64_t *>(_sg);  free_if_needed<int64_t *>(_ss);
-    free_if_needed<int64_t *>(_so);  free_if_needed<int64_t *>(_sgr);
 }
 
 // grad: out has an extra trailing (D) axis, so stride_out has length n2 = n1+1.
@@ -108,11 +97,11 @@ inline void _grad(
     const int64_t * stride_out, const int64_t * stride_inp, const int64_t * stride_grid,
           intptr_t  stream)
 {
-    const offset_t * _sg  = copy_if_needed<offset_t *>(size_grid,   n1);
-    const offset_t * _ss  = copy_if_needed<offset_t *>(size_splinc, n1);
-    const offset_t * _so  = copy_if_needed<offset_t *>(stride_out,  n1 + 1);
-    const offset_t * _si  = copy_if_needed<offset_t *>(stride_inp,  n1);
-    const offset_t * _sgr = copy_if_needed<offset_t *>(stride_grid, n1);
+    const IndexArray<offset_t> _sg  (size_grid, n1);
+    const IndexArray<offset_t> _ss  (size_splinc, n1);
+    const IndexArray<offset_t> _so  (stride_out, n1 + 1);
+    const IndexArray<offset_t> _si  (stride_inp, n1);
+    const IndexArray<offset_t> _sgr (stride_grid, n1);
           scalar_t * _out  = static_cast<      scalar_t *>(out);
     const scalar_t * _inp  = static_cast<const scalar_t *>(inp);
     const scalar_t * _grid = static_cast<const scalar_t *>(grid);
@@ -125,12 +114,7 @@ inline void _grad(
         pushpull::grad<ndim, false, reduce_t, scalar_t, offset_t, I, B>(
             static_cast<offset_t>(nbatch), extrapolate, _out, _inp, _grid,
             _sg, _ss, _so, _si, _sgr, bvec, svec, stream);
-
-    free_if_needed<int64_t *>(_sg);  free_if_needed<int64_t *>(_ss);
-    free_if_needed<int64_t *>(_so);  free_if_needed<int64_t *>(_si);
-    free_if_needed<int64_t *>(_sgr);
 }
-
 } // anonymous namespace
 
 /***********************************************************************
diff --git a/src/lib-cuda/pushpull_backward.cpp b/src/lib-cuda/pushpull_backward.cpp
index d4eed80..1d340cd 100644
--- a/src/lib-cuda/pushpull_backward.cpp
+++ b/src/lib-cuda/pushpull_backward.cpp
@@ -44,13 +44,13 @@ inline void _pull_backward(
     const int64_t * stride_grid,
           intptr_t  stream)
 {
-    const offset_t * _sg   = copy_if_needed<offset_t *>(size_grid,   n1);
-    const offset_t * _ss   = copy_if_needed<offset_t *>(size_splinc, n1);
-    const offset_t * _so   = copy_if_needed<offset_t *>(stride_out,  n1);
-    const offset_t * _sgo  = copy_if_needed<offset_t *>(stride_gout, n1);
-    const offset_t * _si   = copy_if_needed<offset_t *>(stride_inp,  n1);
-    const offset_t * _sgi  = copy_if_needed<offset_t *>(stride_ginp, n1);
-    const offset_t * _sgr  = copy_if_needed<offset_t *>(stride_grid, n1);
+    const IndexArray<offset_t> _sg   (size_grid, n1);
+    const IndexArray<offset_t> _ss   (size_splinc, n1);
+    const IndexArray<offset_t> _so   (stride_out, n1);
+    const IndexArray<offset_t> _sgo  (stride_gout, n1);
+    const IndexArray<offset_t> _si   (stride_inp, n1);
+    const IndexArray<offset_t> _sgi  (stride_ginp, n1);
+    const IndexArray<offset_t> _sgr  (stride_grid, n1);
           scalar_t * _out  = static_cast<      scalar_t *>(out);
           scalar_t * _gout = static_cast<      scalar_t *>(gout);
     const scalar_t * _inp  = static_cast<const scalar_t *>(inp);
@@ -61,11 +61,6 @@ inline void _pull_backward(
         static_cast<offset_t>(nbatch), extrapolate, _out, _gout,
         _inp, _ginp, _grid,
         _sg, _ss, _so, _sgo, _si, _sgi, _sgr, bvec, svec, stream);
-
-    free_if_needed<int64_t *>(_sg);   free_if_needed<int64_t *>(_ss);
-    free_if_needed<int64_t *>(_so);   free_if_needed<int64_t *>(_sgo);
-    free_if_needed<int64_t *>(_si);   free_if_needed<int64_t *>(_sgi);
-    free_if_needed<int64_t *>(_sgr);
 }
 
 template <int ndim, spline::type I, bound::type B,
@@ -81,13 +76,13 @@ inline void _push_backward(
     const int64_t * stride_grid,
           intptr_t  stream)
 {
-    const offset_t * _sg   = copy_if_needed<offset_t *>(size_grid,   n1);
-    const offset_t * _ss   = copy_if_needed<offset_t *>(size_splinc, n1);
-    const offset_t * _so   = copy_if_needed<offset_t *>(stride_out,  n1);
-    const offset_t * _sgo  = copy_if_needed<offset_t *>(stride_gout, n1);
-    const offset_t * _si   = copy_if_needed<offset_t *>(stride_inp,  n1);
-    const offset_t * _sgi  = copy_if_needed<offset_t *>(stride_ginp, n1);
-    const offset_t * _sgr  = copy_if_needed<offset_t *>(stride_grid, n1);
+    const IndexArray<offset_t> _sg   (size_grid, n1);
+    const IndexArray<offset_t> _ss   (size_splinc, n1);
+    const IndexArray<offset_t> _so   (stride_out, n1);
+    const IndexArray<offset_t> _sgo  (stride_gout, n1);
+    const IndexArray<offset_t> _si   (stride_inp, n1);
+    const IndexArray<offset_t> _sgi  (stride_ginp, n1);
+    const IndexArray<offset_t> _sgr  (stride_grid, n1);
           scalar_t * _out  = static_cast<      scalar_t *>(out);
           scalar_t * _gout = static_cast<      scalar_t *>(gout);
     const scalar_t * _inp  = static_cast<const scalar_t *>(inp);
@@ -98,11 +93,6 @@ inline void _push_backward(
         static_cast<offset_t>(nbatch), extrapolate, _out, _gout,
         _inp, _ginp, _grid,
         _sg, _ss, _so, _sgo, _si, _sgi, _sgr, bvec, svec, stream);
-
-    free_if_needed<int64_t *>(_sg);   free_if_needed<int64_t *>(_ss);
-    free_if_needed<int64_t *>(_so);   free_if_needed<int64_t *>(_sgo);
-    free_if_needed<int64_t *>(_si);   free_if_needed<int64_t *>(_sgi);
-    free_if_needed<int64_t *>(_sgr);
 }
 
 template <int ndim, spline::type I, bound::type B,
@@ -116,11 +106,11 @@ inline void _count_backward(
     const int64_t * stride_grid,
           intptr_t  stream)
 {
-    const offset_t * _sg   = copy_if_needed<offset_t *>(size_grid,   n1);
-    const offset_t * _ss   = copy_if_needed<offset_t *>(size_splinc, n1);
-    const offset_t * _sgo  = copy_if_needed<offset_t *>(stride_gout, n1);
-    const offset_t * _sgi  = copy_if_needed<offset_t *>(stride_ginp, n1);
-    const offset_t * _sgr  = copy_if_needed<offset_t *>(stride_grid, n1);
+    const IndexArray<offset_t> _sg   (size_grid, n1);
+    const IndexArray<offset_t> _ss   (size_splinc, n1);
+    const IndexArray<offset_t> _sgo  (stride_gout, n1);
+    const IndexArray<offset_t> _sgi  (stride_ginp, n1);
+    const IndexArray<offset_t> _sgr  (stride_grid, n1);
           scalar_t * _gout = static_cast<      scalar_t *>(gout);
     const scalar_t * _ginp = static_cast<const scalar_t *>(ginp);
     const scalar_t * _grid = static_cast<const scalar_t *>(grid);
@@ -128,10 +118,6 @@ inline void _count_backward(
     pushpull::count_backward<ndim, false, reduce_t, scalar_t, offset_t, I, B>(
         static_cast<offset_t>(nbatch), extrapolate, _gout, _ginp, _grid,
         _sg, _ss, _sgo, _sgi, _sgr, bvec, svec, stream);
-
-    free_if_needed<int64_t *>(_sg);   free_if_needed<int64_t *>(_ss);
-    free_if_needed<int64_t *>(_sgo);  free_if_needed<int64_t *>(_sgi);
-    free_if_needed<int64_t *>(_sgr);
 }
 
 // grad_backward: `ginp` carries the extra trailing (D) axis of `grad`'s
@@ -149,13 +135,13 @@ inline void _grad_backward(
     const int64_t * stride_grid,
           intptr_t  stream)
 {
-    const offset_t * _sg   = copy_if_needed<offset_t *>(size_grid,   n1);
-    const offset_t * _ss   = copy_if_needed<offset_t *>(size_splinc, n1);
-    const offset_t * _so   = copy_if_needed<offset_t *>(stride_out,  n1);
-    const offset_t * _sgo  = copy_if_needed<offset_t *>(stride_gout, n1);
-    const offset_t * _si   = copy_if_needed<offset_t *>(stride_inp,  n1);
-    const offset_t * _sgi  = copy_if_needed<offset_t *>(stride_ginp, n1 + 1);
-    const offset_t * _sgr  = copy_if_needed<offset_t *>(stride_grid, n1);
+    const IndexArray<offset_t> _sg   (size_grid, n1);
+    const IndexArray<offset_t> _ss   (size_splinc, n1);
+    const IndexArray<offset_t> _so   (stride_out, n1);
+    const IndexArray<offset_t> _sgo  (stride_gout, n1);
+    const IndexArray<offset_t> _si   (stride_inp, n1);
+    const IndexArray<offset_t> _sgi  (stride_ginp, n1 + 1);
+    const IndexArray<offset_t> _sgr  (stride_grid, n1);
           scalar_t * _out  = static_cast<      scalar_t *>(out);
           scalar_t * _gout = static_cast<      scalar_t *>(gout);
     const scalar_t * _inp  = static_cast<const scalar_t *>(inp);
@@ -172,13 +158,7 @@ inline void _grad_backward(
             static_cast<offset_t>(nbatch), extrapolate, _out, _gout,
             _inp, _ginp, _grid,
             _sg, _ss, _so, _sgo, _si, _sgi, _sgr, bvec, svec, stream);
-
-    free_if_needed<int64_t *>(_sg);   free_if_needed<int64_t *>(_ss);
-    free_if_needed<int64_t *>(_so);   free_if_needed<int64_t *>(_sgo);
-    free_if_needed<int64_t *>(_si);   free_if_needed<int64_t *>(_sgi);
-    free_if_needed<int64_t *>(_sgr);
 }
-
 } // anonymous namespace
 
 /***********************************************************************
diff --git a/src/lib-cuda/resize.cpp b/src/lib-cuda/resize.cpp
index f0073ec..4fad6c6 100644
--- a/src/lib-cuda/resize.cpp
+++ b/src/lib-cuda/resize.cpp
@@ -29,10 +29,10 @@ inline void _resample(
     const int64_t * stride_inp )   // [nall] input strides
 {
     const int64_t nall = nbatch + ndim;   // == out.ndim == inp.ndim
-    const offset_t * _size_out   = copy_if_needed<offset_t *>(size_out,   nall);
-    const offset_t * _size_inp   = copy_if_needed<offset_t *>(size_inp,   nall);
-    const offset_t * _stride_out = copy_if_needed<offset_t *>(stride_out, nall);
-    const offset_t * _stride_inp = copy_if_needed<offset_t *>(stride_inp, nall);
+    const IndexArray<offset_t> _size_out   (size_out, nall);
+    const IndexArray<offset_t> _size_inp   (size_inp, nall);
+    const IndexArray<offset_t> _stride_out (stride_out, nall);
+    const IndexArray<offset_t> _stride_inp (stride_inp, nall);
           scalar_t * _out        = static_cast<      scalar_t *>(out);
     const scalar_t * _inp        = static_cast<const scalar_t *>(inp);
 
@@ -50,11 +50,6 @@ inline void _resample(
     resize::loop<ndim, scalar_t, offset_t, double, I, B>(
         static_cast<offset_t>(nbatch), _out, _inp, shift, _scale,
         _size_out, _size_inp, _stride_out, _stride_inp);
-
-    free_if_needed<int64_t *>(_size_out);
-    free_if_needed<int64_t *>(_size_inp);
-    free_if_needed<int64_t *>(_stride_out);
-    free_if_needed<int64_t *>(_stride_inp);
 }
 } // anonymous namespace
 
diff --git a/src/lib-cuda/restrict.cpp b/src/lib-cuda/restrict.cpp
index bedaf49..a684fc9 100644
--- a/src/lib-cuda/restrict.cpp
+++ b/src/lib-cuda/restrict.cpp
@@ -29,10 +29,10 @@ inline void _restriction(
     const int64_t * stride_inp )   // [nall] input strides
 {
     const int64_t nall = nbatch + ndim;   // == out.ndim == inp.ndim
-    const offset_t * _size_out   = copy_if_needed<offset_t *>(size_out,   nall);
-    const offset_t * _size_inp   = copy_if_needed<offset_t *>(size_inp,   nall);
-    const offset_t * _stride_out = copy_if_needed<offset_t *>(stride_out, nall);
-    const offset_t * _stride_inp = copy_if_needed<offset_t *>(stride_inp, nall);
+    const IndexArray<offset_t> _size_out   (size_out, nall);
+    const IndexArray<offset_t> _size_inp   (size_inp, nall);
+    const IndexArray<offset_t> _stride_out (stride_out, nall);
+    const IndexArray<offset_t> _stride_inp (stride_inp, nall);
           scalar_t * _out        = static_cast<      scalar_t *>(out);
     const scalar_t * _inp        = static_cast<const scalar_t *>(inp);
 
@@ -50,11 +50,6 @@ inline void _restriction(
     restrict::loop<ndim, scalar_t, offset_t, double, I, B>(
         static_cast<offset_t>(nbatch), _out, _inp, shift, _scale,
         _size_out, _size_inp, _stride_out, _stride_inp);
-
-    free_if_needed<int64_t *>(_size_out);
-    free_if_needed<int64_t *>(_size_inp);
-    free_if_needed<int64_t *>(_stride_out);
-    free_if_needed<int64_t *>(_stride_inp);
 }
 } // anonymous namespace
 
diff --git a/src/lib-cuda/splinc.cpp b/src/lib-cuda/splinc.cpp
index e96e443..c9ec99f 100644
--- a/src/lib-cuda/splinc.cpp
+++ b/src/lib-cuda/splinc.cpp
@@ -76,13 +76,11 @@ inline void _splinc(
     const double  * poles   )   // [npoles] filter poles
 {
     const int64_t    ndim    = nbatch + 1;
-    const offset_t * _size   = copy_if_needed<offset_t *>(size,   ndim);
-    const offset_t * _stride = copy_if_needed<offset_t *>(stride, ndim);
+    const IndexArray<offset_t> _size   (size, ndim);
+    const IndexArray<offset_t> _stride (stride, ndim);
           scalar_t * _inp    = static_cast<scalar_t *>(inp);
     splinc::loop<npoles, B, scalar_t, offset_t, double>(
         static_cast<offset_t>(nbatch), _inp, _size, _stride, poles);
-    free_if_needed<int64_t *>(_size);
-    free_if_needed<int64_t *>(_stride);
 }
 } // anonymous namespace