impelement dpnp.norm (#1746)

* impelement dpnp.norm

* address comments

* add float ty description

* improve test coverage

* make axis for test_norm_ND and ND_complex similar

* mute some tests for on windows

* unmute test on windows

---------

Co-authored-by: Anton <100830759+antonwolfy@users.noreply.github.com>

Author: vtavana

dpnp/backend/extensions/blas/blas_py.cpp

@@ -52,15 +52,15 @@ PYBIND11_MODULE(_blas_impl, m)
 
     {
         m.def("_dot", &blas_ext::dot,
-              "Call `dot` from OneMKL LAPACK library to return "
+              "Call `dot` from OneMKL BLAS library to return "
               "the dot product of two real-valued vectors.",
               py::arg("sycl_queue"), py::arg("vectorA"), py::arg("vectorB"),
               py::arg("result"), py::arg("depends") = py::list());
     }
 
     {
         m.def("_dotc", &blas_ext::dotc,
-              "Call `dotc` from OneMKL LAPACK library to return "
+              "Call `dotc` from OneMKL BLAS library to return "
               "the dot product of two complex vectors, "
               "conjugating the first vector.",
               py::arg("sycl_queue"), py::arg("vectorA"), py::arg("vectorB"),
@@ -69,23 +69,23 @@ PYBIND11_MODULE(_blas_impl, m)
 
     {
         m.def("_dotu", &blas_ext::dotu,
-              "Call `dotu` from OneMKL LAPACK library to return "
+              "Call `dotu` from OneMKL BLAS library to return "
               "the dot product of two complex vectors.",
               py::arg("sycl_queue"), py::arg("vectorA"), py::arg("vectorB"),
               py::arg("result"), py::arg("depends") = py::list());
     }
 
     {
         m.def("_gemm", &blas_ext::gemm,
-              "Call `gemm` from OneMKL LAPACK library to return "
+              "Call `gemm` from OneMKL BLAS library to return "
               "the matrix-matrix product with 2-D matrices.",
               py::arg("sycl_queue"), py::arg("matrixA"), py::arg("matrixB"),
               py::arg("result"), py::arg("depends") = py::list());
     }
 
     {
         m.def("_gemm_batch", &blas_ext::gemm_batch,
-              "Call `gemm_batch` from OneMKL LAPACK library to return "
+              "Call `gemm_batch` from OneMKL BLAS library to return "
               "the matrix-matrix product for a batch of 2-D matrices.",
               py::arg("sycl_queue"), py::arg("matrixA"), py::arg("matrixB"),
               py::arg("result"), py::arg("batch_size"), py::arg("stridea"),

dpnp/linalg/dpnp_algo_linalg.pyx

@@ -48,7 +48,6 @@ __all__ = [
     "dpnp_cond",
     "dpnp_eig",
     "dpnp_eigvals",
-    "dpnp_norm",
 ]
 
 
@@ -171,108 +170,3 @@ cpdef utils.dpnp_descriptor dpnp_eigvals(utils.dpnp_descriptor input):
     c_dpctl.DPCTLEvent_Delete(event_ref)
 
     return res_val
-
-
-cpdef object dpnp_norm(object input, ord=None, axis=None):
-    cdef long size_input = input.size
-    cdef shape_type_c shape_input = input.shape
-
-    dev = input.get_array().sycl_device
-    if input.dtype == dpnp.float32 or not dev.has_aspect_fp64:
-        res_type = dpnp.float32
-    else:
-        res_type = dpnp.float64
-
-    if size_input == 0:
-        return dpnp.array([dpnp.nan], dtype=res_type)
-
-    if isinstance(axis, int):
-        axis_ = tuple([axis])
-    else:
-        axis_ = axis
-
-    ndim = input.ndim
-    if axis is None:
-        if ((ord is None) or
-            (ord in ('f', 'fro') and ndim == 2) or
-                (ord == 2 and ndim == 1)):
-
-            # TODO: change order='K' when support is implemented
-            input = dpnp.ravel(input, order='C')
-            sqnorm = dpnp.dot(input, input)
-            ret = dpnp.sqrt([sqnorm], dtype=res_type)
-            return dpnp.array(ret.reshape(1, *ret.shape), dtype=res_type)
-
-    len_axis = 1 if axis is None else len(axis_)
-    if len_axis == 1:
-        if ord == dpnp.inf:
-            return dpnp.array([dpnp.abs(input).max(axis=axis)])
-        elif ord == -dpnp.inf:
-            return dpnp.array([dpnp.abs(input).min(axis=axis)])
-        elif ord == 0:
-            return input.dtype.type(dpnp.count_nonzero(input, axis=axis))
-        elif ord is None or ord == 2:
-            s = input * input
-            return dpnp.sqrt(dpnp.sum(s, axis=axis), dtype=res_type)
-        elif isinstance(ord, str):
-            raise ValueError(f"Invalid norm order '{ord}' for vectors")
-        else:
-            absx = dpnp.abs(input)
-            absx_size = absx.size
-            absx_power = utils_py.create_output_descriptor_py((absx_size,), absx.dtype, None).get_pyobj()
-
-            absx_flatiter = absx.flat
-
-            for i in range(absx_size):
-                absx_elem = absx_flatiter[i]
-                absx_power[i] = absx_elem ** ord
-            absx_ = dpnp.reshape(absx_power, absx.shape)
-            ret = dpnp.sum(absx_, axis=axis)
-            ret_size = ret.size
-            ret_power = utils_py.create_output_descriptor_py((ret_size,), None, None).get_pyobj()
-
-            ret_flatiter = ret.flat
-
-            for i in range(ret_size):
-                ret_elem = ret_flatiter[i]
-                ret_power[i] = ret_elem ** (1 / ord)
-            ret_ = dpnp.reshape(ret_power, ret.shape)
-            return ret_
-    elif len_axis == 2:
-        row_axis, col_axis = axis_
-        if row_axis == col_axis:
-            raise ValueError('Duplicate axes given.')
-        # if ord == 2:
-        #     ret =  _multi_svd_norm(input, row_axis, col_axis, amax)
-        # elif ord == -2:
-        #     ret = _multi_svd_norm(input, row_axis, col_axis, amin)
-        elif ord == 1:
-            if col_axis > row_axis:
-                col_axis -= 1
-            dpnp_sum_val = dpnp.sum(dpnp.abs(input), axis=row_axis)
-            ret = dpnp_sum_val.min(axis=col_axis)
-        elif ord == dpnp.inf:
-            if row_axis > col_axis:
-                row_axis -= 1
-            dpnp_sum_val = dpnp.sum(dpnp.abs(input), axis=col_axis)
-            ret = dpnp_sum_val.max(axis=row_axis)
-        elif ord == -1:
-            if col_axis > row_axis:
-                col_axis -= 1
-            dpnp_sum_val = dpnp.sum(dpnp.abs(input), axis=row_axis)
-            ret = dpnp_sum_val.min(axis=col_axis)
-        elif ord == -dpnp.inf:
-            if row_axis > col_axis:
-                row_axis -= 1
-            dpnp_sum_val = dpnp.sum(dpnp.abs(input), axis=col_axis)
-            ret = dpnp_sum_val.min(axis=row_axis)
-        elif ord in [None, 'fro', 'f']:
-            ret = dpnp.sqrt(dpnp.sum(input * input, axis=axis))
-        # elif ord == 'nuc':
-        #     ret = _multi_svd_norm(input, row_axis, col_axis, sum)
-        else:
-            raise ValueError("Invalid norm order for matrices.")
-
-        return ret
-    else:
-        raise ValueError("Improper number of dimensions to norm.")

dpnp/linalg/dpnp_iface_linalg.py

@@ -54,6 +54,7 @@
     dpnp_matrix_power,
     dpnp_matrix_rank,
     dpnp_multi_dot,
+    dpnp_norm,
     dpnp_pinv,
     dpnp_qr,
     dpnp_slogdet,
@@ -491,7 +492,7 @@ def multi_dot(arrays, *, out=None):
     """
     Compute the dot product of two or more arrays in a single function call.
 
-    For full documentation refer to :obj:`numpy.multi_dot`.
+    For full documentation refer to :obj:`numpy.linalg.multi_dot`.
 
     Parameters
     ----------
@@ -602,60 +603,109 @@ def pinv(a, rcond=1e-15, hermitian=False):
     return dpnp_pinv(a, rcond=rcond, hermitian=hermitian)
 
 
-def norm(x1, ord=None, axis=None, keepdims=False):
+def norm(x, ord=None, axis=None, keepdims=False):
     """
     Matrix or vector norm.
 
-    This function is able to return one of eight different matrix norms,
-    or one of an infinite number of vector norms (described below), depending
-    on the value of the ``ord`` parameter.
+    For full documentation refer to :obj:`numpy.linalg.norm`.
 
     Parameters
     ----------
-    input : array_like
-        Input array.  If `axis` is None, `x` must be 1-D or 2-D, unless `ord`
-        is None. If both `axis` and `ord` are None, the 2-norm of
-        ``x.ravel`` will be returned.
-    ord : optional
-        Order of the norm (see table under ``Notes``). inf means numpy's
-        `inf` object. The default is None.
-    axis : optional.
+    x : {dpnp.ndarray, usm_ndarray}
+        Input array.  If `axis` is ``None``, `x` must be 1-D or 2-D, unless
+        `ord` is ``None``. If both `axis` and `ord` are ``None``, the 2-norm
+        of ``x.ravel`` will be returned.
+    ord : {int, float, inf, -inf, "fro", "nuc"}, optional
+        Norm type. inf means dpnp's `inf` object. The default is ``None``.
+    axis : {None, int, 2-tuple of ints}, optional
         If `axis` is an integer, it specifies the axis of `x` along which to
         compute the vector norms.  If `axis` is a 2-tuple, it specifies the
         axes that hold 2-D matrices, and the matrix norms of these matrices
-        are computed.  If `axis` is None then either a vector norm (when `x`
-        is 1-D) or a matrix norm (when `x` is 2-D) is returned. The default
-        is None.
+        are computed.  If `axis` is ``None`` then either a vector norm (when
+        `x` is 1-D) or a matrix norm (when `x` is 2-D) is returned.
+        The default is ``None``.
     keepdims : bool, optional
-        If this is set to True, the axes which are normed over are left in the
-        result as dimensions with size one.  With this option the result will
-        broadcast correctly against the original `x`.
+        If this is set to ``True``, the axes which are normed over are left in
+        the result as dimensions with size one.  With this option the result
+        will broadcast correctly against the original `x`.
 
     Returns
     -------
-    n : float or ndarray
+    out : dpnp.ndarray
         Norm of the matrix or vector(s).
-    """
 
-    x1_desc = dpnp.get_dpnp_descriptor(x1, copy_when_nondefault_queue=False)
-    if x1_desc:
-        if (
-            not isinstance(axis, int)
-            and not isinstance(axis, tuple)
-            and axis is not None
-        ):
-            pass
-        elif keepdims is not False:
-            pass
-        elif ord not in [None, 0, 3, "fro", "f"]:
-            pass
-        else:
-            result_obj = dpnp_norm(x1, ord=ord, axis=axis)
-            result = dpnp.convert_single_elem_array_to_scalar(result_obj)
+    Examples
+    --------
+    >>> import dpnp as np
+    >>> a = np.arange(9) - 4
+    >>> a
+    array([-4, -3, -2, -1,  0,  1,  2,  3,  4])
+    >>> b = a.reshape((3, 3))
+    >>> b
+    array([[-4, -3, -2],
+           [-1,  0,  1],
+           [ 2,  3,  4]])
+
+    >>> np.linalg.norm(a)
+    array(7.745966692414834)
+    >>> np.linalg.norm(b)
+    array(7.745966692414834)
+    >>> np.linalg.norm(b, 'fro')
+    array(7.745966692414834)
+    >>> np.linalg.norm(a, np.inf)
+    array(4.)
+    >>> np.linalg.norm(b, np.inf)
+    array(9.)
+    >>> np.linalg.norm(a, -np.inf)
+    array(0.)
+    >>> np.linalg.norm(b, -np.inf)
+    array(2.)
+
+    >>> np.linalg.norm(a, 1)
+    array(20.)
+    >>> np.linalg.norm(b, 1)
+    array(7.)
+    >>> np.linalg.norm(a, -1)
+    array(0.)
+    >>> np.linalg.norm(b, -1)
+    array(6.)
+    >>> np.linalg.norm(a, 2)
+    array(7.745966692414834)
+    >>> np.linalg.norm(b, 2)
+    array(7.3484692283495345)
+
+    >>> np.linalg.norm(a, -2)
+    array(0.)
+    >>> np.linalg.norm(b, -2)
+    array(1.8570331885190563e-016) # may vary
+    >>> np.linalg.norm(a, 3)
+    array(5.8480354764257312) # may vary
+    >>> np.linalg.norm(a, -3)
+    array(0.)
+
+    Using the `axis` argument to compute vector norms:
+
+    >>> c = np.array([[ 1, 2, 3],
+    ...               [-1, 1, 4]])
+    >>> np.linalg.norm(c, axis=0)
+    array([ 1.41421356,  2.23606798,  5.        ])
+    >>> np.linalg.norm(c, axis=1)
+    array([ 3.74165739,  4.24264069])
+    >>> np.linalg.norm(c, ord=1, axis=1)
+    array([ 6.,  6.])
+
+    Using the `axis` argument to compute matrix norms:
+
+    >>> m = np.arange(8).reshape(2,2,2)
+    >>> np.linalg.norm(m, axis=(1,2))
+    array([  3.74165739,  11.22497216])
+    >>> np.linalg.norm(m[0, :, :]), np.linalg.norm(m[1, :, :])
+    (array(3.7416573867739413), array(11.224972160321824))
 
-            return result
+    """
 
-    return call_origin(numpy.linalg.norm, x1, ord, axis, keepdims)
+    dpnp.check_supported_arrays_type(x)
+    return dpnp_norm(x, ord, axis, keepdims)
 
 
 def qr(a, mode="reduced"):