From 0de14d2c640af8714f60b1ec9a90bd43ecf2078a Mon Sep 17 00:00:00 2001 From: mloubout Date: Tue, 1 Sep 2026 10:00:31 -0400 Subject: [PATCH 1/9] 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 dd33047aab6e619459d7bfe8d8108fc2228fd345 Mon Sep 17 00:00:00 2001 From: mloubout Date: Tue, 1 Sep 2026 10:00:32 -0400 Subject: [PATCH 2/9] 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 621348d89f6279185226c99d76405d6e97ed1ec5 Mon Sep 17 00:00:00 2001 From: mloubout Date: Tue, 1 Sep 2026 10:18:24 -0400 Subject: [PATCH 3/9] 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 09e4a600a6c62933d1bc26ac771d15e4145d97cc Mon Sep 17 00:00:00 2001 From: mloubout Date: Tue, 1 Sep 2026 22:42:30 -0400 Subject: [PATCH 4/9] 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 53bcee56dbbd93297634730a3bdd1d2d23962e78 Mon Sep 17 00:00:00 2001 From: mloubout Date: Wed, 2 Sep 2026 21:56:20 -0400 Subject: [PATCH 5/9] 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 db5467aec1337587598f16f47dcd2d09619a6338 Mon Sep 17 00:00:00 2001 From: mloubout Date: Wed, 2 Sep 2026 22:14:30 -0400 Subject: [PATCH 6/9] 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 7bf7effbe8..2f5bd8766b 100644 --- a/devito/finite_differences/differentiable.py +++ b/devito/finite_differences/differentiable.py @@ -994,12 +994,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 a9d11a07068187083ae9389e75206bc85d200671 Mon Sep 17 00:00:00 2001 From: mloubout Date: Wed, 2 Sep 2026 22:22:06 -0400 Subject: [PATCH 7/9] 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 From 7597beebcf5f31c76cb1d75e3bbd3e91c53e4dbc Mon Sep 17 00:00:00 2001 From: mloubout Date: Fri, 4 Sep 2026 10:01:39 -0400 Subject: [PATCH 8/9] compiler: Print an Array initializer through the ordinary printer path An Array initializer already reached the printer, via `ccode` on a `ListInitializer`; what it did not do was reach it at the Array's own precision. The elements were printed with the Operator's settings, so a narrow Array sitting in an Operator whose arithmetic is left at the default width had its entries emitted at the wider type. Pass the Array's dtype to `ccode`, as `Expression` already does, rather than route the initializer around the printer through an `initvalue` hook of its own. A target that needs to spell its literals differently overrides `_print_ListInitializer`, which is the ordinary extension point. With the precision now correct at the point of printing, `_prec` no longer needs to be told whether the arithmetic was narrowed on purpose: a real literal takes the precision it is being printed at, and the `float32` floor applies only where that is not itself a float. That is the same value as before for every dtype other than `float16`. --- devito/core/operator.py | 14 +++++------ devito/finite_differences/differentiable.py | 3 --- devito/ir/cgen/printer.py | 27 ++++----------------- devito/ir/iet/visitors.py | 12 ++------- devito/operator/operator.py | 4 +-- 5 files changed, 15 insertions(+), 45 deletions(-) diff --git a/devito/core/operator.py b/devito/core/operator.py index 58feabdb6f..6cf20767b2 100644 --- a/devito/core/operator.py +++ b/devito/core/operator.py @@ -186,14 +186,6 @@ 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 @@ -210,6 +202,12 @@ class BasicOperator(Operator): See `examples/userapi/08_staggered_interp.ipynb` for a worked example. """ + HALF_ARITH = False + """ + Whether an Operator working in half precision carries out the arithmetic there + too, rounding its literals and its FD weights to half. Off by default. + """ + @classmethod def _normalize_kwargs(cls, **kwargs): # Will be populated with dummy values; this method is actually overridden diff --git a/devito/finite_differences/differentiable.py b/devito/finite_differences/differentiable.py index 2f5bd8766b..c38505414d 100644 --- a/devito/finite_differences/differentiable.py +++ b/devito/finite_differences/differentiable.py @@ -1000,9 +1000,6 @@ def __eq__(self, other): __hash__ = sympy.Basic.__hash__ def _hashable_content(self): - # 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) diff --git a/devito/ir/cgen/printer.py b/devito/ir/cgen/printer.py index 7c15e6df45..1e86eeae51 100644 --- a/devito/ir/cgen/printer.py +++ b/devito/ir/cgen/printer.py @@ -43,12 +43,6 @@ 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', @@ -86,11 +80,11 @@ def _prec(self, expr): 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 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): + # expression takes the precision it is being printed at. The + # `float32` floor applies only where that precision is not + # itself a float, so that an integer default doesn't silently + # degrade the literal + if np.issubdtype(self.dtype, np.floating): return self.dtype try: return np.promote_types(self.dtype, np.float32).type @@ -386,17 +380,6 @@ 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 89a74e92d7..1e2e9d924b 100644 --- a/devito/ir/iet/visitors.py +++ b/devito/ir/iet/visitors.py @@ -360,21 +360,13 @@ 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._gen_initvalue(obj, init)) + # printed at the Array's own precision, not the Operator's + value = c.Initializer(value, self.ccode(init, dtype=obj.dtype)) 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 diff --git a/devito/operator/operator.py b/devito/operator/operator.py index 2b996245cf..174714d0bf 100644 --- a/devito/operator/operator.py +++ b/devito/operator/operator.py @@ -811,8 +811,8 @@ 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 + # A Target may offer a second printer e.g. for Operators + # using half-precision arithmetic if self._sym_options.get('half-arith'): with suppress(AttributeError): return self._Target.HalfArithPrinter From b90c696f5c1e0c099d114db4d02f088dbb1f76bf Mon Sep 17 00:00:00 2001 From: mloubout Date: Tue, 8 Sep 2026 11:27:12 -0400 Subject: [PATCH 9/9] compiler: Floor a literal's precision except where the type is fixed Printing an initializer at the Array's dtype was right, but dropping the `float32` floor from `_prec` to get there was not: `Expression` also passes a dtype, and there it is the width the expression operates at, which does not constrain the width of a literal within it. A half Operator's updates were narrowing their literals with no opt-in. Restore the floor and let a caller waive it with `exact_prec`, which the Array initializer sets: its element type is fixed, so a literal that does not fit is ill-typed rather than merely less accurate. --- devito/ir/cgen/printer.py | 32 +++++++++++++++++++++++--------- devito/ir/iet/visitors.py | 9 ++++++--- 2 files changed, 29 insertions(+), 12 deletions(-) diff --git a/devito/ir/cgen/printer.py b/devito/ir/cgen/printer.py index 1e86eeae51..a7cb7f2cb8 100644 --- a/devito/ir/cgen/printer.py +++ b/devito/ir/cgen/printer.py @@ -37,6 +37,7 @@ class BasePrinter(CodePrinter): Options for code printing. """ _default_settings = {'compiler': None, 'dtype': np.float32, + 'exact_prec': False, **CodePrinter._default_settings} _func_prefix = {} @@ -75,16 +76,28 @@ def doprint(self, expr, assign_to=None): """ return self._print(expr) + @property + def _exact_prec(self): + """ + Whether the printing precision is a hard constraint rather than a + floor. An Array initializer sets it: the element type is fixed, so a + literal that does not fit is not merely less accurate but ill-typed. + """ + return self._settings['exact_prec'] + 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 takes the precision it is being printed at. The - # `float32` floor applies only where that precision is not - # itself a float, so that an integer default doesn't silently - # degrade the literal - if np.issubdtype(self.dtype, np.floating): + # expression is emitted at the precision it is printed at, + # floored at `float32`. The dtype an expression operates at + # does not constrain the width of a literal within it, and + # narrowing one loses accuracy for nothing; an aggregate + # initializer, whose element type is fixed, opts out via + # `exact_prec` + if self._exact_prec and \ + np.issubdtype(self.dtype, np.floating): return self.dtype try: return np.promote_types(self.dtype, np.float32).type @@ -459,8 +472,8 @@ def _print_Fallback(self, expr): @memoized_func -def get_printer(printer, dtype): - return printer(settings={'dtype': dtype}) +def get_printer(printer, dtype, exact_prec=False): + return printer(settings={'dtype': dtype, 'exact_prec': exact_prec}) def ccode(expr, printer=None, dtype=None): @@ -482,5 +495,6 @@ def ccode(expr, printer=None, dtype=None): if printer is None: from devito.passes.iet.languages.C import CPrinter printer = CPrinter - dtype = printer._default_settings['dtype'] if dtype is None else dtype - return get_printer(printer, dtype).doprint(expr, None) + defaults = printer._default_settings + dtype = defaults['dtype'] if dtype is None else dtype + return get_printer(printer, dtype, defaults['exact_prec']).doprint(expr, None) diff --git a/devito/ir/iet/visitors.py b/devito/ir/iet/visitors.py index 1e2e9d924b..d3fe308ccd 100644 --- a/devito/ir/iet/visitors.py +++ b/devito/ir/iet/visitors.py @@ -257,9 +257,12 @@ def __init__(self, *args, printer=None, **kwargs): printer = CPrinter self.printer = printer - def ccode(self, expr, dtype=None): - dtype = self.printer._default_settings['dtype'] if dtype is None else dtype - return get_printer(self.printer, dtype).doprint(expr, None) + def ccode(self, expr, dtype=None, exact_prec=None): + defaults = self.printer._default_settings + dtype = defaults['dtype'] if dtype is None else dtype + if exact_prec is None: + exact_prec = defaults['exact_prec'] + return get_printer(self.printer, dtype, exact_prec).doprint(expr, None) @property def _qualifiers_mapper(self):