From 2d66d92f55c4f50f09e903a207408dc86142d840 Mon Sep 17 00:00:00 2001 From: Max Freedom Pollard <272618364+MaxFreedomPollard@users.noreply.github.com> Date: Mon, 14 Sep 2026 20:07:25 -0400 Subject: [PATCH] misc: Promote to a complex double when mixing a complex float and a double `infer_dtype` in devito/tools/dtypes_lowering.py ranked the floating candidates by `np.dtype(i).itemsize` and returned the widest one. Since `np.complex64` and `np.float64` are both 8 bytes wide, that comparison was a tie, broken arbitrarily by set iteration order, and `np.float64` won. An expression mixing a single precision complex Function with a double one was therefore inferred to be real: compiler generated temporaries were declared `double` instead of `double _Complex`, `csin` degraded to `sin`, and the imaginary part was dropped at runtime with no warning. Ranked the candidates with `np.result_type` instead, which is numpy's own promotion rule and already the answer the docstring describes. Among float16, float32, float64, complex64 and complex128, `np.complex64` with `np.float64` is the only pair whose result changes. The integer branch and the single floating dtype branch were left alone. --- devito/tools/dtypes_lowering.py | 5 ++++- tests/test_dtypes.py | 32 +++++++++++++++++++++++++++++++- 2 files changed, 35 insertions(+), 2 deletions(-) diff --git a/devito/tools/dtypes_lowering.py b/devito/tools/dtypes_lowering.py index f3033aa6e9..4ae7c3d9d0 100644 --- a/devito/tools/dtypes_lowering.py +++ b/devito/tools/dtypes_lowering.py @@ -373,7 +373,10 @@ def infer_dtype(dtypes): fdtypes = {i for i in dtypes if np.issubdtype(i, np.floating) or np.issubdtype(i, np.complexfloating)} if len(fdtypes) > 1: - return max(fdtypes, key=lambda i: np.dtype(i).itemsize) + # NOTE: ranking by itemsize alone is not enough, since e.g. np.complex64 + # and np.float64 are both 8 bytes wide, yet only np.complex128 can hold + # the result of an operation between the two + return np.result_type(*fdtypes).type elif len(fdtypes) == 1: return fdtypes.pop() elif len(dtypes) == 1: diff --git a/tests/test_dtypes.py b/tests/test_dtypes.py index 861ceb17d7..d1c25d0f65 100644 --- a/tests/test_dtypes.py +++ b/tests/test_dtypes.py @@ -16,7 +16,7 @@ from devito.passes.iet.languages.openacc import AccBB, AccPrinter from devito.passes.iet.languages.openmp import OmpBB from devito.symbolics.extended_dtypes import ctypes_vector_mapper -from devito.tools import dtype_to_cstr +from devito.tools import dtype_to_cstr, infer_dtype from devito.types.basic import Basic, Scalar, Symbol from devito.types.dense import TimeFunction from devito.types.sparse import SparseTimeFunction @@ -344,3 +344,33 @@ def test_complex_reduction(dtypeu: np.dtype[np.complexfloating]) -> None: assert f'{ustr} += r0' in str(op) assert np.isclose(u.data[0, 5, 5], expected) + + +def test_complex_double_promotion() -> None: + """ + Tests that an expression mixing a single precision complex with a double + is promoted to a double precision complex, rather than to a double. + """ + # np.complex64 and np.float64 have the same itemsize, so ranking by size + # alone picks one of the two arbitrarily + assert infer_dtype({np.complex64, np.float64}) is np.complex128 + assert infer_dtype({np.float32, np.float64}) is np.float64 + assert infer_dtype({np.complex64, np.complex128}) is np.complex128 + + grid = Grid(shape=(4, 4)) + f = Function(name='f', grid=grid, dtype=np.complex64) + g = Function(name='g', grid=grid, dtype=np.float64) + h = Function(name='h', grid=grid, dtype=np.complex128) + + assert (f * g).dtype is np.complex128 + + f.data[:] = 0.3 + 0.7j + g.data[:] = 2.0 + + # The repeated subexpression is captured by a temporary, which must be + # complex too, or the imaginary part is dropped + op = Operator(Eq(h, sin(f * g) * sin(f * g) + sin(f * g))) + op.apply() + + z = np.sin(np.complex128(f.data[0, 0]) * g.data[0, 0]) + assert np.allclose(h.data, z * z + z)