diff --git a/dpnp/backend/include/dpnp_iface_fft.hpp b/dpnp/backend/include/dpnp_iface_fft.hpp index 5663f14ace1c..47d9a7b952b0 100644 --- a/dpnp/backend/include/dpnp_iface_fft.hpp +++ b/dpnp/backend/include/dpnp_iface_fft.hpp @@ -53,23 +53,23 @@ * Compute the one-dimensional discrete Fourier Transform. * * @param[in] q_ref Reference to SYCL queue. - * @param[in] array_in Input array. - * @param[out] result Output array. + * @param[in] array1_in Input array. + * @param[out] result_out Output array. * @param[in] input_shape Array with shape information for input array. - * @param[in] output_shape Array with shape information for output array. - * @param[in] shape_size Number of elements in @ref input_shape or @ref output_shape arrays. + * @param[in] result_shape Array with shape information for result array. + * @param[in] shape_size Number of elements in @ref input_shape or @ref result_shape arrays. * @param[in] axis Axis ID to compute by. * @param[in] input_boundarie Limit number of elements for @ref axis. * @param[in] inverse Using inverse algorithm. * @param[in] norm Normalization mode. 0 - backward, 1 - forward, 2 - ortho. * @param[in] dep_event_vec_ref Reference to vector of SYCL events. */ -template +template INP_DLLEXPORT DPCTLSyclEventRef dpnp_fft_fft_c(DPCTLSyclQueueRef q_ref, - const void* array_in, - void* result, + const void* array1_in, + void* result_out, const shape_elem_type* input_shape, - const shape_elem_type* output_shape, + const shape_elem_type* result_shape, size_t shape_size, long axis, long input_boundarie, @@ -77,9 +77,52 @@ INP_DLLEXPORT DPCTLSyclEventRef dpnp_fft_fft_c(DPCTLSyclQueueRef q_ref, const size_t norm, const DPCTLEventVectorRef dep_event_vec_ref); -template -INP_DLLEXPORT void dpnp_fft_fft_c(const void* array_in, - void* result, +template +INP_DLLEXPORT void dpnp_fft_fft_c(const void* array1_in, + void* result_out, + const shape_elem_type* input_shape, + const shape_elem_type* output_shape, + size_t shape_size, + long axis, + long input_boundarie, + size_t inverse, + const size_t norm); + + +/** + * @ingroup BACKEND_FFT_API + * @brief 1D discrete Fourier Transform. + * + * Compute the one-dimensional discrete Fourier Transform for real input. + * + * @param[in] q_ref Reference to SYCL queue. + * @param[in] array1_in Input array. + * @param[out] result_out Output array. + * @param[in] input_shape Array with shape information for input array. + * @param[in] result_shape Array with shape information for result array. + * @param[in] shape_size Number of elements in @ref input_shape or @ref result_shape arrays. + * @param[in] axis Axis ID to compute by. + * @param[in] input_boundarie Limit number of elements for @ref axis. + * @param[in] inverse Using inverse algorithm. + * @param[in] norm Normalization mode. 0 - backward, 1 - forward, 2 - ortho. + * @param[in] dep_event_vec_ref Reference to vector of SYCL events. + */ +template +INP_DLLEXPORT DPCTLSyclEventRef dpnp_fft_rfft_c(DPCTLSyclQueueRef q_ref, + const void* array1_in, + void* result_out, + const shape_elem_type* input_shape, + const shape_elem_type* result_shape, + size_t shape_size, + long axis, + long input_boundarie, + size_t inverse, + const size_t norm, + const DPCTLEventVectorRef dep_event_vec_ref); + +template +INP_DLLEXPORT void dpnp_fft_fft_c(const void* array1_in, + void* result_out, const shape_elem_type* input_shape, const shape_elem_type* output_shape, size_t shape_size, diff --git a/dpnp/backend/include/dpnp_iface_fptr.hpp b/dpnp/backend/include/dpnp_iface_fptr.hpp index e79e5d2a1c24..63e2fd750a98 100644 --- a/dpnp/backend/include/dpnp_iface_fptr.hpp +++ b/dpnp/backend/include/dpnp_iface_fptr.hpp @@ -166,6 +166,8 @@ enum class DPNPFuncName : size_t DPNP_FN_FABS_EXT, /**< Used in numpy.fabs() impl, requires extra parameters */ DPNP_FN_FFT_FFT, /**< Used in numpy.fft.fft() impl */ DPNP_FN_FFT_FFT_EXT, /**< Used in numpy.fft.fft() impl, requires extra parameters */ + DPNP_FN_FFT_RFFT, /**< Used in numpy.fft.rfft() impl */ + DPNP_FN_FFT_RFFT_EXT, /**< Used in numpy.fft.rfft() impl, requires extra parameters */ DPNP_FN_FILL_DIAGONAL, /**< Used in numpy.fill_diagonal() impl */ DPNP_FN_FILL_DIAGONAL_EXT, /**< Used in numpy.fill_diagonal() impl, requires extra parameters */ DPNP_FN_FLATTEN, /**< Used in numpy.flatten() impl */ diff --git a/dpnp/backend/kernels/dpnp_krnl_fft.cpp b/dpnp/backend/kernels/dpnp_krnl_fft.cpp index 8c9aeb6e0be6..8a79a35e844b 100644 --- a/dpnp/backend/kernels/dpnp_krnl_fft.cpp +++ b/dpnp/backend/kernels/dpnp_krnl_fft.cpp @@ -50,7 +50,7 @@ class dpnp_fft_fft_c_kernel; template void dpnp_fft_fft_sycl_c(DPCTLSyclQueueRef q_ref, const void* array1_in, - void* result1, + void* result_out, const shape_elem_type* input_shape, const shape_elem_type* output_shape, size_t shape_size, @@ -73,7 +73,7 @@ void dpnp_fft_fft_sycl_c(DPCTLSyclQueueRef q_ref, DPNPC_ptr_adapter<_DataType_input> input1_ptr(q_ref, array1_in, input_size); const _DataType_input* array_1 = input1_ptr.get_ptr(); - _DataType_output* result = reinterpret_cast<_DataType_output*>(result1); + _DataType_output* result = reinterpret_cast<_DataType_output*>(result_out); // kernel specific temporal data shape_elem_type* output_shape_offsets = @@ -171,38 +171,44 @@ void dpnp_fft_fft_sycl_c(DPCTLSyclQueueRef q_ref, } template -void dpnp_fft_fft_mathlib_compute_c(DPCTLSyclQueueRef q_ref, - const void* array1_in, - void* result1, - const shape_elem_type* input_shape, - const size_t shape_size, - const size_t result_size, - _Descriptor_type& desc, - const size_t norm) +void dpnp_fft_fft_mathlib_cmplx_to_cmplx_c(DPCTLSyclQueueRef q_ref, + const void* array1_in, + void* result_out, + const shape_elem_type* input_shape, + const shape_elem_type* result_shape, + const size_t shape_size, + const size_t input_size, + const size_t result_size, + _Descriptor_type& desc, + const size_t norm) { if (!shape_size) { return; } - + sycl::queue queue = *(reinterpret_cast(q_ref)); - - DPNPC_ptr_adapter<_DataType_input> input1_ptr(q_ref, array1_in, result_size); - DPNPC_ptr_adapter<_DataType_output> result_ptr(q_ref, result1, result_size); + + DPNPC_ptr_adapter<_DataType_input> input1_ptr(q_ref, array1_in, input_size); + DPNPC_ptr_adapter<_DataType_output> result_ptr(q_ref, result_out, result_size); _DataType_input* array_1 = input1_ptr.get_ptr(); _DataType_output* result = result_ptr.get_ptr(); - desc.set_value(mkl_dft::config_param::BACKWARD_SCALE, (1.0 / result_size)); - // enum value from math library C interface - // instead of mkl_dft::config_value::NOT_INPLACE - desc.set_value(mkl_dft::config_param::PLACEMENT, DFTI_NOT_INPLACE); - desc.commit(queue); - const size_t n_iter = std::accumulate(input_shape, input_shape + shape_size - 1, 1, std::multiplies()); const size_t shift = input_shape[shape_size - 1]; + double forward_scale = 1.0; + double backward_scale = 1.0 / shift; + + desc.set_value(mkl_dft::config_param::BACKWARD_SCALE, backward_scale); + desc.set_value(mkl_dft::config_param::FORWARD_SCALE, forward_scale); + // enum value from math library C interface + // instead of mkl_dft::config_value::NOT_INPLACE + desc.set_value(mkl_dft::config_param::PLACEMENT, DFTI_NOT_INPLACE); + desc.commit(queue); + std::vector fft_events; fft_events.reserve(n_iter); @@ -215,48 +221,94 @@ void dpnp_fft_fft_mathlib_compute_c(DPCTLSyclQueueRef q_ref, return; } -// norm: backward - 0, forward is 1 -template -void dpnp_fft_fft_mathlib_c(DPCTLSyclQueueRef q_ref, - const void* array1_in, - void* result1, - const shape_elem_type* input_shape, - const size_t shape_size, - const size_t result_size, - const size_t norm) +template +class dpnp_fft_fft_mathlib_real_to_cmplx_c_kernel; + +template +void dpnp_fft_fft_mathlib_real_to_cmplx_c(DPCTLSyclQueueRef q_ref, + const void* array1_in, + void* result_out, + const shape_elem_type* input_shape, + const shape_elem_type* result_shape, + const size_t shape_size, + const size_t input_size, + const size_t result_size, + _Descriptor_type& desc, + const size_t norm, + const size_t real) { - if (!shape_size || !result_size || !array1_in || !result1) + if (!shape_size) { return; } - //will be used with strides - //std::vector dimensions(input_shape, input_shape + shape_size); - if constexpr (std::is_same<_DataType_input, std::complex>::value && - std::is_same<_DataType_output, std::complex>::value) - { - desc_dp_cmplx_t desc(input_shape[shape_size - 1]); + DPNPC_ptr_adapter<_DataType_input> input1_ptr(q_ref, array1_in, input_size); + DPNPC_ptr_adapter<_DataType_output> result_ptr(q_ref, result_out, result_size * 2, true, true); + _DataType_input* array_1 = input1_ptr.get_ptr(); + _DataType_output* result = result_ptr.get_ptr(); + + sycl::queue queue = *(reinterpret_cast(q_ref)); + + const size_t n_iter = + std::accumulate(input_shape, input_shape + shape_size - 1, 1, std::multiplies()); + + const size_t input_shift = input_shape[shape_size - 1]; + const size_t result_shift = result_shape[shape_size - 1];; + + double forward_scale = 1.0; + double backward_scale = 1.0 / input_shift; - dpnp_fft_fft_mathlib_compute_c<_DataType_input, _DataType_output, desc_dp_cmplx_t>( - q_ref, array1_in, result1, input_shape, shape_size, result_size, desc, norm); + desc.set_value(mkl_dft::config_param::BACKWARD_SCALE, backward_scale); + desc.set_value(mkl_dft::config_param::FORWARD_SCALE, forward_scale); + + desc.commit(queue); + + std::vector fft_events; + fft_events.reserve(n_iter); + + for (size_t i = 0; i < n_iter; ++i) { + fft_events.push_back(mkl_dft::compute_forward(desc, array_1 + i * input_shift, result + i * result_shift * 2)); } - else if (std::is_same<_DataType_input, std::complex>::value && - std::is_same<_DataType_output, std::complex>::value) - { - desc_sp_cmplx_t desc(input_shape[shape_size - 1]); - dpnp_fft_fft_mathlib_compute_c<_DataType_input, _DataType_output, desc_sp_cmplx_t>( - q_ref, array1_in, result1, input_shape, shape_size, result_size, desc, norm); + sycl::event::wait(fft_events); + + if (real) { // the output size of the rfft function is input_size/2 + 1 so we don't need to fill the second half of the output + return; } + + size_t n_conj = result_shift % 2 == 0 ? result_shift / 2 - 1 : result_shift / 2; + + sycl::event event; + + sycl::range<2> gws(n_iter, n_conj); + + auto kernel_parallel_for_func = [=](sycl::id<2> global_id) { + size_t i = global_id[0]; + { + size_t j = global_id[1]; + { + *(reinterpret_cast*>(result) + result_shift * (i + 1) - (j + 1)) = std::conj(*(reinterpret_cast*>(result) + result_shift * i + (j + 1))); + } + } + }; + + auto kernel_func = [&](sycl::handler& cgh) { + cgh.parallel_for>( + gws, kernel_parallel_for_func); + }; + + event = queue.submit(kernel_func); + event.wait(); + return; } template DPCTLSyclEventRef dpnp_fft_fft_c(DPCTLSyclQueueRef q_ref, const void* array1_in, - void* result1, + void* result_out, const shape_elem_type* input_shape, - const shape_elem_type* output_shape, + const shape_elem_type* result_shape, size_t shape_size, long axis, long input_boundarie, @@ -264,49 +316,90 @@ DPCTLSyclEventRef dpnp_fft_fft_c(DPCTLSyclQueueRef q_ref, const size_t norm, const DPCTLEventVectorRef dep_event_vec_ref) { - // avoid warning unused variable - (void)dep_event_vec_ref; - DPCTLSyclEventRef event_ref = nullptr; - if (!shape_size) + if (!shape_size || !array1_in || !result_out) { return event_ref; } const size_t result_size = - std::accumulate(output_shape, output_shape + shape_size, 1, std::multiplies()); + std::accumulate(result_shape, result_shape + shape_size, 1, std::multiplies()); const size_t input_size = std::accumulate(input_shape, input_shape + shape_size, 1, std::multiplies()); - if (!input_size || !result_size || !array1_in || !result1) - { - return event_ref; - } - - sycl::queue q = *(reinterpret_cast(q_ref)); + size_t dim = input_shape[shape_size - 1]; - if ((std::is_same<_DataType_input, std::complex>::value && - std::is_same<_DataType_output, std::complex>::value) || - (std::is_same<_DataType_input, std::complex>::value && - std::is_same<_DataType_output, std::complex>::value)) - { - dpnp_fft_fft_mathlib_c<_DataType_input, _DataType_output>( - q_ref, array1_in, result1, input_shape, shape_size, result_size, norm); - } - else + if constexpr (std::is_same<_DataType_output, std::complex>::value || + std::is_same<_DataType_output, std::complex>::value) { - dpnp_fft_fft_sycl_c<_DataType_input, _DataType_output>(q_ref, - array1_in, - result1, - input_shape, - output_shape, - shape_size, - result_size, - input_size, - axis, - input_boundarie, - inverse); + if constexpr (std::is_same<_DataType_input, std::complex>::value && + std::is_same<_DataType_output, std::complex>::value) + { + desc_dp_cmplx_t desc(dim); + dpnp_fft_fft_mathlib_cmplx_to_cmplx_c<_DataType_input, _DataType_output, desc_dp_cmplx_t>( + q_ref, array1_in, result_out, input_shape, result_shape, shape_size, input_size, result_size, desc, norm); + } + /* complex-to-complex, single precision */ + else if constexpr (std::is_same<_DataType_input, std::complex>::value && + std::is_same<_DataType_output, std::complex>::value) + { + desc_sp_cmplx_t desc(dim); + dpnp_fft_fft_mathlib_cmplx_to_cmplx_c<_DataType_input, _DataType_output, desc_sp_cmplx_t>( + q_ref, array1_in, result_out, input_shape, result_shape, shape_size, input_size, result_size, desc, norm); + } + /* real-to-complex, double precision */ + else if constexpr (std::is_same<_DataType_input, double>::value && + std::is_same<_DataType_output, std::complex>::value) + { + desc_dp_real_t desc(dim); + + dpnp_fft_fft_mathlib_real_to_cmplx_c<_DataType_input, double, desc_dp_real_t>( + q_ref, array1_in, result_out, input_shape, result_shape, shape_size, input_size, result_size, desc, norm, 0); + } + /* real-to-complex, single precision */ + else if constexpr (std::is_same<_DataType_input, float>::value && + std::is_same<_DataType_output, std::complex>::value) + { + desc_sp_real_t desc(dim); // try: 2 * result_size + dpnp_fft_fft_mathlib_real_to_cmplx_c<_DataType_input, float, desc_sp_real_t>( + q_ref, array1_in, result_out, input_shape, result_shape, shape_size, input_size, result_size, desc, norm, 0); + } + else if constexpr (std::is_same<_DataType_input, int32_t>::value || + std::is_same<_DataType_input, int64_t>::value) + { + double* array1_copy = reinterpret_cast(dpnp_memory_alloc_c(input_size * sizeof(double))); + + shape_elem_type* copy_strides = reinterpret_cast(dpnp_memory_alloc_c(q_ref, sizeof(shape_elem_type))); + *copy_strides = 1; + shape_elem_type* copy_shape = reinterpret_cast(dpnp_memory_alloc_c(q_ref, sizeof(shape_elem_type))); + *copy_shape = input_size; + shape_elem_type copy_shape_size = 1; + dpnp_copyto_c<_DataType_input, double>(q_ref, array1_copy, input_size, copy_shape_size, copy_shape, copy_strides, + array1_in, input_size, copy_shape_size, copy_shape, copy_strides, NULL, dep_event_vec_ref); + + desc_dp_real_t desc(dim); + dpnp_fft_fft_mathlib_real_to_cmplx_c( + q_ref, array1_copy, result_out, input_shape, result_shape, shape_size, input_size, result_size, desc, norm, 0); + + dpnp_memory_free_c(q_ref, array1_copy); + dpnp_memory_free_c(q_ref, copy_strides); + dpnp_memory_free_c(q_ref, copy_shape); + } + else + { + dpnp_fft_fft_sycl_c<_DataType_input, _DataType_output>(q_ref, + array1_in, + result_out, + input_shape, + result_shape, + shape_size, + result_size, + input_size, + axis, + input_boundarie, + inverse); + } } return event_ref; @@ -363,6 +456,131 @@ DPCTLSyclEventRef (*dpnp_fft_fft_ext_c)(DPCTLSyclQueueRef, const size_t, const DPCTLEventVectorRef) = dpnp_fft_fft_c<_DataType_input, _DataType_output>; + +template +DPCTLSyclEventRef dpnp_fft_rfft_c(DPCTLSyclQueueRef q_ref, + const void* array1_in, + void* result_out, + const shape_elem_type* input_shape, + const shape_elem_type* result_shape, + size_t shape_size, + long axis, + long input_boundarie, + size_t inverse, + const size_t norm, + const DPCTLEventVectorRef dep_event_vec_ref) +{ + DPCTLSyclEventRef event_ref = nullptr; + + if (!shape_size || !array1_in || !result_out) + { + return event_ref; + } + + const size_t result_size = + std::accumulate(result_shape, result_shape + shape_size, 1, std::multiplies()); + const size_t input_size = + std::accumulate(input_shape, input_shape + shape_size, 1, std::multiplies()); + + size_t dim = input_shape[shape_size - 1]; + + if constexpr (std::is_same<_DataType_output, std::complex>::value || + std::is_same<_DataType_output, std::complex>::value) + { + if constexpr (std::is_same<_DataType_input, double>::value && + std::is_same<_DataType_output, std::complex>::value) + { + desc_dp_real_t desc(dim); + + dpnp_fft_fft_mathlib_real_to_cmplx_c<_DataType_input, double, desc_dp_real_t>( + q_ref, array1_in, result_out, input_shape, result_shape, shape_size, input_size, result_size, desc, norm, 1l); + } + /* real-to-complex, single precision */ + else if constexpr (std::is_same<_DataType_input, float>::value && + std::is_same<_DataType_output, std::complex>::value) + { + desc_sp_real_t desc(dim); // try: 2 * result_size + dpnp_fft_fft_mathlib_real_to_cmplx_c<_DataType_input, float, desc_sp_real_t>( + q_ref, array1_in, result_out, input_shape, result_shape, shape_size, input_size, result_size, desc, norm, 1); + } + else if constexpr (std::is_same<_DataType_input, int32_t>::value || + std::is_same<_DataType_input, int64_t>::value) + { + double* array1_copy = reinterpret_cast(dpnp_memory_alloc_c(input_size * sizeof(double))); + + shape_elem_type* copy_strides = reinterpret_cast(dpnp_memory_alloc_c(q_ref, sizeof(shape_elem_type))); + *copy_strides = 1; + shape_elem_type* copy_shape = reinterpret_cast(dpnp_memory_alloc_c(q_ref, sizeof(shape_elem_type))); + *copy_shape = input_size; + shape_elem_type copy_shape_size = 1; + dpnp_copyto_c<_DataType_input, double>(q_ref, array1_copy, input_size, copy_shape_size, copy_shape, copy_strides, + array1_in, input_size, copy_shape_size, copy_shape, copy_strides, NULL, dep_event_vec_ref); + + desc_dp_real_t desc(dim); + dpnp_fft_fft_mathlib_real_to_cmplx_c( + q_ref, array1_copy, result_out, input_shape, result_shape, shape_size, input_size, result_size, desc, norm, 1); + + dpnp_memory_free_c(q_ref, array1_copy); + dpnp_memory_free_c(q_ref, copy_strides); + dpnp_memory_free_c(q_ref, copy_shape); + } + } + + return event_ref; +} + + +template +void dpnp_fft_rfft_c(const void* array1_in, + void* result1, + const shape_elem_type* input_shape, + const shape_elem_type* output_shape, + size_t shape_size, + long axis, + long input_boundarie, + size_t inverse, + const size_t norm) +{ + DPCTLSyclQueueRef q_ref = reinterpret_cast(&DPNP_QUEUE); + DPCTLEventVectorRef dep_event_vec_ref = nullptr; + DPCTLSyclEventRef event_ref = dpnp_fft_rfft_c<_DataType_input, _DataType_output>(q_ref, + array1_in, + result1, + input_shape, + output_shape, + shape_size, + axis, + input_boundarie, + inverse, + norm, + dep_event_vec_ref); + DPCTLEvent_WaitAndThrow(event_ref); +} + +template +void (*dpnp_fft_rfft_default_c)(const void*, + void*, + const shape_elem_type*, + const shape_elem_type*, + size_t, + long, + long, + size_t, + const size_t) = dpnp_fft_rfft_c<_DataType_input, _DataType_output>; + +template +DPCTLSyclEventRef (*dpnp_fft_rfft_ext_c)(DPCTLSyclQueueRef, + const void*, + void*, + const shape_elem_type*, + const shape_elem_type*, + size_t, + long, + long, + size_t, + const size_t, + const DPCTLEventVectorRef) = dpnp_fft_rfft_c<_DataType_input, _DataType_output>; + void func_map_init_fft_func(func_map_t& fmap) { fmap[DPNPFuncName::DPNP_FN_FFT_FFT][eft_INT][eft_INT] = { @@ -377,5 +595,13 @@ void func_map_init_fft_func(func_map_t& fmap) eft_C64, (void*)dpnp_fft_fft_default_c, std::complex>}; fmap[DPNPFuncName::DPNP_FN_FFT_FFT][eft_C128][eft_C128] = { eft_C128, (void*)dpnp_fft_fft_default_c, std::complex>}; + fmap[DPNPFuncName::DPNP_FN_FFT_RFFT][eft_INT][eft_INT] = { + eft_C128, (void*)dpnp_fft_rfft_default_c>}; + fmap[DPNPFuncName::DPNP_FN_FFT_RFFT][eft_LNG][eft_LNG] = { + eft_C128, (void*)dpnp_fft_rfft_default_c>}; + fmap[DPNPFuncName::DPNP_FN_FFT_RFFT][eft_FLT][eft_FLT] = { + eft_C64, (void*)dpnp_fft_rfft_default_c>}; + fmap[DPNPFuncName::DPNP_FN_FFT_RFFT][eft_DBL][eft_DBL] = { + eft_C128, (void*)dpnp_fft_rfft_default_c>}; return; } diff --git a/dpnp/dpnp_algo/dpnp_algo.pxd b/dpnp/dpnp_algo/dpnp_algo.pxd index 1cc650bdaaa8..45453e9f47a1 100644 --- a/dpnp/dpnp_algo/dpnp_algo.pxd +++ b/dpnp/dpnp_algo/dpnp_algo.pxd @@ -87,6 +87,7 @@ cdef extern from "dpnp_iface_fptr.hpp" namespace "DPNPFuncName": # need this na DPNP_FN_EXPM1 DPNP_FN_FABS DPNP_FN_FFT_FFT + DPNP_FN_FFT_RFFT DPNP_FN_FILL_DIAGONAL DPNP_FN_FLATTEN DPNP_FN_FLOOR diff --git a/dpnp/fft/dpnp_algo_fft.pyx b/dpnp/fft/dpnp_algo_fft.pyx index 92bf7fee83b4..d63c7bf9fc68 100644 --- a/dpnp/fft/dpnp_algo_fft.pyx +++ b/dpnp/fft/dpnp_algo_fft.pyx @@ -38,7 +38,8 @@ cimport dpnp.dpnp_utils as utils __all__ = [ - "dpnp_fft" + "dpnp_fft", + "dpnp_rfft" ] ctypedef void(*fptr_dpnp_fft_fft_t)(void *, void * , shape_elem_type * , shape_elem_type * , @@ -73,3 +74,33 @@ cpdef utils.dpnp_descriptor dpnp_fft(utils.dpnp_descriptor input, output_shape.data(), input_shape.size(), axis_norm, input_boundarie, inverse, norm) return result + + +cpdef utils.dpnp_descriptor dpnp_rfft(utils.dpnp_descriptor input, + size_t input_boundarie, + size_t output_boundarie, + long axis, + size_t inverse, + size_t norm): + + cdef shape_type_c input_shape = input.shape + cdef shape_type_c output_shape = input_shape + + cdef long axis_norm = utils.normalize_axis((axis,), input_shape.size())[0] + output_shape[axis_norm] = output_boundarie + + # convert string type names (dtype) to C enum DPNPFuncType + cdef DPNPFuncType param1_type = dpnp_dtype_to_DPNPFuncType(input.dtype) + + # get the FPTR data structure + cdef DPNPFuncData kernel_data = get_dpnp_function_ptr(DPNP_FN_FFT_RFFT, param1_type, param1_type) + + # ceate result array with type given by FPTR data + cdef utils.dpnp_descriptor result = utils.create_output_descriptor(output_shape, kernel_data.return_type, None) + + cdef fptr_dpnp_fft_fft_t func = kernel_data.ptr + # call FPTR function + func(input.get_data(), result.get_data(), input_shape.data(), + output_shape.data(), input_shape.size(), axis_norm, input_boundarie, inverse, norm) + + return result diff --git a/dpnp/fft/dpnp_iface_fft.py b/dpnp/fft/dpnp_iface_fft.py index 61034efc4398..d245de85054a 100644 --- a/dpnp/fft/dpnp_iface_fft.py +++ b/dpnp/fft/dpnp_iface_fft.py @@ -120,9 +120,9 @@ def fft(x1, n=None, axis=-1, norm=None): pass # let fallback to handle exception elif norm is not None: pass - elif axis != -1: + elif n is not None: pass - elif x1_desc.dtype not in (numpy.complex128, numpy.complex64): + elif axis != -1: pass else: output_boundarie = input_boundarie @@ -327,6 +327,8 @@ def ifft(x1, n=None, axis=-1, norm=None): pass # let fallback to handle exception elif norm is not None: pass + elif n is not None: + pass else: output_boundarie = input_boundarie @@ -476,6 +478,8 @@ def ihfft(x1, n=None, axis=-1, norm=None): pass # let fallback to handle exception elif norm is not None: pass + elif n is not None: + pass else: output_boundarie = input_boundarie @@ -518,10 +522,12 @@ def irfft(x1, n=None, axis=-1, norm=None): pass # let fallback to handle exception elif norm is not None: pass + elif n is not None: + pass else: output_boundarie = 2 * (input_boundarie - 1) - result = dpnp_fft(x1_desc, input_boundarie, output_boundarie, axis_param, True, norm_.value).get_pyobj() + result = dpnp_rfft(x1_desc, input_boundarie, output_boundarie, axis_param, True, norm_.value).get_pyobj() # TODO tmp = utils.create_output_array(result_shape, result_c_type, out) # tmp = dparray(result.shape, dtype=dpnp.float64) # for it in range(tmp.size): @@ -638,16 +644,17 @@ def rfft(x1, n=None, axis=-1, norm=None): pass # let fallback to handle exception elif input_boundarie < 1: pass # let fallback to handle exception + elif axis != -1: + pass elif norm is not None: pass - elif x1_desc.ndim > 1: + elif n is not None: pass - elif x1_desc.dtype not in (numpy.complex128, numpy.complex64): + elif x1_desc.dtype in (numpy.complex128, numpy.complex64): pass else: output_boundarie = input_boundarie // 2 + 1 # rfft specific requirenment - - return dpnp_fft(x1_desc, input_boundarie, output_boundarie, axis_param, False, norm_.value).get_pyobj() + return dpnp_rfft(x1_desc, input_boundarie, output_boundarie, axis_param, False, norm_.value).get_pyobj() return call_origin(numpy.fft.rfft, x1, n, axis, norm) diff --git a/tests/test_fft.py b/tests/test_fft.py index 6e7b43f0d569..703495f3b585 100644 --- a/tests/test_fft.py +++ b/tests/test_fft.py @@ -32,3 +32,16 @@ def test_fft_ndim(type, shape): numpy.testing.assert_allclose(dpnp_res, np_res, rtol=1e-4, atol=1e-7) assert dpnp_res.dtype == np_res.dtype + + +@pytest.mark.parametrize("type", ['float32', 'float64', 'int32', 'int64']) +@pytest.mark.parametrize("shape", [(64, ), (8, 8), (4, 16), (4, 4, 4), (2, 4, 4, 2)]) +def test_fft_rfft(type, shape): + np_data = numpy.arange(64, dtype=numpy.dtype(type)).reshape(shape) + dpnp_data = dpnp.arange(64, dtype=numpy.dtype(type)).reshape(shape) + + np_res = numpy.fft.rfft(np_data) + dpnp_res = dpnp.fft.rfft(dpnp_data) + + numpy.testing.assert_allclose(dpnp_res, np_res, rtol=1e-4, atol=1e-7) + assert dpnp_res.dtype == np_res.dtype