# HG changeset patch
# User Markus Meisinger
# Date 1640293701 -3600
#      Thu Dec 23 22:08:21 2021 +0100
# Node ID f184ebd6e502b42d44bac9b5f75fab89507207c5
# Parent  b65d423b12a6a391867663ec6b51c27c8a6f99bb
Do not push to HG! Runtime measurement: Use fftw_execute_dft_c2r for ComplexNDArray when the symflag is given in ifft().

* fft.cc: Peel symflag.
* CNDArray.cc: Implement ifouriersym() to transport symflag to lower software level.
* CNDArray.h: Implement ifouriersym() to transport symflag to lower software level.
* oct-fftw.cc: FFTW planner for c2r transformation and new ifftsym() method.
* oct-fftw.h: FFTW planner for c2r transformation and new ifftsym() method.

diff -r b65d423b12a6 -r f184ebd6e502 libinterp/corefcn/fft.cc
--- a/libinterp/corefcn/fft.cc	Sun Dec 19 20:42:45 2021 +0100
+++ b/libinterp/corefcn/fft.cc	Thu Dec 23 22:08:21 2021 +0100
@@ -41,8 +41,9 @@
 do_fft (const octave_value_list& args, const char *fcn, int type)
 {
   int nargin = args.length ();
+  int numeric_argin = nargin;
 
-  if (nargin < 1 || nargin > 3)
+  if (nargin < 1 || nargin > 4)
     print_usage ();
 
   octave_value retval;
@@ -51,8 +52,24 @@
   dim_vector dims = arg.dims ();
   int ndims = dims.ndims ();
   int dim = -1;
+  bool have_symmetric = false;
 
-  if (nargin > 1)
+  // Peel "symflag" if present
+  if ((nargin > 1) && (! args(nargin - 1).isempty ()) && (! args(nargin - 1).isnumeric ()))
+    {
+      std::string strval = args(nargin - 1).xstring_value ("must be a string");
+      if (strval == "symmetric")
+        {
+          have_symmetric = true;
+        }
+      else if (strval != "nonsymmetric")
+        {
+          error ("%s: unknown symflag", fcn);
+        }
+      numeric_argin = nargin - 1;
+    }
+
+  if (numeric_argin > 1)
     {
       if (! args(1).isempty ())
         {
@@ -66,7 +83,7 @@
         }
     }
 
-  if (nargin > 2)
+  if (numeric_argin > 2)
     {
       double dval = args(2).double_value ();
       if (math::isnan (dval))
@@ -146,10 +163,14 @@
         }
       else if (arg.iscomplex ())
         {
-          ComplexNDArray cnda = arg.complex_array_value ();
+          ComplexNDArray in = arg.complex_array_value ();
+
+          in.resize (dims, 0.0);
 
-          cnda.resize (dims, 0.0);
-          retval = (type != 0 ? cnda.ifourier (dim) : cnda.fourier (dim));
+          if (type && have_symmetric)
+            retval = in.ifouriersym (dim);
+          else
+            retval = (type != 0 ? in.ifourier (dim) : in.fourier (dim));
         }
       else
         err_wrong_type_arg (fcn, arg);
diff -r b65d423b12a6 -r f184ebd6e502 liboctave/array/CNDArray.cc
--- a/liboctave/array/CNDArray.cc	Sun Dec 19 20:42:45 2021 +0100
+++ b/liboctave/array/CNDArray.cc	Thu Dec 23 22:08:21 2021 +0100
@@ -116,6 +116,37 @@
   return retval;
 }
 
+NDArray
+ComplexNDArray::ifouriersym (int dim) const
+{
+  dim_vector dv = dims ();
+
+  if (dim > dv.ndims () || dim < 0)
+    return NDArray ();
+
+  octave_idx_type stride = 1;
+  octave_idx_type n = dv(dim);
+
+  for (int i = 0; i < dim; i++)
+    stride *= dv(i);
+
+  octave_idx_type howmany = numel () / dv(dim);
+  howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
+  octave_idx_type nloop = (stride == 1 ? 1 : numel () / dv(dim) / stride);
+  octave_idx_type dist = (stride == 1 ? n : 1);
+
+  const Complex *in (data ());
+  NDArray retval (dv);
+  double *out (retval.fortran_vec ());
+
+  // Need to be careful here about the distance between fft's
+  for (octave_idx_type k = 0; k < nloop; k++)
+    octave::fftw::ifftsym (in + k * stride * n, out + k * stride * n,
+                           n, howmany, stride, dist);
+
+  return retval;
+}
+
 ComplexNDArray
 ComplexNDArray::fourier2d (void) const
 {
@@ -210,6 +241,17 @@
   return ComplexNDArray ();
 }
 
+NDArray
+ComplexNDArray::ifouriersym (int dim) const
+{
+  octave_unused_parameter (dim);
+
+  (*current_liboctave_error_handler)
+    ("support for FFTW was unavailable or disabled when liboctave was built");
+
+  return ComplexNDArray ();
+}
+
 ComplexNDArray
 ComplexNDArray::fourier2d (void) const
 {
diff -r b65d423b12a6 -r f184ebd6e502 liboctave/array/CNDArray.h
--- a/liboctave/array/CNDArray.h	Sun Dec 19 20:42:45 2021 +0100
+++ b/liboctave/array/CNDArray.h	Thu Dec 23 22:08:21 2021 +0100
@@ -121,6 +121,7 @@
 
   OCTAVE_API ComplexNDArray fourier (int dim = 1) const;
   OCTAVE_API ComplexNDArray ifourier (int dim = 1) const;
+  OCTAVE_API NDArray ifouriersym (int dim = 1) const;
 
   OCTAVE_API ComplexNDArray fourier2d (void) const;
   OCTAVE_API ComplexNDArray ifourier2d (void) const;
diff -r b65d423b12a6 -r f184ebd6e502 liboctave/numeric/oct-fftw.cc
--- a/liboctave/numeric/oct-fftw.cc	Sun Dec 19 20:42:45 2021 +0100
+++ b/liboctave/numeric/oct-fftw.cc	Thu Dec 23 22:08:21 2021 +0100
@@ -75,6 +75,11 @@
     m_inplace[0] = m_inplace[1] = false;
     m_n[0] = m_n[1] = dim_vector ();
 
+    m_cplan = nullptr;
+    m_cd = m_cs = m_cr = m_ch = 0;
+    m_cn = dim_vector ();
+    m_csimd_align = false;
+
 #if defined (HAVE_FFTW3_THREADS)
     int init_ret = fftw_init_threads ();
     if (! init_ret)
@@ -100,6 +105,10 @@
     if (*plan_p)
       fftw_destroy_plan (*plan_p);
 
+    plan_p = reinterpret_cast<fftw_plan *> (&m_cplan);
+    if (*plan_p)
+      fftw_destroy_plan (*plan_p);
+
     plan_p = reinterpret_cast<fftw_plan *> (&m_plan[0]);
     if (*plan_p)
       fftw_destroy_plan (*plan_p);
@@ -133,6 +142,7 @@
         fftw_plan_with_nthreads (nt);
         // Clear the current plans.
         s_instance->m_rplan = nullptr;
+        s_instance->m_cplan = nullptr;
         s_instance->m_plan[0] = s_instance->m_plan[1] = nullptr;
       }
 #else
@@ -389,6 +399,126 @@
     return *cur_plan_p;
   }
 
+  void *
+  fftw_planner::do_create_plan (const int rank, const dim_vector& dims,
+                                octave_idx_type howmany,
+                                octave_idx_type stride,
+                                octave_idx_type dist,
+                                const Complex *in, double *out)
+  {
+    fftw_plan *cur_plan_p = reinterpret_cast<fftw_plan *> (&m_cplan);
+    bool create_new_plan = false;
+    bool ioalign = CHECK_SIMD_ALIGNMENT (in) && CHECK_SIMD_ALIGNMENT (out);
+
+    // Don't create a new plan if we have a non SIMD plan already but
+    // can do SIMD.  This prevents endlessly recreating plans if we
+    // change the alignment.
+
+    if (m_cplan == nullptr || m_cd != dist || m_cs != stride || m_cr != rank
+        || m_ch != howmany || ((ioalign != m_csimd_align) ? ! ioalign : false))
+      create_new_plan = true;
+    else
+      {
+        // We still might not have the same shape of array.
+
+        for (int i = 0; i < rank; i++)
+          if (dims(i) != m_cn(i))
+            {
+              create_new_plan = true;
+              break;
+            }
+      }
+
+    if (create_new_plan)
+      {
+        m_cd = dist;
+        m_cs = stride;
+        m_cr = rank;
+        m_ch = howmany;
+        m_csimd_align = ioalign;
+        m_cn = dims;
+
+        // Note reversal of dimensions for column major storage in FFTW.
+        octave_idx_type nn = 1;
+        OCTAVE_LOCAL_BUFFER (int, tmp, rank);
+
+        for (int i = 0, j = rank-1; i < rank; i++, j--)
+          {
+            tmp[i] = dims(j);
+            nn *= dims(j);
+          }
+
+        int plan_flags = FFTW_PRESERVE_INPUT;
+        bool plan_destroys_in = true;
+
+        switch (m_meth)
+          {
+          case UNKNOWN:
+          case ESTIMATE:
+            plan_flags |= FFTW_ESTIMATE;
+            plan_destroys_in = false;
+            break;
+          case MEASURE:
+            plan_flags |= FFTW_MEASURE;
+            break;
+          case PATIENT:
+            plan_flags |= FFTW_PATIENT;
+            break;
+          case EXHAUSTIVE:
+            plan_flags |= FFTW_EXHAUSTIVE;
+            break;
+          case HYBRID:
+            if (nn < 8193)
+              plan_flags |= FFTW_MEASURE;
+            else
+              {
+                plan_flags |= FFTW_ESTIMATE;
+                plan_destroys_in = false;
+              }
+            break;
+          }
+
+        if (ioalign)
+          plan_flags &= ~FFTW_UNALIGNED;
+        else
+          plan_flags |= FFTW_UNALIGNED;
+
+        if (*cur_plan_p)
+          fftw_destroy_plan (*cur_plan_p);
+
+        if (plan_destroys_in)
+          {
+            // Create matrix with the same size and 16-byte alignment as input
+            OCTAVE_LOCAL_BUFFER (Complex, itmp, nn * howmany + 32);
+            itmp = reinterpret_cast<Complex *>
+                   (((reinterpret_cast<std::ptrdiff_t> (itmp) + 15) & ~ 0xF) +
+                    ((reinterpret_cast<std::ptrdiff_t> (in)) & 0xF));
+
+            *cur_plan_p
+              = fftw_plan_many_dft_c2r (rank, tmp, howmany,
+                                        reinterpret_cast<fftw_complex *> (itmp),
+                                        nullptr, stride, dist,
+                                        out,
+                                        nullptr, stride, dist, plan_flags);
+          }
+        else
+          {
+            *cur_plan_p
+              = fftw_plan_many_dft_c2r (rank, tmp, howmany,
+                                        reinterpret_cast<fftw_complex *> (const_cast<Complex *> (in)),
+                                        nullptr, stride, dist,
+                                        out,
+                                        nullptr, stride, dist, plan_flags);
+          }
+
+        if (*cur_plan_p == nullptr)
+          (*current_liboctave_error_handler) ("Error creating FFTW plan");
+      }
+
+    return *cur_plan_p;
+  }
+
+
   fftw_planner::FftwMethod
   fftw_planner::do_method (void)
   {
@@ -408,11 +538,13 @@
             m_meth = _meth;
             if (m_rplan)
               fftw_destroy_plan (reinterpret_cast<fftw_plan> (m_rplan));
+            if (m_cplan)
+              fftw_destroy_plan (reinterpret_cast<fftw_plan> (m_cplan));
             if (m_plan[0])
               fftw_destroy_plan (reinterpret_cast<fftw_plan> (m_plan[0]));
             if (m_plan[1])
               fftw_destroy_plan (reinterpret_cast<fftw_plan> (m_plan[1]));
-            m_rplan = m_plan[0] = m_plan[1] = nullptr;
+            m_rplan = m_cplan = m_plan[0] = m_plan[1] = nullptr;
           }
       }
     else
@@ -921,6 +1053,30 @@
   }
 
   int
+  fftw::ifftsym (const Complex *in, double *out, std::size_t npts,
+                 std::size_t nsamples, octave_idx_type stride,
+                 octave_idx_type dist)
+  {
+    dist = (dist < 0 ? npts : dist);
+
+    dim_vector dv (npts, 1);
+    void *vplan = fftw_planner::create_plan (1, dv, nsamples,
+                                             stride, dist, in, out);
+    fftw_plan m_plan = reinterpret_cast<fftw_plan> (vplan);
+
+    fftw_execute_dft_c2r (m_plan,
+                          reinterpret_cast<fftw_complex *> (const_cast<Complex *>(in)),
+                          reinterpret_cast<double *>(out));
+
+    const double scale = npts;
+    for (std::size_t j = 0; j < nsamples; j++)
+      for (std::size_t i = 0; i < npts; i++)
+        out[i*stride + j*dist] /= scale;
+
+    return 0;
+  }
+
+  int
   fftw::fftNd (const double *in, Complex *out, const int rank,
                const dim_vector& dv)
   {
diff -r b65d423b12a6 -r f184ebd6e502 liboctave/numeric/oct-fftw.h
--- a/liboctave/numeric/oct-fftw.h	Sun Dec 19 20:42:45 2021 +0100
+++ b/liboctave/numeric/oct-fftw.h	Thu Dec 23 22:08:21 2021 +0100
@@ -90,6 +90,17 @@
         : nullptr;
     }
 
+    static void *
+    create_plan (const int rank, const dim_vector& dims,
+                 octave_idx_type howmany, octave_idx_type stride,
+                 octave_idx_type dist, const Complex *in, double *out)
+    {
+      return instance_ok ()
+        ? s_instance->do_create_plan (rank, dims, howmany, stride, dist,
+                                      in, out)
+        : nullptr;
+    }
+
     static FftwMethod method (void)
     {
       static FftwMethod dummy;
@@ -129,6 +140,11 @@
                     octave_idx_type howmany, octave_idx_type stride,
                     octave_idx_type dist, const double *in, Complex *out);
 
+    void *
+    do_create_plan (const int rank, const dim_vector& dims,
+                    octave_idx_type howmany, octave_idx_type stride,
+                    octave_idx_type dist, const Complex *in, double *out);
+
     FftwMethod do_method (void);
 
     FftwMethod do_method (FftwMethod meth);
@@ -158,6 +174,26 @@
     bool m_simd_align[2];
     bool m_inplace[2];
 
+    // Plan for ifft symmetric
+    void *m_cplan;
+
+    // dist
+    octave_idx_type m_cd;
+
+    // stride
+    octave_idx_type m_cs;
+
+    // rank
+    int m_cr;
+
+    // howmany
+    octave_idx_type m_ch;
+
+    // dims
+    dim_vector m_cn;
+
+    bool m_csimd_align;
+
     // Plan for fft of real values
     void *m_rplan;
 
@@ -353,6 +389,9 @@
     static int ifft (const Complex *in, Complex *out, std::size_t npts,
                      std::size_t nsamples = 1, octave_idx_type stride = 1,
                      octave_idx_type dist = -1);
+    static int ifftsym (const Complex *in, double *out, std::size_t npts,
+                        std::size_t nsamples = 1, octave_idx_type stride = 1,
+                        octave_idx_type dist = -1);
 
     static int fftNd (const double *, Complex *, const int, const dim_vector&);
     static int fftNd (const Complex *, Complex *, const int,
