From 52f5fe5c300531e9543b2ecbc7c1318bfb0cff15 Mon Sep 17 00:00:00 2001 From: mloubout Date: Tue, 1 Sep 2026 10:00:31 -0400 Subject: [PATCH 1/7] dsl: Type the FD weights after the expression they differentiate The Weights of a non-expanded derivative were always built at the default precision, so a `float16` stencil got `float` coefficients. Every wavefield*weight product then bound to the mixed-precision operators and was promoted, defeating the point of the half-precision wavefield. --- devito/finite_differences/finite_difference.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/devito/finite_differences/finite_difference.py b/devito/finite_differences/finite_difference.py index de2e92898d..73becb8b84 100644 --- a/devito/finite_differences/finite_difference.py +++ b/devito/finite_differences/finite_difference.py @@ -223,7 +223,8 @@ def make_derivative(expr, dim, fd_order, deriv_order, side, matvec, x0, coeffici expand = expand(dim) if not expand and indices.expr is not None: - weights = Weights(name='w', dimensions=indices.free_dim, initvalue=weights) + weights = Weights(name='w', dimensions=indices.free_dim, + initvalue=weights, dtype=expr.dtype) # Inject the StencilDimension # E.g. `x + i*h_x` into `f(x)` s.t. `f(x + i*h_x)` From 9c30eb4c28ce782afea874bfc92302f22bde4c37 Mon Sep 17 00:00:00 2001 From: mloubout Date: Tue, 1 Sep 2026 10:00:32 -0400 Subject: [PATCH 2/7] compiler: Print an Array initializer at the Array's own precision `_gen_value` printed the initializer with the printer's default dtype rather than the Array's, which stamped a `float` suffix onto the entries of a `double` Array and silently rounded them to single precision. Route it through a new `initvalue` printer hook, which also gives the targets a place to specialize an initializer whose type cannot be built from a plain literal. --- devito/ir/cgen/printer.py | 11 +++++++++++ devito/ir/iet/visitors.py | 11 ++++++++++- 2 files changed, 21 insertions(+), 1 deletion(-) diff --git a/devito/ir/cgen/printer.py b/devito/ir/cgen/printer.py index 96ae8c56ae..b49848482f 100644 --- a/devito/ir/cgen/printer.py +++ b/devito/ir/cgen/printer.py @@ -373,6 +373,17 @@ def _print_FieldFromComposite(self, expr): def _print_ListInitializer(self, expr): return f"{{{', '.join(self._print(i) for i in expr.params)}}}" + def initvalue(self, init, dtype): + """ + Print the aggregate initializer `init` of an Array of type `dtype`. + + Kept separate from `_print_ListInitializer` because a static + initializer, unlike an expression, cannot rely on implicit conversions: + some types (e.g. CUDA's `__half`) are only constructible from a literal + via a runtime call, which is illegal in that position. + """ + return self._print(init) + def _print_IndexedPointer(self, expr): base = self._print(expr.base) return f"{base}{''.join(f'[{self._print(i)}]' for i in expr.index)}" diff --git a/devito/ir/iet/visitors.py b/devito/ir/iet/visitors.py index 19a4604454..89a74e92d7 100644 --- a/devito/ir/iet/visitors.py +++ b/devito/ir/iet/visitors.py @@ -360,12 +360,21 @@ def _gen_value(self, obj, mode=1, masked=()): if obj.is_Array and obj.initvalue is not None and mode == 1: init = ListInitializer(obj.initvalue) if not obj._mem_constant or init.is_numeric: - value = c.Initializer(value, self.ccode(init)) + value = c.Initializer(value, self._gen_initvalue(obj, init)) elif obj.is_LocalObject and obj.initvalue is not None and mode == 1: value = c.Initializer(value, self.ccode(obj.initvalue)) return value + def _gen_initvalue(self, obj, init): + """ + Convert the aggregate initializer `init` of the Array `obj` into a C + string, delegating to the printer so that languages whose types cannot + be built from plain literals (e.g. CUDA's `__half`) can specialize it. + """ + printer = get_printer(self.printer, obj.dtype) + return printer.initvalue(init, obj.dtype) + def _gen_rettype(self, obj): try: return self._gen_value(obj, 0).typename From 764e7fe6a206b1dbe9313a663f85dd230c8b71f9 Mon Sep 17 00:00:00 2001 From: mloubout Date: Tue, 1 Sep 2026 10:18:24 -0400 Subject: [PATCH 3/7] compiler: Keep a real literal at a half-precision Operator's precision `_prec` floors an untyped real literal at `float32` so that an integer default doesn't degrade the arithmetic around it. That floor also caught `float16`, which is never a fallback but an explicit request, so every literal in a half-precision Operator printed one type too wide. Only apply the floor when the default is not already a real type. --- devito/ir/cgen/printer.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/devito/ir/cgen/printer.py b/devito/ir/cgen/printer.py index b49848482f..941c10e661 100644 --- a/devito/ir/cgen/printer.py +++ b/devito/ir/cgen/printer.py @@ -79,6 +79,13 @@ def _prec(self, expr): dtype = sympy_dtype(expr, default=self.dtype) if dtype is None or np.issubdtype(dtype, np.integer): if any(isinstance(i, Float) for i in expr.atoms()): + # A real literal in an otherwise integer (or untyped) + # expression is emitted at the Operator's precision, floored at + # `float32` so that an integer default doesn't degrade it. + # A `float16` default is a deliberate choice though, so leave + # it alone rather than silently widening the arithmetic + if np.issubdtype(self.dtype, np.floating): + return self.dtype try: return np.promote_types(self.dtype, np.float32).type except np.exceptions.DTypePromotionError: From 81ff00e8286949f3fbf65bbcaf2904bba38e40f4 Mon Sep 17 00:00:00 2001 From: mloubout Date: Tue, 1 Sep 2026 22:42:30 -0400 Subject: [PATCH 4/7] compiler: Check stability with an accumulator wider than the field The stability check sums the whole field and asks whether the result is finite. The accumulator took the field's own dtype, so in half precision it overflowed within a few thousand points and reported an instability that wasn't there -- making `errctl=max`, the very option one reaches for to diagnose a suspected instability, unusable exactly where it is needed. Give it at least single precision. --- devito/passes/iet/errors.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/devito/passes/iet/errors.py b/devito/passes/iet/errors.py index 85bf3b93a8..d9f8be1eef 100644 --- a/devito/passes/iet/errors.py +++ b/devito/passes/iet/errors.py @@ -56,7 +56,12 @@ def _check_stability(iet, wmovs=(), rcompile=None, sregistry=None): else: continue - accumulator = Symbol(name='accumulator', dtype=f.dtype) + # The accumulator sums the whole field, so it is given at least single + # precision: in half precision it would overflow within a few thousand + # points and report an instability that isn't there + dtype = np.promote_types(f.dtype, np.float32).type + + accumulator = Symbol(name='accumulator', dtype=dtype) eqns = [Eq(accumulator, 0.0), Inc(accumulator, f.subs(f.time_dim, 0))] irs, byproduct = rcompile(eqns) From 6181932174398b0ba1436e8399d2873097332575 Mon Sep 17 00:00:00 2001 From: mloubout Date: Wed, 2 Sep 2026 21:56:20 -0400 Subject: [PATCH 5/7] compiler: Let an Operator opt into arithmetic at its own precision A real literal in an otherwise integer expression is emitted at the Operator's precision, floored at `float32` so that an integer default does not degrade it. An Operator working in half wants that floor most of the time -- half is a storage format, and the accuracy of the literals is worth more than the width of the multiply -- but not always. Give the printer a flag for it, off by default, and have `_printer` pick up a Target's second printer where one is offered. --- devito/ir/cgen/printer.py | 12 +++++++++--- devito/operator/operator.py | 5 +++++ 2 files changed, 14 insertions(+), 3 deletions(-) diff --git a/devito/ir/cgen/printer.py b/devito/ir/cgen/printer.py index 941c10e661..7c15e6df45 100644 --- a/devito/ir/cgen/printer.py +++ b/devito/ir/cgen/printer.py @@ -43,6 +43,12 @@ class BasePrinter(CodePrinter): _func_literals = {} _prec_literals = {np.float32: 'F', np.complex64: 'F'} + # Whether the arithmetic is carried at the Operator's own precision, even + # where that is narrower than `float32`. Off by default: a narrow dtype is + # a storage choice, and it takes a deliberate one to also give up the + # accuracy of the literals + _half_arith = False + _qualifiers_mapper = { 'is_extern': 'extern', 'is_const': 'const', @@ -82,9 +88,9 @@ def _prec(self, expr): # A real literal in an otherwise integer (or untyped) # expression is emitted at the Operator's precision, floored at # `float32` so that an integer default doesn't degrade it. - # A `float16` default is a deliberate choice though, so leave - # it alone rather than silently widening the arithmetic - if np.issubdtype(self.dtype, np.floating): + # A printer that has opted into narrow arithmetic keeps its own + # precision instead, rather than have the literal widen it + if self._half_arith and np.issubdtype(self.dtype, np.floating): return self.dtype try: return np.promote_types(self.dtype, np.float32).type diff --git a/devito/operator/operator.py b/devito/operator/operator.py index a57ce5bd04..26eac29196 100644 --- a/devito/operator/operator.py +++ b/devito/operator/operator.py @@ -811,6 +811,11 @@ def _soname(self): @cached_property def _printer(self): + # A Target may offer a second printer for Operators that have opted + # into carrying their precision into the arithmetic + if self._options.get('half-arith'): + with suppress(AttributeError): + return self._Target.HalfArithPrinter return self._Target.Printer @cached_property From 7912a7fa076073949060b6d4fca220188644c2a5 Mon Sep 17 00:00:00 2001 From: mloubout Date: Wed, 2 Sep 2026 22:14:30 -0400 Subject: [PATCH 6/7] dsl: Tell two Weights of different precision apart The same coefficients at two precisions are two different arrays, but neither `__eq__` nor `_hashable_content` looked at the dtype, so the first one built answered for both. An Operator asking for its weights in one precision would be handed whichever an earlier Operator had cached. Compare and hash on it. The name goes in rather than the type itself, which does not order and so cannot be sorted alongside the rest. --- devito/finite_differences/differentiable.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/devito/finite_differences/differentiable.py b/devito/finite_differences/differentiable.py index 88e3cb214e..cc3beedd07 100644 --- a/devito/finite_differences/differentiable.py +++ b/devito/finite_differences/differentiable.py @@ -979,12 +979,17 @@ def __eq__(self, other): self.name == other.name and self.dimension == other.dimension and self.indices == other.indices and + self.dtype is other.dtype and self.weights == other.weights) __hash__ = sympy.Basic.__hash__ def _hashable_content(self): - return (self.name, self.dimension, str(self.weights), self.scope) + # NOTE: `dtype` belongs here. The same coefficients at two precisions + # are two different arrays, and leaving it out has one of them fetched + # from the cache in place of the other + return (self.name, self.dimension, str(self.weights), self.scope, + np.dtype(self.dtype).name) @property def dimension(self): From 975918866f656c1bdc6419d5374c5d186f2b03bf Mon Sep 17 00:00:00 2001 From: mloubout Date: Wed, 2 Sep 2026 22:22:06 -0400 Subject: [PATCH 7/7] compiler: Make half arithmetic a symbolic option Whether an Operator working in half also computes in half decides what is calculated, not how quickly: the literals and the FD coefficients are rounded to three decimal digits. That is a mathematical choice, so it belongs with `interp-mode` in `sym_opt` rather than among the codegen options, and is validated and defaulted alongside it. --- devito/core/operator.py | 11 ++++++++++- devito/operator/operator.py | 2 +- 2 files changed, 11 insertions(+), 2 deletions(-) diff --git a/devito/core/operator.py b/devito/core/operator.py index f8f34b0a24..58feabdb6f 100644 --- a/devito/core/operator.py +++ b/devito/core/operator.py @@ -186,6 +186,14 @@ class BasicOperator(Operator): # ------------------------------------------------------------------ INTERP_MODE = 'direct' + + HALF_ARITH = False + """ + Whether an Operator working in half precision carries the arithmetic there + too, rounding its literals and its FD weights to half. Off by default: half + is a storage format, and giving up the accuracy of the coefficients as well + is a mathematical choice rather than a consequence of it. + """ """ Default for the `sym_opt={'interp-mode': ...}` option. Controls how a product of fields living at different staggered locations is mapped @@ -230,7 +238,8 @@ def _normalize_sym_kwargs(cls, **kwargs): the Operator. Returns the normalized `sym_options` dict. """ so = dict(kwargs.get('sym_options', {})) - out = {'interp-mode': so.pop('interp-mode', cls.INTERP_MODE)} + out = {'interp-mode': so.pop('interp-mode', cls.INTERP_MODE), + 'half-arith': so.pop('half-arith', cls.HALF_ARITH)} if so: raise InvalidOperator( diff --git a/devito/operator/operator.py b/devito/operator/operator.py index 26eac29196..2b996245cf 100644 --- a/devito/operator/operator.py +++ b/devito/operator/operator.py @@ -813,7 +813,7 @@ def _soname(self): def _printer(self): # A Target may offer a second printer for Operators that have opted # into carrying their precision into the arithmetic - if self._options.get('half-arith'): + if self._sym_options.get('half-arith'): with suppress(AttributeError): return self._Target.HalfArithPrinter return self._Target.Printer