Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

impelement dpnp.norm #1746

Merged
merged 16 commits into from
Apr 3, 2024
Merged
Show file tree
Hide file tree
Changes from 2 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 5 additions & 5 deletions dpnp/backend/extensions/blas/blas_py.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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"),
Expand All @@ -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"),
Expand Down
106 changes: 0 additions & 106 deletions dpnp/linalg/dpnp_algo_linalg.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,6 @@ __all__ = [
"dpnp_cond",
"dpnp_eig",
"dpnp_eigvals",
"dpnp_norm",
]


Expand Down Expand Up @@ -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.")
130 changes: 90 additions & 40 deletions dpnp/linalg/dpnp_iface_linalg.py
Original file line number Diff line number Diff line change
Expand Up @@ -53,6 +53,7 @@
dpnp_inv,
dpnp_matrix_rank,
dpnp_multi_dot,
dpnp_norm,
dpnp_pinv,
dpnp_qr,
dpnp_slogdet,
Expand Down Expand Up @@ -456,7 +457,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
----------
Expand Down Expand Up @@ -567,60 +568,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 : {non-zero int, inf, -inf, "fro", "nuc"}, optional
vtavana marked this conversation as resolved.
Show resolved Hide resolved
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)

return result

return call_origin(numpy.linalg.norm, x1, ord, axis, keepdims)
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))

"""

dpnp.check_supported_arrays_type(x)
return dpnp_norm(x, ord, axis, keepdims)


def qr(a, mode="reduced"):
Expand Down
Loading
Loading