antonwolfy · Oct 25, 2026
@@ -4,6 +4,22 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [0.22.0] - 2027-MM-DD
+
+### Added
+
+### Changed
+
+### Deprecated
+
+* Deprecated `dpnp.dot` and `dpnp.ndarray.dot` when `a.ndim >= 2` and `b.ndim > 2`, which place the leading dimensions of `b` in the middle of the result; use `dpnp.tensordot(a, b, axes=(-1, -2))` or `dpnp.matmul` instead [#3124](https://github.com/IntelPython/dpnp/pull/3124)
+
+### Removed
+
+### Fixed
+
+### Security
+
 ## [0.21.0] - 2026-MM-DD
 
 This release achieves 100% compliance with Python Array API specification (revision [2025.12](https://data-apis.org/array-api/2025.12/)).

@@ -1300,7 +1300,12 @@ def dot(self, b, out=None):
                [8., 8.]])
         """
 
-        return dpnp.dot(self, b, out)
+        # lazy import avoids circular imports
+        # pylint: disable=import-outside-toplevel
+        from .dpnp_iface_linearalgebra import _dot
+
+        # calling _dot keeps the deprecation warning pointing to the caller
+        return _dot(self, b, out=out)
 
     @property
     def dtype(self):

@@ -40,6 +40,8 @@
 """
 
 # pylint: disable=no-name-in-module
+import warnings
+
 import numpy
 
 import dpnp
@@ -111,6 +113,14 @@ def dot(a, b, out=None):
         Returns the dot product of `a` and `b`.
         If `out` is given, then it is returned.
 
+    Warnings
+    --------
+    Calling the function when `a` has 2 or more dimensions and `b` has more
+    than 2 dimensions is deprecated and will eventually raise an error,
+    because the leading dimensions of `b` end up in the middle of the result.
+    Use ``dpnp.tensordot(a, b, axes=(-1, -2))`` for the same result, or
+    :obj:`dpnp.matmul` for batched matrix multiplication.
+
     See Also
     --------
     :obj:`dpnp.ndarray.dot` : Equivalent method.
@@ -141,15 +151,14 @@ def dot(a, b, out=None):
     array([[4, 1],
            [2, 2]])
 
-    >>> a = np.arange(3 * 4 * 5 * 6).reshape((3, 4, 5, 6))
-    >>> b = np.arange(3 * 4 * 5 * 6)[::-1].reshape((5, 4, 6, 3))
-    >>> np.dot(a, b)[2, 3, 2, 1, 2, 2]
-    array(499128)
-    >>> sum(a[2, 3, 2, :] * b[1, 2, :, 2])
-    array(499128)
-
     """
 
+    return _dot(a, b, out=out)
+
+
+def _dot(a, b, out=None):
+    """Implementation of `dpnp.dot` shared with `dpnp.ndarray.dot`."""
+
     dpnp.check_supported_arrays_type(a, b, scalar_type=True)
 
     if out is not None:
@@ -180,6 +189,18 @@ def dot(a, b, out=None):
     # TODO: investigate usage of matmul for some possible
     # use cases instead of dpnp.tensordot
     result = dpnp.tensordot(a, b, axes=(-1, -2))
+
+    # a_ndim >= 2 here; warn once tensordot validated shapes
+    if b_ndim > 2:
+        # stacklevel=3 points to the caller of dpnp.dot/dpnp.ndarray.dot
+        warnings.warn(
+            "dpnp.dot received arrays with a.ndim >= 2 and b.ndim > 2. "
+            "This places the leading dimensions of b in the middle of the "
+            "result, which is deprecated and will eventually raise an error. "
+            "Use dpnp.tensordot(a, b, axes=(-1, -2)) instead.",
+            DeprecationWarning,
+            stacklevel=3,
+        )
     return dpnp.get_result_array(result, out, casting="no")
 
 

@@ -128,7 +128,22 @@ def test_usm_ndarray_linalg_batch(func, gen_kwargs, func_kwargs):
 # check linear algebra functions from dpnp
 # with multidimensional usm_ndarray as input
 @pytest.mark.parametrize(
-    "func", ["dot", "inner", "kron", "matmul", "outer", "tensordot", "vdot"]
+    "func",
+    [
+        # N-D dot is deprecated, the warning is checked in TestDot
+        pytest.param(
+            "dot",
+            marks=pytest.mark.filterwarnings(
+                "ignore:dpnp.dot received arrays:DeprecationWarning"
+            ),
+        ),
+        "inner",
+        "kron",
+        "matmul",
+        "outer",
+        "tensordot",
+        "vdot",
+    ],
 )
 def test_usm_ndarray_linearalgebra_batch(func):
     shape = (2, 2, 2, 2)

@@ -236,7 +236,7 @@ def test_scalar(self, dtype):
             ((4,), (4, 2)),
             ((5, 4, 3), (3,)),
             ((4,), (5, 4, 3)),
-            ((5, 3, 4), (6, 4, 2)),
+            ((5, 3, 4), (4, 2)),
         ],
         ids=[
             "0d_1d",
@@ -248,7 +248,7 @@ def test_scalar(self, dtype):
             "1d_2d",
             "3d_1d",
             "1d_3d",
-            "3d_3d",
+            "3d_2d",
         ],
     )
     def test_basic(self, dtype, shape1, shape2):
@@ -307,7 +307,7 @@ def test_out_scalar(self, dtype):
             ((4,), (4, 2), (2,)),
             ((5, 4, 3), (3,), (5, 4)),
             ((4,), (5, 4, 3), (5, 3)),
-            ((5, 3, 4), (6, 4, 2), (5, 3, 6, 2)),
+            ((5, 3, 4), (4, 2), (5, 3, 2)),
         ],
         ids=[
             "0d_1d",
@@ -319,7 +319,7 @@ def test_out_scalar(self, dtype):
             "1d_2d",
             "3d_1d",
             "1d_3d",
-            "3d_3d",
+            "3d_2d",
         ],
     )
     def test_out(self, dtype, shape1, shape2, out_shape):
@@ -394,9 +394,9 @@ def test_out_error_scalar(self, ia):
             ((3, 4), (4,), (3,)),
             ((5, 4, 3), (3,), (5, 4)),
             ((4,), (3, 4, 2), (3, 2)),
-            ((5, 3, 4), (6, 4, 2), (5, 3, 6, 2)),
+            ((5, 3, 4), (4, 2), (5, 3, 2)),
         ],
-        ids=["1d_1d", "2d_2d", "2d_1d", "3d_1d", "1d_3d", "3d_3d"],
+        ids=["1d_1d", "2d_2d", "2d_1d", "3d_1d", "1d_3d", "3d_2d"],
     )
     def test_out_error(self, shape1, shape2, out_shape):
         a = numpy.ones(shape1, dtype=numpy.int32)
@@ -430,6 +430,49 @@ def test_out_error(self, shape1, shape2, out_shape):
             with pytest.raises(ValueError):
                 numpy.dot(a, b, out=np_out)
 
+    @pytest.mark.parametrize("dtype", get_all_dtypes(no_none=True))
+    @pytest.mark.parametrize(
+        "shape1, shape2",
+        [((4, 3), (5, 3, 2)), ((5, 3, 4), (6, 4, 2))],
+        ids=["2d_3d", "3d_3d"],
+    )
+    def test_nd_deprecated(self, dtype, shape1, shape2):
+        a = generate_random_numpy_array(shape1, dtype, low=-5, high=5)
+        b = generate_random_numpy_array(shape2, dtype, low=-5, high=5)
+        ia = dpnp.array(a)
+        ib = dpnp.array(b)
+
+        if numpy_version() >= "2.6.0.dev0":
+            with pytest.warns(DeprecationWarning, match="b.ndim > 2"):
+                expected = numpy.dot(a, b)
+        else:
+            expected = numpy.dot(a, b)
+
+        with pytest.warns(DeprecationWarning, match="b.ndim > 2"):
+            result = dpnp.dot(ia, ib)
+        assert_dtype_allclose(result, expected)
+
+        # ndarray
+        with pytest.warns(DeprecationWarning, match="b.ndim > 2"):
+            result = ia.dot(ib)
+        assert_dtype_allclose(result, expected)
+
+        dp_out = dpnp.empty(expected.shape, dtype=result.dtype)
+        with pytest.warns(DeprecationWarning, match="b.ndim > 2"):
+            result = dpnp.dot(ia, ib, out=dp_out)
+        assert result is dp_out
+        assert_dtype_allclose(result, expected)
+
+    def test_nd_deprecated_stacklevel(self):
+        # the warning points to the caller of the function and the method
+        ia = dpnp.ones((4, 3))
+        ib = dpnp.ones((5, 3, 2))
+
+        with pytest.warns(DeprecationWarning) as record:
+            dpnp.dot(ia, ib)
+            ia.dot(ib)
+        assert [w.filename for w in record] == [__file__] * 2
+
 
 class TestInner:
     @pytest.mark.parametrize("dtype", get_all_dtypes(no_none=True))

@@ -13,6 +13,11 @@
 from dpnp.tests.third_party.cupy import testing
 from dpnp.tests.third_party.cupy.testing._helper import skip_if_after_baseline
 
+# both numpy.dot (NumPy 2.6+) and dpnp.dot warn for a.ndim >= 2 and b.ndim > 2
+_ignore_dot_nd_deprecation = pytest.mark.filterwarnings(
+    "ignore:.*dot received arrays with a.ndim >= 2:DeprecationWarning"
+)
+
 
 @testing.parameterize(
     *testing.product(
@@ -45,6 +50,7 @@
 class TestDot(unittest.TestCase):
     # no_int8=True is added to avoid overflow for shape=((2, 3, 4), (3, 4, 2))
     # and ((2, 3), (3, 4)) cases on cpu
+    @_ignore_dot_nd_deprecation
     @testing.for_all_dtypes_combination(["dtype_a", "dtype_b"], no_int8=True)
     @testing.numpy_cupy_allclose(type_check=has_support_aspect64())
     def test_dot(self, xp, dtype_a, dtype_b):
@@ -59,6 +65,7 @@ def test_dot(self, xp, dtype_a, dtype_b):
             b = testing.shaped_arange(shape_b, xp, dtype_b)
         return xp.dot(a, b)
 
+    @_ignore_dot_nd_deprecation
     @testing.for_float_dtypes(name="dtype_a")
     @testing.for_float_dtypes(name="dtype_b")
     @testing.for_float_dtypes(name="dtype_c")
@@ -245,13 +252,15 @@ def test_dot_vec3(self, xp, dtype):
         b = testing.shaped_arange((2,), xp, dtype)
         return xp.dot(a, b)
 
+    @_ignore_dot_nd_deprecation
     @testing.for_all_dtypes(no_int8=True)
     @testing.numpy_cupy_allclose()
     def test_transposed_dot(self, xp, dtype):
         a = testing.shaped_arange((2, 3, 4), xp, dtype).transpose(1, 0, 2)
         b = testing.shaped_arange((2, 3, 4), xp, dtype).transpose(0, 2, 1)
         return xp.dot(a, b)
 
+    @_ignore_dot_nd_deprecation
     @testing.for_all_dtypes(no_int8=True)
     @testing.numpy_cupy_allclose()
     def test_transposed_dot_with_out(self, xp, dtype):
@@ -261,6 +270,7 @@ def test_transposed_dot_with_out(self, xp, dtype):
         xp.dot(a, b, out=c)
         return c
 
+    @_ignore_dot_nd_deprecation
     @testing.for_all_dtypes()
     def test_transposed_dot_with_out_f_contiguous(self, dtype):
         for xp in (numpy, cupy):