Update dpnp.linalg.matrix_rank() implementation (#1717)

* 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

Author: vlad-perevezentsev

dpnp/backend/include/dpnp_iface_fptr.hpp

@@ -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  */

dpnp/backend/kernels/dpnp_krnl_linalg.cpp

@@ -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,
@@ -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] = {

dpnp/dpnp_algo/dpnp_algo.pxd

@@ -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

dpnp/linalg/dpnp_algo_linalg.pyx

@@ -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)
@@ -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

dpnp/linalg/dpnp_iface_linalg.py

@@ -51,6 +51,7 @@
     dpnp_det,
     dpnp_eigh,
     dpnp_inv,
+    dpnp_matrix_rank,
     dpnp_pinv,
     dpnp_qr,
     dpnp_slogdet,
@@ -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):

dpnp/linalg/dpnp_utils_linalg.py

@@ -39,6 +39,7 @@
     "dpnp_det",
     "dpnp_eigh",
     "dpnp_inv",
+    "dpnp_matrix_rank",
     "dpnp_pinv",
     "dpnp_qr",
     "dpnp_slogdet",
@@ -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):