FFT. Inverse version for fft(), fft2() and fftn() (#350)

Author: shssf

dpnp/backend/backend_iface_fft.hpp

@@ -69,6 +69,7 @@
  * @param[in]  shape_size      Number of elements in @ref input_shape or @ref output_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.
  */
 template <typename _DataType>
 INP_DLLEXPORT void dpnp_fft_fft_c(const void* array_in,
@@ -77,6 +78,7 @@ INP_DLLEXPORT void dpnp_fft_fft_c(const void* array_in,
                                   const long* output_shape,
                                   size_t shape_size,
                                   long axis,
-                                  long input_boundarie);
+                                  long input_boundarie,
+                                  size_t inverse);
 
 #endif // BACKEND_IFACE_FFT_H

dpnp/backend/dpnp_kernels_fft.cpp

@@ -42,7 +42,8 @@ void dpnp_fft_fft_c(const void* array1_in,
                     const long* output_shape,
                     size_t shape_size,
                     long axis,
-                    long input_boundarie)
+                    long input_boundarie,
+                    size_t inverse)
 {
     const size_t result_size = std::accumulate(output_shape, output_shape + shape_size, 1, std::multiplies<size_t>());
     if (!(result_size && shape_size))
@@ -51,6 +52,7 @@ void dpnp_fft_fft_c(const void* array1_in,
     }
 
     cl::sycl::event event;
+    const double kernel_pi = inverse ? -M_PI : M_PI;
 
     const _DataType_input* array_1 = reinterpret_cast<const _DataType_input*>(array1_in);
     _DataType_output* result = reinterpret_cast<_DataType_output*>(result1);
@@ -107,7 +109,7 @@ void dpnp_fft_fft_c(const void* array1_in,
             }
 
             const size_t output_local_id = xyz_thread[axis];
-            const double angle = 2.0 * M_PI * it * output_local_id / axis_length;
+            const double angle = 2.0 * kernel_pi * it * output_local_id / axis_length;
 
             const double angle_cos = cl::sycl::cos(angle);
             const double angle_sin = cl::sycl::sin(angle);
@@ -116,6 +118,12 @@ void dpnp_fft_fft_c(const void* array1_in,
             sum_imag += -in_real * angle_sin + in_imag * angle_cos;
         }
 
+        if (inverse)
+        {
+            sum_real = sum_real / input_boundarie;
+            sum_imag = sum_imag / input_boundarie;
+        }
+
         result[output_id] = _DataType_output(sum_real, sum_imag);
     };
 

dpnp/fft/dpnp_algo_fft.pyx

@@ -42,10 +42,10 @@ __all__ = [
     "dpnp_fft"
 ]
 
-ctypedef void(*fptr_dpnp_fft_fft_t)(void *, void * , long * , long * , size_t, long, long)
+ctypedef void(*fptr_dpnp_fft_fft_t)(void *, void * , long * , long * , size_t, long, long, size_t)
 
 
-cpdef dparray dpnp_fft(dparray input, size_t input_boundarie, size_t output_boundarie, long axis):
+cpdef dparray dpnp_fft(dparray input, size_t input_boundarie, size_t output_boundarie, long axis, size_t inverse):
 
     cdef dparray_shape_type input_shape = input.shape
     cdef dparray_shape_type output_shape = input_shape
@@ -66,6 +66,6 @@ cpdef dparray dpnp_fft(dparray input, size_t input_boundarie, size_t output_boun
     cdef fptr_dpnp_fft_fft_t func = <fptr_dpnp_fft_fft_t > 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)
+         output_shape.data(), input_shape.size(), axis_norm, input_boundarie, inverse)
 
     return result

dpnp/fft/dpnp_iface_fft.py

@@ -100,7 +100,7 @@ def fft(x1, n=None, axis=-1, norm=None):
         else:
             output_boundarie = input_boundarie
 
-            return dpnp_fft(x1, input_boundarie, output_boundarie, axis_param)
+            return dpnp_fft(x1, input_boundarie, output_boundarie, axis_param, False)
 
     return call_origin(numpy.fft.fft, x1, n, axis, norm)
 
@@ -197,7 +197,7 @@ def ifft(x1, n=None, axis=-1, norm=None):
 
     is_x1_dparray = isinstance(x1, dparray)
 
-    if (not use_origin_backend(x1) and is_x1_dparray and 0):
+    if (not use_origin_backend(x1) and is_x1_dparray):
         if axis is None:
             axis_param = -1      # the most right dimension (default value)
         else:
@@ -217,7 +217,7 @@ def ifft(x1, n=None, axis=-1, norm=None):
         else:
             output_boundarie = input_boundarie
 
-            return dpnp_fft(x1, input_boundarie, output_boundarie, axis_param)
+            return dpnp_fft(x1, input_boundarie, output_boundarie, axis_param, True)
 
     return call_origin(numpy.fft.ifft, x1, n, axis, norm)
 
@@ -240,11 +240,11 @@ def ifft2(x1, s=None, axes=(-2, -1), norm=None):
 
     is_x1_dparray = isinstance(x1, dparray)
 
-    if (not use_origin_backend(x1) and is_x1_dparray and 0):
+    if (not use_origin_backend(x1) and is_x1_dparray):
         if norm is not None:
             pass
         else:
-            return fftn(x1, s, axes, norm)
+            return ifftn(x1, s, axes, norm)
 
     return call_origin(numpy.fft.ifft2, x1, s, axes, norm)
 
@@ -267,7 +267,7 @@ def ifftn(x1, s=None, axes=None, norm=None):
 
     is_x1_dparray = isinstance(x1, dparray)
 
-    if (not use_origin_backend(x1) and is_x1_dparray and 0):
+    if (not use_origin_backend(x1) and is_x1_dparray):
         if s is None:
             boundaries = tuple([x1.shape[i] for i in range(x1.ndim)])
         else:
@@ -291,7 +291,7 @@ def ifftn(x1, s=None, axes=None, norm=None):
                 except IndexError:
                     checker_throw_axis_error("fft.ifftn", "is out of bounds", param_axis, f"< {len(boundaries)}")
 
-                x1_iter = fft(x1_iter, n=param_n, axis=param_axis, norm=norm)
+                x1_iter = ifft(x1_iter, n=param_n, axis=param_axis, norm=norm)
 
             return x1_iter
 
@@ -332,9 +332,13 @@ def irfft(x1, n=None, axis=-1, norm=None):
         elif norm is not None:
             pass
         else:
-            output_boundarie = input_boundarie
+            output_boundarie = 2 * (input_boundarie - 1)
 
-            return dpnp_fft(x1, input_boundarie, output_boundarie, axis_param)
+            result = dpnp_fft(x1, input_boundarie, output_boundarie, axis_param, True)
+            tmp = dparray(result.shape, dtype=dpnp.float64)
+            for it in range(tmp.size):
+                tmp[it] = result[it].real
+            return tmp
 
     return call_origin(numpy.fft.irfft, x1, n, axis, norm)
 
@@ -357,11 +361,11 @@ def irfft2(x1, s=None, axes=(-2, -1), norm=None):
 
     is_x1_dparray = isinstance(x1, dparray)
 
-    if (not use_origin_backend(x1) and is_x1_dparray and 0):
+    if (not use_origin_backend(x1) and is_x1_dparray):
         if norm is not None:
             pass
         else:
-            return fftn(x1, s, axes, norm)
+            return irfftn(x1, s, axes, norm)
 
     return call_origin(numpy.fft.irfft2, x1, s, axes, norm)
 
@@ -408,7 +412,7 @@ def irfftn(x1, s=None, axes=None, norm=None):
                 except IndexError:
                     checker_throw_axis_error("fft.irfftn", "is out of bounds", param_axis, f"< {len(boundaries)}")
 
-                x1_iter = fft(x1_iter, n=param_n, axis=param_axis, norm=norm)
+                x1_iter = irfft(x1_iter, n=param_n, axis=param_axis, norm=norm)
 
             return x1_iter
 
@@ -451,7 +455,7 @@ def rfft(x1, n=None, axis=-1, norm=None):
         else:
             output_boundarie = input_boundarie // 2 + 1  # rfft specific requirenment
 
-            return dpnp_fft(x1, input_boundarie, output_boundarie, axis_param)
+            return dpnp_fft(x1, input_boundarie, output_boundarie, axis_param, False)
 
     return call_origin(numpy.fft.rfft, x1, n, axis, norm)