MULTIPLY enable broadcasting version 0.1 (#620)

* MULTIPLY enable broadcasting version 0.1

Author: densmirn

dpnp/backend/include/dpnp_gen_2arg_3type_tbl.hpp

@@ -66,7 +66,6 @@ MACRO_2ARG_3TYPES_OP(dpnp_fmod_c, cl::sycl::fmod((double)input_elem1, (double)in
 MACRO_2ARG_3TYPES_OP(dpnp_hypot_c, cl::sycl::hypot((double)input_elem1, (double)input_elem2), oneapi::mkl::vm::hypot)
 MACRO_2ARG_3TYPES_OP(dpnp_maximum_c, cl::sycl::max(input_elem1, input_elem2), oneapi::mkl::vm::fmax)
 MACRO_2ARG_3TYPES_OP(dpnp_minimum_c, cl::sycl::min(input_elem1, input_elem2), oneapi::mkl::vm::fmin)
-MACRO_2ARG_3TYPES_OP(dpnp_multiply_c, input_elem1* input_elem2, oneapi::mkl::vm::mul)
 MACRO_2ARG_3TYPES_OP(dpnp_power_c, cl::sycl::pow((double)input_elem1, (double)input_elem2), oneapi::mkl::vm::pow)
 MACRO_2ARG_3TYPES_OP(dpnp_subtract_c, input_elem1 - input_elem2, oneapi::mkl::vm::sub)
 

dpnp/backend/include/dpnp_iface.hpp

@@ -725,6 +725,33 @@ INP_DLLEXPORT void dpnp_floor_divide_c(void* array1_in, void* array2_in, void* r
 template <typename _DataType_input, typename _DataType_output>
 INP_DLLEXPORT void dpnp_modf_c(void* array1_in, void* result1_out, void* result2_out, size_t size);
 
+/**
+ * @ingroup BACKEND_API
+ * @brief multiply function.
+ *
+ * @param [out] result_out         Output array.
+ * @param [in]  input1_in          Input 1 either array or scalar.
+ * @param [in]  input1_size        Number of elements in input 1.
+ * @param [in]  input1_shape       Shape of input 1.
+ * @param [in]  input1_shape_ndim  Size of shape 1.
+ * @param [in]  input2_in          Input 2 either array or scalar.
+ * @param [in]  input2_size        Number of elements in input 2.
+ * @param [in]  input2_shape       Shape of input 2.
+ * @param [in]  input2_shape_ndim  Size of shape 2.
+ * @param [in]  where              Mask array.
+ */
+template <typename _DataType_output, typename _DataType_input1, typename _DataType_input2>
+INP_DLLEXPORT void dpnp_multiply_c(void* result_out,
+                                   const void* input1_in,
+                                   const size_t input1_size,
+                                   const size_t* input1_shape,
+                                   const size_t input1_shape_ndim,
+                                   const void* input2_in,
+                                   const size_t input2_size,
+                                   const size_t* input2_shape,
+                                   const size_t input2_shape_ndim,
+                                   const size_t* where);
+
 /**
  * @ingroup BACKEND_API
  * @brief Implementation of ones function

dpnp/backend/kernels/dpnp_krnl_elemwise.cpp

@@ -510,23 +510,6 @@ static void func_map_init_elemwise_2arg_3type(func_map_t& fmap)
     fmap[DPNPFuncName::DPNP_FN_MINIMUM][eft_DBL][eft_FLT] = {eft_DBL, (void*)dpnp_minimum_c<double, float, double>};
     fmap[DPNPFuncName::DPNP_FN_MINIMUM][eft_DBL][eft_DBL] = {eft_DBL, (void*)dpnp_minimum_c<double, double, double>};
 
-    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_INT][eft_INT] = {eft_INT, (void*)dpnp_multiply_c<int, int, int>};
-    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_INT][eft_LNG] = {eft_LNG, (void*)dpnp_multiply_c<int, long, long>};
-    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_INT][eft_FLT] = {eft_DBL, (void*)dpnp_multiply_c<int, float, double>};
-    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_INT][eft_DBL] = {eft_DBL, (void*)dpnp_multiply_c<int, double, double>};
-    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_LNG][eft_INT] = {eft_LNG, (void*)dpnp_multiply_c<long, int, long>};
-    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_LNG][eft_LNG] = {eft_LNG, (void*)dpnp_multiply_c<long, long, long>};
-    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_LNG][eft_FLT] = {eft_DBL, (void*)dpnp_multiply_c<long, float, double>};
-    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_LNG][eft_DBL] = {eft_DBL, (void*)dpnp_multiply_c<long, double, double>};
-    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_FLT][eft_INT] = {eft_DBL, (void*)dpnp_multiply_c<float, int, double>};
-    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_FLT][eft_LNG] = {eft_DBL, (void*)dpnp_multiply_c<float, long, double>};
-    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_FLT][eft_FLT] = {eft_FLT, (void*)dpnp_multiply_c<float, float, float>};
-    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_FLT][eft_DBL] = {eft_DBL, (void*)dpnp_multiply_c<float, double, double>};
-    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_DBL][eft_INT] = {eft_DBL, (void*)dpnp_multiply_c<double, int, double>};
-    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_DBL][eft_LNG] = {eft_DBL, (void*)dpnp_multiply_c<double, long, double>};
-    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_DBL][eft_FLT] = {eft_DBL, (void*)dpnp_multiply_c<double, float, double>};
-    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_DBL][eft_DBL] = {eft_DBL, (void*)dpnp_multiply_c<double, double, double>};
-
     fmap[DPNPFuncName::DPNP_FN_POWER][eft_INT][eft_INT] = {eft_INT, (void*)dpnp_power_c<int, int, int>};
     fmap[DPNPFuncName::DPNP_FN_POWER][eft_INT][eft_LNG] = {eft_LNG, (void*)dpnp_power_c<int, long, long>};
     fmap[DPNPFuncName::DPNP_FN_POWER][eft_INT][eft_FLT] = {eft_DBL, (void*)dpnp_power_c<int, float, double>};

dpnp/backend/kernels/dpnp_krnl_mathematical.cpp

@@ -232,6 +232,79 @@ void dpnp_modf_c(void* array1_in, void* result1_out, void* result2_out, size_t s
     event.wait();
 }
 
+template <typename _DataType_output, typename _DataType_input1, typename _DataType_input2>
+class dpnp_multiply_c_kernel;
+
+template <typename _DataType_output, typename _DataType_input1, typename _DataType_input2>
+void dpnp_multiply_c(void* result_out,
+                     const void* input1_in,
+                     const size_t input1_size,
+                     const size_t* input1_shape,
+                     const size_t input1_shape_ndim,
+                     const void* input2_in,
+                     const size_t input2_size,
+                     const size_t* input2_shape,
+                     const size_t input2_shape_ndim,
+                     const size_t* where)
+{
+    // avoid warning unused variable
+    (void)input1_shape;
+    (void)input1_shape_ndim;
+    (void)input2_shape;
+    (void)input2_shape_ndim;
+    (void)where;
+
+    if (!input1_size || !input2_size)
+    {
+        return;
+    }
+
+    const size_t result_size = (input2_size > input1_size) ? input2_size : input1_size;
+
+    const _DataType_input1* input1_data = reinterpret_cast<const _DataType_input1*>(input1_in);
+    const _DataType_input2* input2_data = reinterpret_cast<const _DataType_input2*>(input2_in);
+    _DataType_output* result = reinterpret_cast<_DataType_output*>(result_out);
+
+    cl::sycl::range<1> gws(result_size);
+    auto kernel_parallel_for_func = [=](cl::sycl::id<1> global_id) {
+        size_t i = global_id[0]; /*for (size_t i = 0; i < result_size; ++i)*/
+        {
+            const _DataType_input1 input1_elem = (input1_size == 1) ? input1_data[0] : input1_data[i];
+            const _DataType_input2 input2_elem = (input2_size == 1) ? input2_data[0] : input2_data[i];
+            result[i] = input1_elem * input2_elem;
+        }
+    };
+    auto kernel_func = [&](cl::sycl::handler& cgh) {
+        cgh.parallel_for<class dpnp_multiply_c_kernel<_DataType_output, _DataType_input1,
+                                                      _DataType_input2>>(gws, kernel_parallel_for_func);
+    };
+
+    cl::sycl::event event;
+
+    if (input1_size == input2_size)
+    {
+        if constexpr ((std::is_same<_DataType_input1, double>::value ||
+                       std::is_same<_DataType_input1, float>::value) &&
+                      std::is_same<_DataType_input2, _DataType_input1>::value)
+        {
+            _DataType_input1* input1 = const_cast<_DataType_input1*>(input1_data);
+            _DataType_input2* input2 = const_cast<_DataType_input2*>(input2_data);
+            // https://docs.oneapi.com/versions/latest/onemkl/mul.html
+            event = oneapi::mkl::vm::mul(DPNP_QUEUE, result_size, input1, input2, result);
+        }
+        else
+        {
+            event = DPNP_QUEUE.submit(kernel_func);
+        }
+    }
+    else
+    {
+        event = DPNP_QUEUE.submit(kernel_func);
+    }
+
+    event.wait();
+}
+
 template <typename _KernelNameSpecialization1, typename _KernelNameSpecialization2, typename _KernelNameSpecialization3>
 class dpnp_remainder_c_kernel;
 
@@ -411,6 +484,34 @@ void func_map_init_mathematical(func_map_t& fmap)
     fmap[DPNPFuncName::DPNP_FN_MODF][eft_FLT][eft_FLT] = {eft_FLT, (void*)dpnp_modf_c<float, float>};
     fmap[DPNPFuncName::DPNP_FN_MODF][eft_DBL][eft_DBL] = {eft_DBL, (void*)dpnp_modf_c<double, double>};
 
+    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_BLN][eft_BLN] = {eft_BLN, (void*)dpnp_multiply_c<bool, bool, bool>};
+    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_BLN][eft_INT] = {eft_INT, (void*)dpnp_multiply_c<int, bool, int>};
+    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_BLN][eft_LNG] = {eft_LNG, (void*)dpnp_multiply_c<long, bool, long>};
+    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_BLN][eft_FLT] = {eft_FLT, (void*)dpnp_multiply_c<float, bool, float>};
+    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_BLN][eft_DBL] = {eft_DBL, (void*)dpnp_multiply_c<double, bool, double>};
+    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_INT][eft_BLN] = {eft_INT, (void*)dpnp_multiply_c<int, int, bool>};
+    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_INT][eft_INT] = {eft_INT, (void*)dpnp_multiply_c<int, int, int>};
+    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_INT][eft_LNG] = {eft_LNG, (void*)dpnp_multiply_c<long, int, long>};
+    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_INT][eft_FLT] = {eft_DBL, (void*)dpnp_multiply_c<double, int, float>};
+    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_INT][eft_DBL] = {eft_DBL, (void*)dpnp_multiply_c<double, int, double>};
+    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_LNG][eft_BLN] = {eft_LNG, (void*)dpnp_multiply_c<long, long, bool>};
+    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_LNG][eft_INT] = {eft_LNG, (void*)dpnp_multiply_c<long, long, int>};
+    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_LNG][eft_LNG] = {eft_LNG, (void*)dpnp_multiply_c<long, long, long>};
+    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_LNG][eft_FLT] = {eft_DBL, (void*)dpnp_multiply_c<double, long, float>};
+    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_LNG][eft_DBL] = {eft_DBL, (void*)dpnp_multiply_c<double, long, double>};
+    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_FLT][eft_BLN] = {eft_FLT, (void*)dpnp_multiply_c<float, float, bool>};
+    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_FLT][eft_INT] = {eft_DBL, (void*)dpnp_multiply_c<double, float, int>};
+    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_FLT][eft_LNG] = {eft_DBL, (void*)dpnp_multiply_c<double, float, long>};
+    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_FLT][eft_FLT] = {eft_FLT, (void*)dpnp_multiply_c<float, float, float>};
+    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_FLT][eft_DBL] = {eft_DBL, (void*)dpnp_multiply_c<double, float, double>};
+    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_DBL][eft_BLN] = {eft_DBL, (void*)dpnp_multiply_c<double, double, bool>};
+    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_DBL][eft_INT] = {eft_DBL, (void*)dpnp_multiply_c<double, double, int>};
+    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_DBL][eft_LNG] = {eft_DBL, (void*)dpnp_multiply_c<double, double, long>};
+    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_DBL][eft_FLT] = {eft_DBL, (void*)dpnp_multiply_c<double, double, float>};
+    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_DBL][eft_DBL] = {eft_DBL, (void*)dpnp_multiply_c<double, double, double>};
+    fmap[DPNPFuncName::DPNP_FN_MULTIPLY][eft_C128][eft_C128] = {
+        eft_C128, (void*)dpnp_multiply_c<std::complex<double>, std::complex<double>, std::complex<double>>};
+
     fmap[DPNPFuncName::DPNP_FN_REMAINDER][eft_INT][eft_INT] = {eft_INT, (void*)dpnp_remainder_c<int, int, int>};
     fmap[DPNPFuncName::DPNP_FN_REMAINDER][eft_INT][eft_LNG] = {eft_LNG, (void*)dpnp_remainder_c<int, long, long>};
     fmap[DPNPFuncName::DPNP_FN_REMAINDER][eft_INT][eft_FLT] = {eft_DBL, (void*)dpnp_remainder_c<int, float, double>};

dpnp/dpnp_algo/dpnp_algo.pxd

@@ -212,6 +212,8 @@ ctypedef void(*fptr_1out_t)(void *, size_t)
 ctypedef void(*fptr_1in_1out_t)(void * , void * , size_t)
 ctypedef void(*fptr_2in_1out_t)(void * , void*, void*, size_t)
 ctypedef void(*fptr_2in_1out_new_t)(void * , void*, size_t, void*, size_t) # to be fused with fptr_2in_1out_t
+ctypedef void(*fptr_2in_1out_full_t)(void *, const void *, const size_t, const long*, const size_t,
+                                     const void *, const size_t, const long*, const size_t, const long*)
 ctypedef void(*fptr_blas_gemm_2in_1out_t)(void * , void * , void * , size_t, size_t, size_t)
 ctypedef void(*dpnp_reduction_c_t)(void * , const void * , const size_t*, const size_t, const long*, const size_t, const void * , const long*)
 
@@ -282,7 +284,7 @@ cpdef dparray dpnp_divide(dparray array1, dparray array2)
 cpdef dparray dpnp_hypot(dparray array1, dparray array2)
 cpdef dparray dpnp_maximum(dparray array1, dparray array2)
 cpdef dparray dpnp_minimum(dparray array1, dparray array2)
-cpdef dparray dpnp_multiply(dparray array1, array2)
+cpdef dparray dpnp_multiply(object x1_obj, object x2_obj)
 cpdef dparray dpnp_negative(dparray array1)
 cpdef dparray dpnp_power(dparray array1, array2)
 cpdef dparray dpnp_remainder(dparray array1, dparray array2)

dpnp/dpnp_algo/dpnp_algo_mathematical.pyx

@@ -32,7 +32,6 @@ and the rest of the library
 
 """
 
-
 from dpnp.dpnp_utils cimport *
 import dpnp
 import numpy
@@ -273,23 +272,50 @@ cpdef tuple dpnp_modf(dparray x1):
     return result1, result2
 
 
-cpdef dparray dpnp_multiply(dparray x1, x2):
-    cdef dparray result
-    if dpnp.isscalar(x2):
-        x2_ = dpnp.array([x2])
+cpdef dparray dpnp_multiply(object x1_obj, object x2_obj, dparray out=None, object where=True):
+    cdef dparray_shape_type x1_shape, x2_shape, result_shape
 
-        types_map = {
-            (dpnp.int32, dpnp.float64): dpnp.float64,
-            (dpnp.int64, dpnp.float64): dpnp.float64,
-        }
+    cdef bint x1_obj_is_dparray = isinstance(x1_obj, dparray)
+    cdef bint x2_obj_is_dparray = isinstance(x2_obj, dparray)
 
-        res_type = types_map.get((x1.dtype.type, x2_.dtype.type), x1.dtype)
-        result = dparray(x1.shape, dtype=res_type)
-        for i in range(x1.size):
-            result[i] = x1[i] * x2
-        return result.reshape(x1.shape)
+    cdef dparray x1_dparray, x2_dparray
+
+    common_type = find_common_type(x1_obj, x2_obj)
+
+    if x1_obj_is_dparray:
+        x1_dparray = x1_obj
     else:
-        return call_fptr_2in_1out(DPNP_FN_MULTIPLY, x1, x2, x1.shape)
+        x1_dparray = dparray((1,), dtype=common_type)
+        copy_values_to_dparray(x1_dparray, (x1_obj,))
+
+    if x2_obj_is_dparray:
+        x2_dparray = x2_obj
+    else:
+        x2_dparray = dparray((1,), dtype=common_type)
+        copy_values_to_dparray(x2_dparray, (x2_obj,))
+
+    x1_shape = x1_dparray.shape
+    x2_shape = x2_dparray.shape
+    result_shape = get_common_shape(x1_shape, x2_shape)
+
+    # Convert string type names (dparray.dtype) to C enum DPNPFuncType
+    cdef DPNPFuncType x1_c_type = dpnp_dtype_to_DPNPFuncType(x1_dparray.dtype)
+    cdef DPNPFuncType x2_c_type = dpnp_dtype_to_DPNPFuncType(x2_dparray.dtype)
+
+    # get the FPTR data structure
+    cdef DPNPFuncData kernel_data = get_dpnp_function_ptr(DPNP_FN_MULTIPLY, x1_c_type, x2_c_type)
+
+    cdef DPNPFuncType result_c_type = get_output_c_type(DPNP_FN_MULTIPLY, kernel_data.return_type, out, None)
+
+    # Create result array
+    cdef dparray result = create_output_array(result_shape, result_c_type, out)
+
+    cdef fptr_2in_1out_full_t func = <fptr_2in_1out_full_t > kernel_data.ptr
+    # Call FPTR function
+    func(result.get_data(), x1_dparray.get_data(), x1_dparray.size, x1_shape.data(), x1_shape.size(),
+         x2_dparray.get_data(), x2_dparray.size, x2_shape.data(), x2_shape.size(), NULL)
+
+    return result
 
 
 cpdef dparray dpnp_nancumprod(dparray x1):

dpnp/dpnp_iface_mathematical.py

@@ -972,7 +972,7 @@ def modf(x, **kwargs):
     return call_origin(numpy.modf, x, **kwargs)
 
 
-def multiply(x1, x2, **kwargs):
+def multiply(x1, x2, out=None, where=True, **kwargs):
     """
     Multiply arguments element-wise.
 
@@ -995,33 +995,32 @@ def multiply(x1, x2, **kwargs):
     [1, 4, 9, 16, 25]
 
     """
+    x1_is_scalar, x2_is_scalar = dpnp.isscalar(x1), dpnp.isscalar(x2)
+    x1_is_dparray, x2_is_dparray = isinstance(x1, dparray), isinstance(x2, dparray)
 
-    is_x1_dparray = isinstance(x1, dparray)
-    is_x2_dparray = isinstance(x2, dparray)
-
-    is_x1_scalar = dpnp.isscalar(x1)
-    is_x2_scalar = dpnp.isscalar(x2)
-
-    if not use_origin_backend(x1):
-        if kwargs:
+    if not use_origin_backend(x1) and not kwargs:
+        if not x1_is_dparray and not x1_is_scalar:
+            pass
+        elif not x2_is_dparray and not x2_is_scalar:
+            pass
+        elif x1_is_scalar and x2_is_scalar:
             pass
-        elif not (is_x1_dparray or is_x1_scalar):
+        elif x1_is_dparray and x1.ndim == 0:
             pass
-        elif not (is_x2_dparray or is_x2_scalar):
+        elif x2_is_dparray and x2.ndim == 0:
             pass
-        elif is_x1_scalar and is_x2_scalar:
+        elif x1_is_dparray and x2_is_dparray and x1.size != x2.size:
             pass
-        elif (is_x1_dparray and is_x2_dparray) and (x1.size != x2.size):
+        elif x1_is_dparray and x2_is_dparray and x1.shape != x2.shape:
+            pass
+        elif out is not None and not isinstance(out, dparray):
             pass
-        elif (is_x1_dparray and is_x2_dparray) and (x1.shape != x2.shape):
+        elif not where:
             pass
         else:
-            if is_x1_scalar:
-                return dpnp_multiply(x2, x1)
-            else:
-                return dpnp_multiply(x1, x2)
+            return dpnp_multiply(x1, x2, out, where)
 
-    return call_origin(numpy.multiply, x1, x2, **kwargs)
+    return call_origin(numpy.multiply, x1, x2, out=out, where=where, **kwargs)
 
 
 def nancumprod(x1, **kwargs):

dpnp/dpnp_utils/dpnp_algo_utils.pxd

@@ -95,6 +95,11 @@ Returns a tuple of:
 2. dtype
 """
 
+cpdef find_common_type(object x1_obj, object x2_obj)
+"""
+Find common type of 2 input objects
+"""
+
 cdef long copy_values_to_dparray(dparray dst, input_obj, size_t dst_idx=*) except -1
 """
 Copy values to `dst` by iterating element by element in `input_obj`
@@ -125,6 +130,11 @@ cpdef nd2dp_array(arr)
 Convert ndarray to dparray
 """
 
+cdef dparray_shape_type get_common_shape(dparray_shape_type input1_shape, dparray_shape_type input2_shape)
+"""
+Calculate common shape from input shapes
+"""
+
 cdef dparray_shape_type get_reduction_output_shape(dparray_shape_type input_shape, object axis, cpp_bool keepdims)
 """
 Calculate output array shape in reduction functions