Linalg svd func (#492)

* linalg.svd() impl

Author: Alexander-Makaryev

dpnp/backend/include/dpnp_iface_fptr.hpp

@@ -153,6 +153,7 @@ enum class DPNPFuncName : size_t
     DPNP_FN_STD,                      /**< Used in numpy.std() implementation  */
     DPNP_FN_SUBTRACT,                 /**< Used in numpy.subtract() implementation  */
     DPNP_FN_SUM,                      /**< Used in numpy.sum() implementation  */
+    DPNP_FN_SVD,                      /**< Used in numpy.linalg.svd() implementation  */
     DPNP_FN_TAN,                      /**< Used in numpy.tan() implementation  */
     DPNP_FN_TANH,                     /**< Used in numpy.tanh() implementation  */
     DPNP_FN_TRANSPOSE,                /**< Used in numpy.transpose() implementation  */

dpnp/backend/kernels/dpnp_krnl_linalg.cpp

@@ -32,6 +32,7 @@
 #include "queue_sycl.hpp"
 
 namespace mkl_blas = oneapi::mkl::blas::row_major;
+namespace mkl_lapack = oneapi::mkl::lapack;
 
 template <typename _DataType>
 class dpnp_cholesky_c_kernel;
@@ -309,6 +310,58 @@ void dpnp_matrix_rank_c(void* array1_in, void* result1, size_t* shape, size_t nd
     return;
 }
 
+template <typename _InputDT, typename _ComputeDT, typename _SVDT>
+void dpnp_svd_c(void* array1_in, void* result1, void* result2, void* result3, size_t size_m, size_t size_n)
+{
+    cl::sycl::event event;
+
+    _InputDT* in_array = reinterpret_cast<_InputDT*>(array1_in);
+
+    // math lib gesvd func overrides input
+    _ComputeDT* in_a = reinterpret_cast<_ComputeDT*>(dpnp_memory_alloc_c(size_m * size_n * sizeof(_ComputeDT)));
+    for (size_t it = 0; it < size_m * size_n; ++it)
+    {
+        in_a[it] = in_array[it];
+    }
+
+    _ComputeDT* res_u = reinterpret_cast<_ComputeDT*>(result1);
+    _SVDT* res_s = reinterpret_cast<_SVDT*>(result2);
+    _ComputeDT* res_vt = reinterpret_cast<_ComputeDT*>(result3);
+
+    const std::int64_t m = size_m;
+    const std::int64_t n = size_n;
+
+    const std::int64_t lda = std::max<size_t>(1UL, n);
+    const std::int64_t ldu = std::max<size_t>(1UL, m);
+    const std::int64_t ldvt = std::max<size_t>(1UL, n);
+
+    const std::int64_t scratchpad_size1 = mkl_lapack::gesvd_scratchpad_size<_ComputeDT>(
+        DPNP_QUEUE, oneapi::mkl::jobsvd::vectors, oneapi::mkl::jobsvd::vectors, n, m, lda, ldvt, ldu);
+
+    const std::int64_t scratchpad_size = scratchpad_size1;
+
+    _ComputeDT* scratchpad = reinterpret_cast<_ComputeDT*>(dpnp_memory_alloc_c(scratchpad_size * sizeof(_ComputeDT)));
+
+    event = mkl_lapack::gesvd(DPNP_QUEUE,
+                              oneapi::mkl::jobsvd::vectors, // onemkl::job jobu,
+                              oneapi::mkl::jobsvd::vectors, // onemkl::job jobvt,
+                              n,
+                              m,
+                              in_a,
+                              lda,
+                              res_s,
+                              res_vt,
+                              ldvt,
+                              res_u,
+                              ldu,
+                              scratchpad,
+                              scratchpad_size);
+
+    event.wait();
+
+    dpnp_memory_free_c(scratchpad);
+}
+
 void func_map_init_linalg_func(func_map_t& fmap)
 {
     fmap[DPNPFuncName::DPNP_FN_CHOLESKY][eft_INT][eft_INT] = {eft_INT, (void*)dpnp_cholesky_c<int>};
@@ -331,5 +384,12 @@ void func_map_init_linalg_func(func_map_t& fmap)
     fmap[DPNPFuncName::DPNP_FN_MATRIX_RANK][eft_FLT][eft_FLT] = {eft_FLT, (void*)dpnp_matrix_rank_c<float>};
     fmap[DPNPFuncName::DPNP_FN_MATRIX_RANK][eft_DBL][eft_DBL] = {eft_DBL, (void*)dpnp_matrix_rank_c<double>};
 
+    fmap[DPNPFuncName::DPNP_FN_SVD][eft_INT][eft_INT] = {eft_DBL, (void*)dpnp_svd_c<int, double, double>};
+    fmap[DPNPFuncName::DPNP_FN_SVD][eft_LNG][eft_LNG] = {eft_DBL, (void*)dpnp_svd_c<long, double, double>};
+    fmap[DPNPFuncName::DPNP_FN_SVD][eft_FLT][eft_FLT] = {eft_FLT, (void*)dpnp_svd_c<float, float, float>};
+    fmap[DPNPFuncName::DPNP_FN_SVD][eft_DBL][eft_DBL] = {eft_DBL, (void*)dpnp_svd_c<double, double, double>};
+    fmap[DPNPFuncName::DPNP_FN_SVD][eft_C128][eft_C128] = {
+        eft_C128, (void*)dpnp_svd_c<std::complex<double>, std::complex<double>, double>};
+
     return;
 }

dpnp/dpnp_algo/dpnp_algo.pxd

@@ -126,6 +126,7 @@ cdef extern from "dpnp_iface_fptr.hpp" namespace "DPNPFuncName":  # need this na
         DPNP_FN_STD
         DPNP_FN_SUBTRACT
         DPNP_FN_SUM
+        DPNP_FN_SVD
         DPNP_FN_TAN
         DPNP_FN_TANH
         DPNP_FN_TRANSPOSE

dpnp/linalg/dpnp_algo_linalg.pyx

@@ -48,20 +48,16 @@ __all__ = [
     "dpnp_eigvals",
     "dpnp_inv",
     "dpnp_matrix_rank",
-    "dpnp_norm"
+    "dpnp_norm",
+    "dpnp_svd",
 ]
 
 
 # C function pointer to the C library template functions
 ctypedef void(*custom_linalg_1in_1out_func_ptr_t)(void *, void * , size_t * , size_t)
-
-
-# C function pointer to the C library template functions
 ctypedef void(*custom_linalg_1in_1out_func_ptr_t_)(void * , void * , size_t * )
-
-
-# C function pointer to the C library template functions
 ctypedef void(*custom_linalg_1in_1out_with_size_func_ptr_t_)(void *, void * , size_t)
+ctypedef void(*custom_linalg_1in_3out_shape_t)(void *, void * , void * , void * , size_t , size_t )
 
 
 cpdef dparray dpnp_cholesky(dparray input):
@@ -300,3 +296,30 @@ cpdef dparray dpnp_norm(dparray input, ord=None, axis=None):
         return ret
     else:
         raise ValueError("Improper number of dimensions to norm.")
+
+
+cpdef tuple dpnp_svd(dparray x1, full_matrices, compute_uv, hermitian):
+    cdef size_t size_m = x1.shape[0]
+    cdef size_t size_n = x1.shape[1]
+
+    cdef DPNPFuncType param1_type = dpnp_dtype_to_DPNPFuncType(x1.dtype)
+    cdef DPNPFuncData kernel_data = get_dpnp_function_ptr(DPNP_FN_SVD, param1_type, param1_type)
+
+    result_type = dpnp_DPNPFuncType_to_dtype(< size_t > kernel_data.return_type)
+
+    if x1.dtype == dpnp.float32:
+        type_s = dpnp.float32
+    else:
+        type_s = dpnp.float64
+
+    size_s = min(size_m, size_n)
+
+    cdef dparray res_u = dparray((size_m, size_m), dtype=result_type)
+    cdef dparray res_s = dparray((size_s, ), dtype=type_s)
+    cdef dparray res_vt = dparray((size_n, size_n), dtype=result_type)
+
+    cdef custom_linalg_1in_3out_shape_t func = < custom_linalg_1in_3out_shape_t > kernel_data.ptr
+
+    func(x1.get_data(), res_u.get_data(), res_s.get_data(), res_vt.get_data(), size_m, size_n)
+
+    return (res_u, res_s, res_vt)

dpnp/linalg/dpnp_iface_linalg.py

@@ -58,7 +58,8 @@
     "matrix_power",
     "matrix_rank",
     "multi_dot",
-    "norm"
+    "norm",
+    "svd",
 ]
 
 
@@ -394,3 +395,77 @@ def norm(input, ord=None, axis=None, keepdims=False):
         return result
 
     return call_origin(numpy.linalg.norm, input, ord, axis, keepdims)
+
+
+def svd(a, full_matrices=True, compute_uv=True, hermitian=False):
+    """
+    Singular Value Decomposition.
+
+    For full documentation refer to :obj:`numpy.linalg.svd`.
+
+    Examples
+    --------
+    >>> import dpnp as np
+    >>> a = np.random.randn(9, 6) + 1j*np.random.randn(9, 6)
+    >>> b = np.random.randn(2, 7, 8, 3) + 1j*np.random.randn(2, 7, 8, 3)
+
+    Reconstruction based on full SVD, 2D case:
+
+    >>> u, s, vh = np.linalg.svd(a, full_matrices=True)
+    >>> u.shape, s.shape, vh.shape
+    ((9, 9), (6,), (6, 6))
+    >>> np.allclose(a, np.dot(u[:, :6] * s, vh))
+    True
+    >>> smat = np.zeros((9, 6), dtype=complex)
+    >>> smat[:6, :6] = np.diag(s)
+    >>> np.allclose(a, np.dot(u, np.dot(smat, vh)))
+    True
+
+    Reconstruction based on reduced SVD, 2D case:
+
+    >>> u, s, vh = np.linalg.svd(a, full_matrices=False)
+    >>> u.shape, s.shape, vh.shape
+    ((9, 6), (6,), (6, 6))
+    >>> np.allclose(a, np.dot(u * s, vh))
+    True
+    >>> smat = np.diag(s)
+    >>> np.allclose(a, np.dot(u, np.dot(smat, vh)))
+    True
+
+    Reconstruction based on full SVD, 4D case:
+
+    >>> u, s, vh = np.linalg.svd(b, full_matrices=True)
+    >>> u.shape, s.shape, vh.shape
+    ((2, 7, 8, 8), (2, 7, 3), (2, 7, 3, 3))
+    >>> np.allclose(b, np.matmul(u[..., :3] * s[..., None, :], vh))
+    True
+    >>> np.allclose(b, np.matmul(u[..., :3], s[..., None] * vh))
+    True
+
+    Reconstruction based on reduced SVD, 4D case:
+
+    >>> u, s, vh = np.linalg.svd(b, full_matrices=False)
+    >>> u.shape, s.shape, vh.shape
+    ((2, 7, 8, 3), (2, 7, 3), (2, 7, 3, 3))
+    >>> np.allclose(b, np.matmul(u * s[..., None, :], vh))
+    True
+    >>> np.allclose(b, np.matmul(u, s[..., None] * vh))
+    True
+
+    """
+
+    if not use_origin_backend(a):
+        if not isinstance(a, dparray):
+            pass
+        if not a.ndim == 2:
+            pass
+        if not full_matrices == True:
+            pass
+        if not compute_uv == True:
+            pass
+        if not hermitian == False:
+            pass
+        else:
+            return dpnp_svd(a, full_matrices, compute_uv, hermitian)
+
+    return call_origin(numpy.linalg.svd, a, full_matrices, compute_uv, hermitian)

tests/skipped_tests.tbl

@@ -137,6 +137,8 @@ tests/test_linalg.py::test_norm3[(1, 2)-2-[[[1, 2], [3, 4]], [[5, 6], [7, 8]]]]
 tests/test_linalg.py::test_norm3[(1, 2)-2-[[[1, 0], [3, 0]], [[5, 0], [7, 0]]]]
 tests/test_linalg.py::test_norm3[(1, 2)-3-[[[1, 2], [3, 4]], [[5, 6], [7, 8]]]]
 tests/test_linalg.py::test_norm3[(1, 2)-3-[[[1, 0], [3, 0]], [[5, 0], [7, 0]]]]
+tests/test_linalg.py::test_svd[(3,4)-complex128]
+tests/test_linalg.py::test_svd[(5,3)-complex128]
 tests/test_random.py::TestDistributionsRayleigh::test_moments
 tests/third_party/cupy/binary_tests/test_elementwise.py::TestElementwise::test_bitwise_and
 tests/third_party/cupy/binary_tests/test_elementwise.py::TestElementwise::test_bitwise_or

tests/skipped_tests_gpu.tbl

@@ -142,6 +142,20 @@ tests/test_linalg.py::test_norm3[None--numpy.Inf-[[[1, 0], [3, 0]], [[5, 0], [7,
 tests/test_linalg.py::test_norm3[None-numpy.Inf-[[[1, 0], [3, 0]], [[5, 0], [7, 0]]]]
 tests/test_linalg.py::test_norm3[None--numpy.Inf-[[[1, 2], [3, 4]], [[5, 6], [7, 8]]]]
 tests/test_linalg.py::test_norm3[None-numpy.Inf-[[[1, 2], [3, 4]], [[5, 6], [7, 8]]]]
+tests/test_linalg.py::test_svd[(3,4)-float64]
+tests/test_linalg.py::test_svd[(3,4)-float32]
+tests/test_linalg.py::test_svd[(3,4)-int64]
+tests/test_linalg.py::test_svd[(3,4)-int32]
+tests/test_linalg.py::test_svd[(3,4)-complex128]
+tests/test_linalg.py::test_svd[(5,3)-float64]
+tests/test_linalg.py::test_svd[(5,3)-float32]
+tests/test_linalg.py::test_svd[(5,3)-int64]
+tests/test_linalg.py::test_svd[(5,3)-int32]
+tests/test_linalg.py::test_svd[(5,3)-complex128]
+tests/test_linalg.py::test_svd[(16,16)-float64]
+tests/test_linalg.py::test_svd[(16,16)-float32]
+tests/test_linalg.py::test_svd[(16,16)-int64]
+tests/test_linalg.py::test_svd[(16,16)-int32]
 tests/test_random.py::TestDistributionsRayleigh::test_moments
 tests/test_statistics.py::test_median[2-float32]
 tests/test_statistics.py::test_median[2-float64]

tests/test_linalg.py

@@ -197,3 +197,49 @@ def test_norm3(array, ord, axis):
     result = inp.linalg.norm(ia, ord=ord, axis=axis)
     expected = numpy.linalg.norm(a, ord=ord, axis=axis)
     numpy.testing.assert_array_equal(expected, result)
+
+
+@pytest.mark.parametrize("type",
+                         [numpy.float64, numpy.float32, numpy.int64, numpy.int32, numpy.complex128],
+                         ids=['float64', 'float32', 'int64', 'int32', 'complex128'])
+@pytest.mark.parametrize("shape",
+                         [(2,2), (3,4), (5,3), (16,16)],
+                         ids=['(2,2)', '(3,4)', '(5,3)', '(16,16)'])
+def test_svd(type, shape):
+    a = numpy.arange(shape[0] * shape[1], dtype=type).reshape(shape)
+    ia = inp.array(a)
+
+    np_u, np_s, np_vt = numpy.linalg.svd(a)
+    dpnp_u, dpnp_s, dpnp_vt = inp.linalg.svd(ia)
+
+    assert (dpnp_u.dtype == np_u.dtype)
+    assert (dpnp_s.dtype == np_s.dtype)
+    assert (dpnp_vt.dtype == np_vt.dtype)
+    assert (dpnp_u.shape == np_u.shape)
+    assert (dpnp_s.shape == np_s.shape)
+    assert (dpnp_vt.shape == np_vt.shape)
+
+    if type == numpy.float32:
+        tol = 1e-03
+    else:
+        tol = 1e-12
+
+    # check decomposition
+    dpnp_diag_s = numpy.zeros(shape, dtype=dpnp_s.dtype)
+    for i in range(len(dpnp_s)):
+        dpnp_diag_s[i, i] = dpnp_s[i]
+    numpy.testing.assert_allclose(ia, numpy.dot(dpnp_u, numpy.dot(dpnp_diag_s, dpnp_vt)), rtol=tol, atol=tol)
+
+    # compare singular values
+    numpy.testing.assert_allclose(dpnp_s, np_s, rtol=tol, atol=tol)
+
+    # change sign of vectors
+    for i in range(min(shape[0], shape[1])):
+        if np_u[0, i] * dpnp_u[0, i] < 0:
+            np_u[:, i] = -np_u[:, i]
+            np_vt[i, :] = -np_vt[i, :]
+
+    # compare vectors for non-zero values
+    for i in range(numpy.count_nonzero(np_s > tol)):
+        numpy.testing.assert_allclose(numpy.array(dpnp_u)[:, i], np_u[:, i], rtol=tol, atol=tol)
+        numpy.testing.assert_allclose(numpy.array(dpnp_vt)[i, :], np_vt[i, :], rtol=tol, atol=tol)