Explicit shape comparison for dpnp and numpy outputs (#2295)

In the test suite, the result arrays from `dpnp` and `numpy` were
compared, but their shapes were not explicitly checked for equality. As
a result, tests could pass even if the shapes differed, such as one
being `(1,)` and the other `()`.

This PR adds an explicit check to ensure that the output shapes of
`numpy` and `dpnp` match. Additionally, tests have been updated
accordingly to pass with this new check.

Author: vtavana

CHANGELOG.md

@@ -28,6 +28,8 @@ This release achieves 100% compliance with Python Array API specification (revis
 * Updated `dpnp.einsum` to add support for `order=None` [#2411](https://github.com/IntelPython/dpnp/pull/2411)
 * Updated Python Array API specification version supported to `2024.12` [#2416](https://github.com/IntelPython/dpnp/pull/2416)
 * Removed `einsum_call` keyword from `dpnp.einsum_path` signature [#2421](https://github.com/IntelPython/dpnp/pull/2421)
+* Updated `dpnp.vdot` to return a 0-D array when one of the inputs is a scalar [#2295](https://github.com/IntelPython/dpnp/pull/2295)
+* Updated `dpnp.outer` to return the same dtype as NumPy when multiplying an array with a scalar [#2295](https://github.com/IntelPython/dpnp/pull/2295)
 * Changed `"max dimensions"` to `None` in array API capabilities [#2432](https://github.com/IntelPython/dpnp/pull/2432)
 * Updated kernel header `i0.hpp` to expose `cyl_bessel_i0` function depending on build target [#2440](https://github.com/IntelPython/dpnp/pull/2440)
 * Added MKL functions `arg`, `copysign`, `i0`, and `inv` from VM namespace to be used by implementation of the appropriate element-wise functions [#2445](https://github.com/IntelPython/dpnp/pull/2445)

dpnp/dpnp_iface_linearalgebra.py

@@ -69,16 +69,30 @@
 
 
 # TODO: implement a specific scalar-array kernel
-def _call_multiply(a, b, out=None):
-    """Call multiply function for special cases of scalar-array dots."""
+def _call_multiply(a, b, out=None, outer_calc=False):
+    """
+    Adjusted multiply function for handling special cases of scalar-array dot
+    products in linear algebra.
+
+    `dpnp.multiply` cannot directly be used for calculating scalar-array dots,
+    because the output dtype of multiply is not the same as the expected dtype
+    for scalar-array dots. For example, if `sc` is an scalar and `a` is an
+    array of type `float32`, then `dpnp.multiply(a, sc).dtype == dpnp.float32`
+    (similar to NumPy). However, for scalar-array dots, such as the dot
+    function, we need `dpnp.dot(a, sc).dtype == dpnp.float64` to align with
+    NumPy. This functions adjusts the behavior of `dpnp.multiply` function to
+    meet this requirement.
+
+    """
 
     sc, arr = (a, b) if dpnp.isscalar(a) else (b, a)
     sc_dtype = map_dtype_to_device(type(sc), arr.sycl_device)
     res_dtype = dpnp.result_type(sc_dtype, arr)
+    multiply_func = dpnp.multiply.outer if outer_calc else dpnp.multiply
     if out is not None and out.dtype == arr.dtype:
-        res = dpnp.multiply(a, b, out=out)
+        res = multiply_func(a, b, out=out)
     else:
-        res = dpnp.multiply(a, b, dtype=res_dtype)
+        res = multiply_func(a, b, dtype=res_dtype)
     return dpnp.get_result_array(res, out, casting="no")
 
 
@@ -1109,16 +1123,15 @@ def outer(a, b, out=None):
 
     dpnp.check_supported_arrays_type(a, b, scalar_type=True, all_scalars=False)
     if dpnp.isscalar(a):
-        x1 = a
         x2 = dpnp.ravel(b)[None, :]
+        result = _call_multiply(a, x2, out=out, outer_calc=True)
     elif dpnp.isscalar(b):
         x1 = dpnp.ravel(a)[:, None]
-        x2 = b
+        result = _call_multiply(x1, b, out=out, outer_calc=True)
     else:
-        x1 = dpnp.ravel(a)
-        x2 = dpnp.ravel(b)
+        result = dpnp.multiply.outer(dpnp.ravel(a), dpnp.ravel(b), out=out)
 
-    return dpnp.multiply.outer(x1, x2, out=out)
+    return result
 
 
 def tensordot(a, b, axes=2):
@@ -1288,13 +1301,13 @@ def vdot(a, b):
         if b.size != 1:
             raise ValueError("The second array should be of size one.")
         a_conj = numpy.conj(a)
-        return _call_multiply(a_conj, b)
+        return dpnp.squeeze(_call_multiply(a_conj, b))
 
     if dpnp.isscalar(b):
         if a.size != 1:
             raise ValueError("The first array should be of size one.")
         a_conj = dpnp.conj(a)
-        return _call_multiply(a_conj, b)
+        return dpnp.squeeze(_call_multiply(a_conj, b))
 
     if a.ndim == 1 and b.ndim == 1:
         return dpnp_dot(a, b, out=None, conjugate=True)

dpnp/dpnp_utils/dpnp_utils_linearalgebra.py

@@ -1108,7 +1108,7 @@ def dpnp_multiplication(
                 result = dpnp.moveaxis(result, (-2, -1), axes_res)
             elif len(axes_res) == 1:
                 result = dpnp.moveaxis(result, (-1,), axes_res)
-            return dpnp.ascontiguousarray(result)
+            return result
 
         return dpnp.asarray(result, order=order)
 

dpnp/tests/helper.py

@@ -10,13 +10,29 @@
 from . import config
 
 
+def _assert_dtype(a_dt, b_dt, check_only_type_kind=False):
+    if check_only_type_kind:
+        assert a_dt.kind == b_dt.kind, f"{a_dt.kind} != {b_dt.kind}"
+    else:
+        assert a_dt == b_dt, f"{a_dt} != {b_dt}"
+
+
+def _assert_shape(a, b):
+    # it is assumed `a` is a `dpnp.ndarray` and so it has shape attribute
+    if hasattr(b, "shape"):
+        assert a.shape == b.shape, f"{a.shape} != {b.shape}"
+    else:
+        # numpy output is scalar, then dpnp is 0-D array
+        assert a.shape == (), f"{a.shape} != ()"
+
+
 def assert_dtype_allclose(
     dpnp_arr,
     numpy_arr,
     check_type=True,
     check_only_type_kind=False,
     factor=8,
-    relative_factor=None,
+    check_shape=True,
 ):
     """
     Assert DPNP and NumPy array based on maximum dtype resolution of input arrays
@@ -37,10 +53,13 @@ def assert_dtype_allclose(
     for all data types supported by DPNP when set to True.
     It is effective only when 'check_type' is also set to True.
     The parameter `factor` scales the resolution used for comparing the arrays.
+    The parameter `check_shape`, when True (default), asserts the shape of input arrays is the same.
 
     """
 
-    list_64bit_types = [numpy.float64, numpy.complex128]
+    if check_shape:
+        _assert_shape(dpnp_arr, numpy_arr)
+
     is_inexact = lambda x: hasattr(x, "dtype") and dpnp.issubdtype(
         x.dtype, dpnp.inexact
     )
@@ -57,34 +76,32 @@ def assert_dtype_allclose(
             else -dpnp.inf
         )
         tol = factor * max(tol_dpnp, tol_numpy)
-        assert_allclose(dpnp_arr.asnumpy(), numpy_arr, atol=tol, rtol=tol)
+        assert_allclose(dpnp_arr, numpy_arr, atol=tol, rtol=tol, strict=False)
         if check_type:
+            list_64bit_types = [numpy.float64, numpy.complex128]
             numpy_arr_dtype = numpy_arr.dtype
             dpnp_arr_dtype = dpnp_arr.dtype
             dpnp_arr_dev = dpnp_arr.sycl_device
 
             if check_only_type_kind:
-                assert dpnp_arr_dtype.kind == numpy_arr_dtype.kind
+                _assert_dtype(dpnp_arr_dtype, numpy_arr_dtype, True)
             else:
                 is_np_arr_f2 = numpy_arr_dtype == numpy.float16
 
                 if is_np_arr_f2:
                     if has_support_aspect16(dpnp_arr_dev):
-                        assert dpnp_arr_dtype == numpy_arr_dtype
+                        _assert_dtype(dpnp_arr_dtype, numpy_arr_dtype)
                 elif (
                     numpy_arr_dtype not in list_64bit_types
                     or has_support_aspect64(dpnp_arr_dev)
                 ):
-                    assert dpnp_arr_dtype == numpy_arr_dtype
+                    _assert_dtype(dpnp_arr_dtype, numpy_arr_dtype)
                 else:
-                    assert dpnp_arr_dtype.kind == numpy_arr_dtype.kind
+                    _assert_dtype(dpnp_arr_dtype, numpy_arr_dtype, True)
     else:
-        assert_array_equal(dpnp_arr.asnumpy(), numpy_arr)
+        assert_array_equal(dpnp_arr, numpy_arr, strict=False)
         if check_type and hasattr(numpy_arr, "dtype"):
-            if check_only_type_kind:
-                assert dpnp_arr.dtype.kind == numpy_arr.dtype.kind
-            else:
-                assert dpnp_arr.dtype == numpy_arr.dtype
+            _assert_dtype(dpnp_arr.dtype, numpy_arr.dtype, check_only_type_kind)
 
 
 def generate_random_numpy_array(

dpnp/tests/test_arraycreation.py

@@ -952,15 +952,13 @@ def test_ascontiguousarray1(data):
     result = dpnp.ascontiguousarray(data)
     expected = numpy.ascontiguousarray(data)
     assert_dtype_allclose(result, expected)
-    assert result.shape == expected.shape
 
 
 @pytest.mark.parametrize("data", [(), 1, (2, 3), [4]])
 def test_ascontiguousarray2(data):
     result = dpnp.ascontiguousarray(dpnp.array(data))
     expected = numpy.ascontiguousarray(numpy.array(data))
     assert_dtype_allclose(result, expected)
-    assert result.shape == expected.shape
 
 
 @pytest.mark.parametrize(
@@ -970,15 +968,13 @@ def test_asfortranarray1(data):
     result = dpnp.asfortranarray(data)
     expected = numpy.asfortranarray(data)
     assert_dtype_allclose(result, expected)
-    assert result.shape == expected.shape
 
 
 @pytest.mark.parametrize("data", [(), 1, (2, 3), [4]])
 def test_asfortranarray2(data):
     result = dpnp.asfortranarray(dpnp.array(data))
     expected = numpy.asfortranarray(numpy.array(data))
     assert_dtype_allclose(result, expected)
-    assert result.shape == expected.shape
 
 
 def test_meshgrid_raise_error():

dpnp/tests/test_arraypad.py

@@ -42,7 +42,6 @@ def test_basic(self, mode):
         result = dpnp.pad(a_dp, (25, 20), mode=mode)
         if mode == "empty":
             # omit uninitialized "empty" boundary from the comparison
-            assert result.shape == expected.shape
             assert_equal(result[25:-20], expected[25:-20])
         else:
             assert_array_equal(result, expected)
@@ -70,7 +69,6 @@ def test_non_contiguous_array(self, mode):
         result = dpnp.pad(a_dp, (2, 3), mode=mode)
         if mode == "empty":
             # omit uninitialized "empty" boundary from the comparison
-            assert result.shape == expected.shape
             assert_equal(result[2:-3, 2:-3], expected[2:-3, 2:-3])
         else:
             assert_array_equal(result, expected)
@@ -287,10 +285,10 @@ def test_linear_ramp_end_values(self):
         """Ensure that end values are exact."""
         a_dp = dpnp.ones(10).reshape(2, 5)
         a = dpnp.pad(a_dp, (223, 123), mode="linear_ramp")
-        assert_equal(a[:, 0], 0.0)
-        assert_equal(a[:, -1], 0.0)
-        assert_equal(a[0, :], 0.0)
-        assert_equal(a[-1, :], 0.0)
+        assert_equal(a[:, 0], 0.0, strict=False)
+        assert_equal(a[:, -1], 0.0, strict=False)
+        assert_equal(a[0, :], 0.0, strict=False)
+        assert_equal(a[-1, :], 0.0, strict=False)
 
     @pytest.mark.parametrize(
         "dtype", [numpy.uint32, numpy.uint64] + get_all_dtypes(no_none=True)
@@ -426,7 +424,6 @@ def test_empty(self):
         expected = numpy.pad(a_np, [(2, 3), (3, 1)], "empty")
         result = dpnp.pad(a_dp, [(2, 3), (3, 1)], "empty")
         # omit uninitialized "empty" boundary from the comparison
-        assert result.shape == expected.shape
         assert_equal(result[2:-3, 3:-1], expected[2:-3, 3:-1])
 
     # Check how padding behaves on arrays with an empty dimension.

dpnp/tests/test_dlpack.py

@@ -40,7 +40,6 @@ def test_dtype_passthrough(self, xp, dt):
         x = xp.arange(5).astype(dt)
         y = xp.from_dlpack(x)
 
-        assert y.dtype == x.dtype
         assert_array_equal(x, y)
 
     @pytest.mark.parametrize("xp", [dpnp, numpy])

dpnp/tests/test_fill.py

@@ -38,7 +38,7 @@ def test_fill_strided_array():
     expected = dpnp.tile(dpnp.asarray([0, 1], dtype=a.dtype), 50)
 
     b.fill(1)
-    assert_array_equal(b, 1)
+    assert_array_equal(b, 1, strict=False)
     assert_array_equal(a, expected)
 
 
@@ -51,7 +51,7 @@ def test_fill_strided_2d_array(order):
     expected[::-2, ::2] = 1
 
     b.fill(1)
-    assert_array_equal(b, 1)
+    assert_array_equal(b, 1, strict=False)
     assert_array_equal(a, expected)
 
 
@@ -60,27 +60,27 @@ def test_fill_memset(order):
     a = dpnp.ones((10, 10), dtype="i4", order=order)
     a.fill(0)
 
-    assert_array_equal(a, 0)
+    assert_array_equal(a, 0, strict=False)
 
 
 def test_fill_float_complex_to_int():
     a = dpnp.ones((10, 10), dtype="i4")
 
     a.fill(complex(2, 0))
-    assert_array_equal(a, 2)
+    assert_array_equal(a, 2, strict=False)
 
     a.fill(float(3))
-    assert_array_equal(a, 3)
+    assert_array_equal(a, 3, strict=False)
 
 
 def test_fill_complex_to_float():
     a = dpnp.ones((10, 10), dtype="f4")
 
     a.fill(complex(2, 0))
-    assert_array_equal(a, 2)
+    assert_array_equal(a, 2, strict=False)
 
 
 def test_fill_bool():
     a = dpnp.full(5, fill_value=7, dtype="i4")
     a.fill(True)
-    assert_array_equal(a, 1)
+    assert_array_equal(a, 1, strict=False)