Skip to content

Commit

Permalink
Update dpnp.linalg.matrix_rank() implementation (IntelPython#1717)
Browse files Browse the repository at this point in the history
* Update dpnp.linalg.matrix_rank impl

* Add cupy tests for dpnp.linalg.matrix_rank

* Add dpnp tests for dpnp.linalg.matrix_rank

* Remove old impl of dpnp_matrix_rank

* Address remarks

* Address remarks
  • Loading branch information
vlad-perevezentsev authored Feb 19, 2024
1 parent a100705 commit d34b8ce
Show file tree
Hide file tree
Showing 10 changed files with 253 additions and 142 deletions.
42 changes: 20 additions & 22 deletions dpnp/backend/include/dpnp_iface_fptr.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -178,28 +178,26 @@ enum class DPNPFuncName : size_t
DPNP_FN_KRON, /**< Used in numpy.kron() impl */
DPNP_FN_KRON_EXT, /**< Used in numpy.kron() impl, requires extra parameters
*/
DPNP_FN_LEFT_SHIFT, /**< Used in numpy.left_shift() impl */
DPNP_FN_LOG, /**< Used in numpy.log() impl */
DPNP_FN_LOG10, /**< Used in numpy.log10() impl */
DPNP_FN_LOG2, /**< Used in numpy.log2() impl */
DPNP_FN_LOG1P, /**< Used in numpy.log1p() impl */
DPNP_FN_MATMUL, /**< Used in numpy.matmul() impl */
DPNP_FN_MATRIX_RANK, /**< Used in numpy.linalg.matrix_rank() impl */
DPNP_FN_MATRIX_RANK_EXT, /**< Used in numpy.linalg.matrix_rank() impl,
requires extra parameters */
DPNP_FN_MAX, /**< Used in numpy.max() impl */
DPNP_FN_MAXIMUM, /**< Used in numpy.fmax() impl */
DPNP_FN_MAXIMUM_EXT, /**< Used in numpy.fmax() impl , requires extra
parameters */
DPNP_FN_MEAN, /**< Used in numpy.mean() impl */
DPNP_FN_MEDIAN, /**< Used in numpy.median() impl */
DPNP_FN_MEDIAN_EXT, /**< Used in numpy.median() impl, requires extra
parameters */
DPNP_FN_MIN, /**< Used in numpy.min() impl */
DPNP_FN_MINIMUM, /**< Used in numpy.fmin() impl */
DPNP_FN_MINIMUM_EXT, /**< Used in numpy.fmax() impl, requires extra
parameters */
DPNP_FN_MODF, /**< Used in numpy.modf() impl */
DPNP_FN_LEFT_SHIFT, /**< Used in numpy.left_shift() impl */
DPNP_FN_LOG, /**< Used in numpy.log() impl */
DPNP_FN_LOG10, /**< Used in numpy.log10() impl */
DPNP_FN_LOG2, /**< Used in numpy.log2() impl */
DPNP_FN_LOG1P, /**< Used in numpy.log1p() impl */
DPNP_FN_MATMUL, /**< Used in numpy.matmul() impl */
DPNP_FN_MATRIX_RANK, /**< Used in numpy.linalg.matrix_rank() impl */
DPNP_FN_MAX, /**< Used in numpy.max() impl */
DPNP_FN_MAXIMUM, /**< Used in numpy.fmax() impl */
DPNP_FN_MAXIMUM_EXT, /**< Used in numpy.fmax() impl , requires extra
parameters */
DPNP_FN_MEAN, /**< Used in numpy.mean() impl */
DPNP_FN_MEDIAN, /**< Used in numpy.median() impl */
DPNP_FN_MEDIAN_EXT, /**< Used in numpy.median() impl, requires extra
parameters */
DPNP_FN_MIN, /**< Used in numpy.min() impl */
DPNP_FN_MINIMUM, /**< Used in numpy.fmin() impl */
DPNP_FN_MINIMUM_EXT, /**< Used in numpy.fmax() impl, requires extra
parameters */
DPNP_FN_MODF, /**< Used in numpy.modf() impl */
DPNP_FN_MODF_EXT, /**< Used in numpy.modf() impl, requires extra parameters
*/
DPNP_FN_MULTIPLY, /**< Used in numpy.multiply() impl */
Expand Down
18 changes: 0 additions & 18 deletions dpnp/backend/kernels/dpnp_krnl_linalg.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -579,15 +579,6 @@ template <typename _DataType>
void (*dpnp_matrix_rank_default_c)(void *, void *, shape_elem_type *, size_t) =
dpnp_matrix_rank_c<_DataType>;

template <typename _DataType>
DPCTLSyclEventRef (*dpnp_matrix_rank_ext_c)(DPCTLSyclQueueRef,
void *,
void *,
shape_elem_type *,
size_t,
const DPCTLEventVectorRef) =
dpnp_matrix_rank_c<_DataType>;

template <typename _InputDT, typename _ComputeDT>
DPCTLSyclEventRef dpnp_qr_c(DPCTLSyclQueueRef q_ref,
void *array1_in,
Expand Down Expand Up @@ -969,15 +960,6 @@ void func_map_init_linalg_func(func_map_t &fmap)
fmap[DPNPFuncName::DPNP_FN_MATRIX_RANK][eft_DBL][eft_DBL] = {
eft_DBL, (void *)dpnp_matrix_rank_default_c<double>};

fmap[DPNPFuncName::DPNP_FN_MATRIX_RANK_EXT][eft_INT][eft_INT] = {
eft_INT, (void *)dpnp_matrix_rank_ext_c<int32_t>};
fmap[DPNPFuncName::DPNP_FN_MATRIX_RANK_EXT][eft_LNG][eft_LNG] = {
eft_LNG, (void *)dpnp_matrix_rank_ext_c<int64_t>};
fmap[DPNPFuncName::DPNP_FN_MATRIX_RANK_EXT][eft_FLT][eft_FLT] = {
eft_FLT, (void *)dpnp_matrix_rank_ext_c<float>};
fmap[DPNPFuncName::DPNP_FN_MATRIX_RANK_EXT][eft_DBL][eft_DBL] = {
eft_DBL, (void *)dpnp_matrix_rank_ext_c<double>};

fmap[DPNPFuncName::DPNP_FN_QR][eft_INT][eft_INT] = {
eft_DBL, (void *)dpnp_qr_default_c<int32_t, double>};
fmap[DPNPFuncName::DPNP_FN_QR][eft_LNG][eft_LNG] = {
Expand Down
2 changes: 0 additions & 2 deletions dpnp/dpnp_algo/dpnp_algo.pxd
Original file line number Diff line number Diff line change
Expand Up @@ -78,8 +78,6 @@ cdef extern from "dpnp_iface_fptr.hpp" namespace "DPNPFuncName": # need this na
DPNP_FN_FULL_LIKE
DPNP_FN_KRON
DPNP_FN_KRON_EXT
DPNP_FN_MATRIX_RANK
DPNP_FN_MATRIX_RANK_EXT
DPNP_FN_MAXIMUM
DPNP_FN_MAXIMUM_EXT
DPNP_FN_MEDIAN
Expand Down
46 changes: 0 additions & 46 deletions dpnp/linalg/dpnp_algo_linalg.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -48,24 +48,14 @@ __all__ = [
"dpnp_cond",
"dpnp_eig",
"dpnp_eigvals",
"dpnp_matrix_rank",
"dpnp_norm",
]


# C function pointer to the C library template functions
ctypedef c_dpctl.DPCTLSyclEventRef(*custom_linalg_1in_1out_func_ptr_t)(c_dpctl.DPCTLSyclQueueRef,
void *, void * ,shape_elem_type * ,
size_t, const c_dpctl.DPCTLEventVectorRef)
ctypedef c_dpctl.DPCTLSyclEventRef(*custom_linalg_1in_1out_func_ptr_t_)(c_dpctl.DPCTLSyclQueueRef,
void * , void * , size_t * ,
const c_dpctl.DPCTLEventVectorRef)
ctypedef c_dpctl.DPCTLSyclEventRef(*custom_linalg_1in_1out_with_size_func_ptr_t_)(c_dpctl.DPCTLSyclQueueRef,
void *, void * , size_t,
const c_dpctl.DPCTLEventVectorRef)
ctypedef c_dpctl.DPCTLSyclEventRef(*custom_linalg_1in_3out_shape_t)(c_dpctl.DPCTLSyclQueueRef,
void *, void * , void * , void * ,
size_t , size_t, const c_dpctl.DPCTLEventVectorRef)
ctypedef c_dpctl.DPCTLSyclEventRef(*custom_linalg_2in_1out_func_ptr_t)(c_dpctl.DPCTLSyclQueueRef,
void *, void * , void * , size_t,
const c_dpctl.DPCTLEventVectorRef)
Expand Down Expand Up @@ -183,42 +173,6 @@ cpdef utils.dpnp_descriptor dpnp_eigvals(utils.dpnp_descriptor input):
return res_val


cpdef utils.dpnp_descriptor dpnp_matrix_rank(utils.dpnp_descriptor input):
cdef shape_type_c input_shape = input.shape
cdef DPNPFuncType param1_type = dpnp_dtype_to_DPNPFuncType(input.dtype)

cdef DPNPFuncData kernel_data = get_dpnp_function_ptr(DPNP_FN_MATRIX_RANK_EXT, param1_type, param1_type)

input_obj = input.get_array()

# create result array with type given by FPTR data
cdef utils.dpnp_descriptor result = utils.create_output_descriptor((1,),
kernel_data.return_type,
None,
device=input_obj.sycl_device,
usm_type=input_obj.usm_type,
sycl_queue=input_obj.sycl_queue)

result_sycl_queue = result.get_array().sycl_queue

cdef c_dpctl.SyclQueue q = <c_dpctl.SyclQueue> result_sycl_queue
cdef c_dpctl.DPCTLSyclQueueRef q_ref = q.get_queue_ref()

cdef custom_linalg_1in_1out_func_ptr_t func = <custom_linalg_1in_1out_func_ptr_t > kernel_data.ptr

cdef c_dpctl.DPCTLSyclEventRef event_ref = func(q_ref,
input.get_data(),
result.get_data(),
input_shape.data(),
input.ndim,
NULL) # dep_events_ref

with nogil: c_dpctl.DPCTLEvent_WaitAndThrow(event_ref)
c_dpctl.DPCTLEvent_Delete(event_ref)

return result


cpdef object dpnp_norm(object input, ord=None, axis=None):
cdef long size_input = input.size
cdef shape_type_c shape_input = input.shape
Expand Down
49 changes: 30 additions & 19 deletions dpnp/linalg/dpnp_iface_linalg.py
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,7 @@
dpnp_det,
dpnp_eigh,
dpnp_inv,
dpnp_matrix_rank,
dpnp_pinv,
dpnp_qr,
dpnp_slogdet,
Expand Down Expand Up @@ -397,47 +398,57 @@ def matrix_power(input, count):
return call_origin(numpy.linalg.matrix_power, input, count)


def matrix_rank(input, tol=None, hermitian=False):
def matrix_rank(A, tol=None, hermitian=False):
"""
Return matrix rank of array.
Return matrix rank of array using SVD method.
Rank of the array is the number of singular values of the array that are
greater than `tol`.
Parameters
----------
M : {(M,), (..., M, N)} array_like
A : {(M,), (..., M, N)} {dpnp.ndarray, usm_ndarray}
Input vector or stack of matrices.
tol : (...) array_like, float, optional
tol : (...) {float, dpnp.ndarray, usm_ndarray}, optional
Threshold below which SVD values are considered zero. If `tol` is
None, and ``S`` is an array with singular values for `M`, and
``eps`` is the epsilon value for datatype of ``S``, then `tol` is
set to ``S.max() * max(M.shape) * eps``.
hermitian : bool, optional
If True, `M` is assumed to be Hermitian (symmetric if real-valued),
If True, `A` is assumed to be Hermitian (symmetric if real-valued),
enabling a more efficient method for finding singular values.
Defaults to False.
Returns
-------
rank : (...) array_like
Rank of M.
rank : (...) dpnp.ndarray
Rank of A.
"""
See Also
--------
:obj:`dpnp.linalg.svd` : Singular Value Decomposition.
x1_desc = dpnp.get_dpnp_descriptor(input, copy_when_nondefault_queue=False)
if x1_desc:
if tol is not None:
pass
elif hermitian:
pass
else:
result_obj = dpnp_matrix_rank(x1_desc).get_pyobj()
result = dpnp.convert_single_elem_array_to_scalar(result_obj)
Examples
--------
>>> import dpnp as np
>>> from dpnp.linalg import matrix_rank
>>> matrix_rank(np.eye(4)) # Full rank matrix
array(4)
>>> I=np.eye(4); I[-1,-1] = 0. # rank deficient matrix
>>> matrix_rank(I)
array(3)
>>> matrix_rank(np.ones((4,))) # 1 dimension - rank 1 unless all 0
array(1)
>>> matrix_rank(np.zeros((4,)))
array(0)
return result
"""

dpnp.check_supported_arrays_type(A)
if tol is not None:
dpnp.check_supported_arrays_type(tol, scalar_type=True)

return call_origin(numpy.linalg.matrix_rank, input, tol, hermitian)
return dpnp_matrix_rank(A, tol=tol, hermitian=hermitian)


def multi_dot(arrays, out=None):
Expand Down
24 changes: 24 additions & 0 deletions dpnp/linalg/dpnp_utils_linalg.py
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,7 @@
"dpnp_det",
"dpnp_eigh",
"dpnp_inv",
"dpnp_matrix_rank",
"dpnp_pinv",
"dpnp_qr",
"dpnp_slogdet",
Expand Down Expand Up @@ -999,6 +1000,29 @@ def dpnp_inv(a):
return b_f


def dpnp_matrix_rank(A, tol=None, hermitian=False):
"""
dpnp_matrix_rank(A, tol=None, hermitian=False)
Return matrix rank of array using SVD method.
"""

if A.ndim < 2:
return (A != 0).any().astype(int)

S = dpnp_svd(A, compute_uv=False, hermitian=hermitian)

if tol is None:
rtol = max(A.shape[-2:]) * dpnp.finfo(S.dtype).eps
tol = S.max(axis=-1, keepdims=True) * rtol
elif not dpnp.isscalar(tol):
# Add a new axis to match Numpy's output
tol = tol[..., None]

return dpnp.count_nonzero(S > tol, axis=-1)


def dpnp_pinv(a, rcond=1e-15, hermitian=False):
"""
dpnp_pinv(a, rcond=1e-15, hermitian=False):
Expand Down
Loading

0 comments on commit d34b8ce

Please sign in to comment.