diff --git a/devito/core/operator.py b/devito/core/operator.py index f8f34b0a24..0a3ead0112 100644 --- a/devito/core/operator.py +++ b/devito/core/operator.py @@ -202,6 +202,14 @@ 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 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. + """ + @classmethod def _normalize_kwargs(cls, **kwargs): # Will be populated with dummy values; this method is actually overridden @@ -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/finite_differences/differentiable.py b/devito/finite_differences/differentiable.py index 88e3cb214e..765c07ab01 100644 --- a/devito/finite_differences/differentiable.py +++ b/devito/finite_differences/differentiable.py @@ -979,12 +979,14 @@ 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) + return (self.name, self.dimension, str(self.weights), self.scope, + np.dtype(self.dtype).name) @property def dimension(self): 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)` diff --git a/devito/ir/cgen/printer.py b/devito/ir/cgen/printer.py index 96ae8c56ae..1e86eeae51 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 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 except np.exceptions.DTypePromotionError: diff --git a/devito/ir/iet/visitors.py b/devito/ir/iet/visitors.py index 19a4604454..c6aa7bfce3 100644 --- a/devito/ir/iet/visitors.py +++ b/devito/ir/iet/visitors.py @@ -360,7 +360,10 @@ 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)) + # NOTE: printed at the Array's own precision, not the + # Operator's: the two differ for a narrow Array, and it is the + # element type the initializer has to be legal against + 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)) diff --git a/devito/operator/operator.py b/devito/operator/operator.py index a57ce5bd04..2b996245cf 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._sym_options.get('half-arith'): + with suppress(AttributeError): + return self._Target.HalfArithPrinter return self._Target.Printer @cached_property 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)