diff --git a/devito/core/operator.py b/devito/core/operator.py index f8f34b0a24..6cf20767b2 100644 --- a/devito/core/operator.py +++ b/devito/core/operator.py @@ -202,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 @@ -230,7 +236,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 7bf7effbe8..c38505414d 100644 --- a/devito/finite_differences/differentiable.py +++ b/devito/finite_differences/differentiable.py @@ -994,12 +994,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..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,10 +76,29 @@ 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 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 except np.exceptions.DTypePromotionError: @@ -452,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): @@ -475,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 19a4604454..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): @@ -360,7 +363,8 @@ 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)) + # 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)) diff --git a/devito/operator/operator.py b/devito/operator/operator.py index a57ce5bd04..174714d0bf 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 e.g. for Operators + # using half-precision 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)