Fp16 constant weights

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(

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):

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)`

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)

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))
 

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

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)