Repository: devitocodes/devito
Issue: improve handling of precision promotion for norm calculation

in devito/builtins/arithmetic.py 31, the accumulator is set to float64 to reduce roundoff, in effect assuming the array type is float32.

If f.type is int this convert to floating arithmetic might be unexpected.

What happens if f.type is already float64? Then you'd really want to promote to float128.

An alternative approach is to add an xarg to allow the caller to specify the type of the accumulator.

Another approach is to add a function that returns the promoted type for a given type, ie
float8 -> float16
float16 -> float32
float32 -> float64
float64 -> float128


Changed files:
tests/test_builtins.py, devito/builtins/arithmetic.py, devito/builtins/utils.py

diff --git a/devito/builtins/arithmetic.py b/devito/builtins/arithmetic.py
index 122218bac..3d55c1e21 100644
--- a/devito/builtins/arithmetic.py
+++ b/devito/builtins/arithmetic.py
@@ -6,6 +6,21 @@ from devito.builtins.utils import MPIReduction
 
 __all__ = ['norm', 'sumall', 'inner', 'mmin', 'mmax']
 
+accumulator_mapper = {
+    # Integer accumulates on Float64
+    np.int8: np.float64, np.uint8: np.float64,
+    np.int16: np.float64, np.uint16: np.float64,
+    np.int32: np.float64, np.uint32: np.float64,
+    np.int64: np.float64, np.uint64: np.float64,
+    # FloatX accumulates on Float2X
+    np.float16: np.float32,
+    np.float32: np.float64,
+    # NOTE: np.float128 isn't really a thing, see for example
+    # https://github.com/numpy/numpy/issues/10288
+    # https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html#1070
+    np.float64: np.float64
+}
+
 
 @dv.switchconfig(log_level='ERROR')
 def norm(f, order=2):
@@ -28,9 +43,10 @@ def norm(f, order=2):
     # otherwise we would eventually be summing more than expected
     p, eqns = f.guard() if f.is_SparseFunction else (f, [])
 
-    s = dv.types.Symbol(name='sum', dtype=np.float64)
+    dtype = accumulator_mapper[f.dtype]
+    s = dv.types.Symbol(name='sum', dtype=dtype)
 
-    with MPIReduction(f) as mr:
+    with MPIReduction(f, dtype=dtype) as mr:
         op = dv.Operator([dv.Eq(s, 0.0)] +
                          eqns +
                          [dv.Inc(s, dv.Abs(Pow(p, order))), dv.Eq(mr.n[0], s)],
@@ -59,9 +75,10 @@ def sumall(f):
     # otherwise we would eventually be summing more than expected
     p, eqns = f.guard() if f.is_SparseFunction else (f, [])
 
-    s = dv.types.Symbol(name='sum', dtype=f.dtype)
+    dtype = accumulator_mapper[f.dtype]
+    s = dv.types.Symbol(name='sum', dtype=dtype)
 
-    with MPIReduction(f) as mr:
+    with MPIReduction(f, dtype=dtype) as mr:
         op = dv.Operator([dv.Eq(s, 0.0)] +
                          eqns +
                          [dv.Inc(s, p), dv.Eq(mr.n[0], s)],
@@ -113,9 +130,10 @@ def inner(f, g):
     # otherwise we would eventually be summing more than expected
     rhs, eqns = f.guard(f*g) if f.is_SparseFunction else (f*g, [])
 
-    s = dv.types.Symbol(name='sum', dtype=f.dtype)
+    dtype = accumulator_mapper[f.dtype]
+    s = dv.types.Symbol(name='sum', dtype=dtype)
 
-    with MPIReduction(f, g) as mr:
+    with MPIReduction(f, g, dtype=dtype) as mr:
         op = dv.Operator([dv.Eq(s, 0.0)] +
                          eqns +
                          [dv.Inc(s, rhs), dv.Eq(mr.n[0], s)],
diff --git a/devito/builtins/utils.py b/devito/builtins/utils.py
index 64f54d24b..c71d31930 100644
--- a/devito/builtins/utils.py
+++ b/devito/builtins/utils.py
@@ -9,7 +9,7 @@ class MPIReduction(object):
     A context manager to build MPI-aware reduction Operators.
     """
 
-    def __init__(self, *functions, op=dv.mpi.MPI.SUM):
+    def __init__(self, *functions, op=dv.mpi.MPI.SUM, dtype=None):
         grids = {f.grid for f in functions}
         if len(grids) == 0:
             self.grid = None
@@ -17,11 +17,14 @@ class MPIReduction(object):
             self.grid = grids.pop()
         else:
             raise ValueError("Multiple Grids found")
-        dtype = {f.dtype for f in functions}
-        if len(dtype) == 1:
-            self.dtype = dtype.pop()
+        if dtype is not None:
+            self.dtype = dtype
         else:
-            raise ValueError("Illegal mixed data types")
+            dtype = {f.dtype for f in functions}
+            if len(dtype) == 1:
+                self.dtype = dtype.pop()
+            else:
+                raise ValueError("Illegal mixed data types")
         self.v = None
         self.op = op
 
diff --git a/tests/test_builtins.py b/tests/test_builtins.py
index 418d0d544..358f54f95 100644
--- a/tests/test_builtins.py
+++ b/tests/test_builtins.py
@@ -377,6 +377,38 @@ class TestBuiltinsResult(object):
 
         f.data[:] = 1.0
 
-        assert np.isclose(norm(f),
-                          switchconfig(language='openmp')(norm)(f),
+        v = norm(f)
+        assert np.isclose(v, switchconfig(language='openmp')(norm)(f),
                           rtol=1e-5)
+        assert type(v) is f.dtype
+
+    def test_issue_1863(self):
+        grid = Grid(shape=(256, 256))
+
+        f0 = Function(name='f0', grid=grid, dtype=np.int16)
+        f1 = Function(name='f1', grid=grid, dtype=np.int32)
+        f2 = Function(name='f2', grid=grid, dtype=np.float32)
+        f3 = Function(name='f3', grid=grid, dtype=np.float64)
+
+        f0.data[:] = 1
+        f1.data[:] = 1
+        f2.data[:] = 1
+        f3.data[:] = 1
+
+        # `f0` is the critical use case here:
+        # Grid 256x256 = 2**16 = 65536 max number an int16 can represent.
+        # But the underlying accumulator is a float64, so we should be good
+        # If the accumulator were not greater than 16 bits, we would cause
+        # an overflow, which is what issue #1863 exposes
+
+        v0 = norm(f0)
+        v1 = norm(f1)
+        v2 = norm(f2)
+        v3 = norm(f3)
+        assert v0 == v1
+        assert v0 == v2
+        assert v0 == v3
+        assert type(v0) == np.int16
+        assert type(v1) == np.int32
+        assert type(v2) == np.float32
+        assert type(v3) == np.float64