26 #if defined (HAVE_CONFIG_H)
30 #if defined (HAVE_FFTW3_H)
40 #if defined (HAVE_FFTW3_THREADS) || defined (HAVE_FFTW3F_THREADS)
46 #if defined (HAVE_FFTW)
69 : m_meth (ESTIMATE), m_rplan (nullptr), m_rd (0), m_rs (0), m_rr (0),
70 m_rh (0), m_rn (), m_rsimd_align (false), m_nthreads (1)
72 m_plan[0] = m_plan[1] =
nullptr;
73 m_d[0] = m_d[1] = m_s[0] = m_s[1] = m_r[0] = m_r[1] = m_h[0] = m_h[1] = 0;
74 m_simd_align[0] = m_simd_align[1] =
false;
75 m_inplace[0] = m_inplace[1] =
false;
78 #if defined (HAVE_FFTW3_THREADS)
79 int init_ret = fftw_init_threads ();
81 (*current_liboctave_error_handler) (
"Error initializing FFTW threads");
93 fftw_plan_with_nthreads (m_nthreads);
97 fftw_import_system_wisdom ();
104 plan_p =
reinterpret_cast<fftw_plan *
> (&m_rplan);
106 fftw_destroy_plan (*plan_p);
108 plan_p =
reinterpret_cast<fftw_plan *
> (&m_plan[0]);
110 fftw_destroy_plan (*plan_p);
112 plan_p =
reinterpret_cast<fftw_plan *
> (&m_plan[1]);
114 fftw_destroy_plan (*plan_p);
134 #if defined (HAVE_FFTW3_THREADS)
137 s_instance->m_nthreads = nt;
138 fftw_plan_with_nthreads (nt);
140 s_instance->m_rplan =
nullptr;
141 s_instance->m_plan[0] = s_instance->m_plan[1] =
nullptr;
144 octave_unused_parameter (nt);
146 (*current_liboctave_warning_handler)
147 (
"unable to change number of threads without FFTW thread support");
151 #define CHECK_SIMD_ALIGNMENT(x) \
152 (((reinterpret_cast<std::ptrdiff_t> (x)) & 0xF) == 0)
155 fftw_planner::do_create_plan (
int dir,
const int rank,
162 int which = (dir == FFTW_FORWARD) ? 0 : 1;
163 fftw_plan *cur_plan_p =
reinterpret_cast<fftw_plan *
> (&m_plan[which]);
164 bool create_new_plan =
false;
166 bool ioinplace = (in == out);
172 if (m_plan[which] ==
nullptr || m_d[which] != dist || m_s[which] != stride
173 || m_r[which] != rank || m_h[which] != howmany
174 || ioinplace != m_inplace[which]
175 || ((ioalign != m_simd_align[which]) ? ! ioalign :
false))
176 create_new_plan =
true;
181 for (
int i = 0; i < rank; i++)
182 if (dims(i) != m_n[which](i))
184 create_new_plan =
true;
194 m_h[which] = howmany;
195 m_simd_align[which] = ioalign;
196 m_inplace[which] = ioinplace;
203 for (
int i = 0, j = rank-1; i < rank; i++, j--)
210 bool plan_destroys_in =
true;
216 plan_flags |= FFTW_ESTIMATE;
217 plan_destroys_in =
false;
220 plan_flags |= FFTW_MEASURE;
223 plan_flags |= FFTW_PATIENT;
226 plan_flags |= FFTW_EXHAUSTIVE;
230 plan_flags |= FFTW_MEASURE;
233 plan_flags |= FFTW_ESTIMATE;
234 plan_destroys_in =
false;
240 plan_flags &= ~FFTW_UNALIGNED;
242 plan_flags |= FFTW_UNALIGNED;
245 fftw_destroy_plan (*cur_plan_p);
248 itmp =
const_cast<Complex *
> (in);
251 if (plan_destroys_in)
255 itmp =
reinterpret_cast<Complex *
>
256 (((
reinterpret_cast<std::ptrdiff_t
> (itmp) + 15) & ~ 0xF)
257 + ((
reinterpret_cast<std::ptrdiff_t
> (in)) & 0xF));
264 = fftw_plan_many_dft (rank, tmp, howmany,
265 reinterpret_cast<fftw_complex *
> (itmp),
266 nullptr, stride, dist,
267 reinterpret_cast<fftw_complex *
> (otmp),
268 nullptr, stride, dist, dir, plan_flags);
270 if (*cur_plan_p ==
nullptr)
271 (*current_liboctave_error_handler) (
"Error creating FFTW plan");
278 fftw_planner::do_create_plan (
const int rank,
const dim_vector& dims,
282 const double *in,
Complex *out)
284 fftw_plan *cur_plan_p =
reinterpret_cast<fftw_plan *
> (&m_rplan);
285 bool create_new_plan =
false;
292 if (m_rplan ==
nullptr || m_rd != dist || m_rs != stride || m_rr != rank
293 || m_rh != howmany || ((ioalign != m_rsimd_align) ? ! ioalign :
false))
294 create_new_plan =
true;
299 for (
int i = 0; i < rank; i++)
300 if (dims(i) != m_rn(i))
302 create_new_plan =
true;
313 m_rsimd_align = ioalign;
320 for (
int i = 0, j = rank-1; i < rank; i++, j--)
327 bool plan_destroys_in =
true;
333 plan_flags |= FFTW_ESTIMATE;
334 plan_destroys_in =
false;
337 plan_flags |= FFTW_MEASURE;
340 plan_flags |= FFTW_PATIENT;
343 plan_flags |= FFTW_EXHAUSTIVE;
347 plan_flags |= FFTW_MEASURE;
350 plan_flags |= FFTW_ESTIMATE;
351 plan_destroys_in =
false;
357 plan_flags &= ~FFTW_UNALIGNED;
359 plan_flags |= FFTW_UNALIGNED;
362 fftw_destroy_plan (*cur_plan_p);
365 itmp =
const_cast<double *
> (in);
368 if (plan_destroys_in)
373 nn * howmany * (in_place + 1) + 32);
374 itmp =
reinterpret_cast<double *
>
375 (((
reinterpret_cast<std::ptrdiff_t
> (itmp) + 15) & ~ 0xF)
376 + ((
reinterpret_cast<std::ptrdiff_t
> (in)) & 0xF));
379 otmp =
reinterpret_cast<Complex *
> (itmp);
383 = fftw_plan_many_dft_r2c (rank, tmp, howmany, itmp,
384 nullptr, stride, dist,
385 reinterpret_cast<fftw_complex *
> (otmp),
386 nullptr, stride, dist, plan_flags);
388 if (*cur_plan_p ==
nullptr)
389 (*current_liboctave_error_handler) (
"Error creating FFTW plan");
396 fftw_planner::do_method ()
402 fftw_planner::do_method (FftwMethod _meth)
413 fftw_destroy_plan (
reinterpret_cast<fftw_plan
> (m_rplan));
415 fftw_destroy_plan (
reinterpret_cast<fftw_plan
> (m_plan[0]));
417 fftw_destroy_plan (
reinterpret_cast<fftw_plan
> (m_plan[1]));
418 m_rplan = m_plan[0] = m_plan[1] =
nullptr;
429 : m_meth (ESTIMATE), m_rplan (nullptr), m_rd (0), m_rs (0), m_rr (0),
430 m_rh (0), m_rn (), m_rsimd_align (false), m_nthreads (1)
432 m_plan[0] = m_plan[1] =
nullptr;
433 m_d[0] = m_d[1] = m_s[0] = m_s[1] = m_r[0] = m_r[1] = m_h[0] = m_h[1] = 0;
434 m_simd_align[0] = m_simd_align[1] =
false;
435 m_inplace[0] = m_inplace[1] =
false;
438 #if defined (HAVE_FFTW3F_THREADS)
439 int init_ret = fftwf_init_threads ();
441 (*current_liboctave_error_handler) (
"Error initializing FFTW3F threads");
448 fftwf_plan_with_nthreads (m_nthreads);
452 fftwf_import_system_wisdom ();
459 plan_p =
reinterpret_cast<fftwf_plan *
> (&m_rplan);
461 fftwf_destroy_plan (*plan_p);
463 plan_p =
reinterpret_cast<fftwf_plan *
> (&m_plan[0]);
465 fftwf_destroy_plan (*plan_p);
467 plan_p =
reinterpret_cast<fftwf_plan *
> (&m_plan[1]);
469 fftwf_destroy_plan (*plan_p);
489 #if defined (HAVE_FFTW3F_THREADS)
492 s_instance->m_nthreads = nt;
493 fftwf_plan_with_nthreads (nt);
495 s_instance->m_rplan =
nullptr;
496 s_instance->m_plan[0] = s_instance->m_plan[1] =
nullptr;
499 octave_unused_parameter (nt);
501 (*current_liboctave_warning_handler)
502 (
"unable to change number of threads without FFTW thread support");
507 float_fftw_planner::do_create_plan (
int dir,
const int rank,
515 int which = (dir == FFTW_FORWARD) ? 0 : 1;
516 fftwf_plan *cur_plan_p =
reinterpret_cast<fftwf_plan *
> (&m_plan[which]);
517 bool create_new_plan =
false;
519 bool ioinplace = (in == out);
525 if (m_plan[which] ==
nullptr || m_d[which] != dist || m_s[which] != stride
526 || m_r[which] != rank || m_h[which] != howmany
527 || ioinplace != m_inplace[which]
528 || ((ioalign != m_simd_align[which]) ? ! ioalign :
false))
529 create_new_plan =
true;
534 for (
int i = 0; i < rank; i++)
535 if (dims(i) != m_n[which](i))
537 create_new_plan =
true;
547 m_h[which] = howmany;
548 m_simd_align[which] = ioalign;
549 m_inplace[which] = ioinplace;
556 for (
int i = 0, j = rank-1; i < rank; i++, j--)
563 bool plan_destroys_in =
true;
569 plan_flags |= FFTW_ESTIMATE;
570 plan_destroys_in =
false;
573 plan_flags |= FFTW_MEASURE;
576 plan_flags |= FFTW_PATIENT;
579 plan_flags |= FFTW_EXHAUSTIVE;
583 plan_flags |= FFTW_MEASURE;
586 plan_flags |= FFTW_ESTIMATE;
587 plan_destroys_in =
false;
593 plan_flags &= ~FFTW_UNALIGNED;
595 plan_flags |= FFTW_UNALIGNED;
598 fftwf_destroy_plan (*cur_plan_p);
604 if (plan_destroys_in)
609 (((
reinterpret_cast<std::ptrdiff_t
> (itmp) + 15) & ~ 0xF)
610 + ((
reinterpret_cast<std::ptrdiff_t
> (in)) & 0xF));
617 = fftwf_plan_many_dft (rank, tmp, howmany,
618 reinterpret_cast<fftwf_complex *
> (itmp),
619 nullptr, stride, dist,
620 reinterpret_cast<fftwf_complex *
> (otmp),
621 nullptr, stride, dist, dir, plan_flags);
623 if (*cur_plan_p ==
nullptr)
624 (*current_liboctave_error_handler) (
"Error creating FFTW plan");
631 float_fftw_planner::do_create_plan (
const int rank,
const dim_vector& dims,
637 fftwf_plan *cur_plan_p =
reinterpret_cast<fftwf_plan *
> (&m_rplan);
638 bool create_new_plan =
false;
645 if (m_rplan ==
nullptr || m_rd != dist || m_rs != stride || m_rr != rank
646 || m_rh != howmany || ((ioalign != m_rsimd_align) ? ! ioalign :
false))
647 create_new_plan =
true;
652 for (
int i = 0; i < rank; i++)
653 if (dims(i) != m_rn(i))
655 create_new_plan =
true;
666 m_rsimd_align = ioalign;
673 for (
int i = 0, j = rank-1; i < rank; i++, j--)
680 bool plan_destroys_in =
true;
686 plan_flags |= FFTW_ESTIMATE;
687 plan_destroys_in =
false;
690 plan_flags |= FFTW_MEASURE;
693 plan_flags |= FFTW_PATIENT;
696 plan_flags |= FFTW_EXHAUSTIVE;
700 plan_flags |= FFTW_MEASURE;
703 plan_flags |= FFTW_ESTIMATE;
704 plan_destroys_in =
false;
710 plan_flags &= ~FFTW_UNALIGNED;
712 plan_flags |= FFTW_UNALIGNED;
715 fftwf_destroy_plan (*cur_plan_p);
718 itmp =
const_cast<float *
> (in);
721 if (plan_destroys_in)
726 nn * howmany * (in_place + 1) + 32);
727 itmp =
reinterpret_cast<float *
>
728 (((
reinterpret_cast<std::ptrdiff_t
> (itmp) + 15) & ~ 0xF)
729 + ((
reinterpret_cast<std::ptrdiff_t
> (in)) & 0xF));
736 = fftwf_plan_many_dft_r2c (rank, tmp, howmany, itmp,
737 nullptr, stride, dist,
738 reinterpret_cast<fftwf_complex *
> (otmp),
739 nullptr, stride, dist, plan_flags);
741 if (*cur_plan_p ==
nullptr)
742 (*current_liboctave_error_handler) (
"Error creating FFTW plan");
749 float_fftw_planner::do_method ()
755 float_fftw_planner::do_method (FftwMethod _meth)
766 fftwf_destroy_plan (
reinterpret_cast<fftwf_plan
> (m_rplan));
768 fftwf_destroy_plan (
reinterpret_cast<fftwf_plan
> (m_plan[0]));
770 fftwf_destroy_plan (
reinterpret_cast<fftwf_plan
> (m_plan[1]));
771 m_rplan = m_plan[0] = m_plan[1] =
nullptr;
779 template <
typename T>
781 convert_packcomplex_1d (T *out, std::size_t nr, std::size_t nc,
788 for (std::size_t i = 0; i < nr; i++)
789 for (std::size_t j = nc/2+1; j < nc; j++)
790 out[j*stride + i*dist] =
conj (out[(nc - j)*stride + i*dist]);
795 template <
typename T>
797 convert_packcomplex_Nd (T *out,
const dim_vector& dv)
799 std::size_t nc = dv(0);
800 std::size_t nr = dv(1);
801 std::size_t np = (dv.
ndims () > 2 ? dv.
numel () / nc / nr : 1);
802 std::size_t nrp = nr * np;
809 for (std::size_t i = 0; i < nrp; i++)
811 ptr1 = out + i * (nc/2 + 1) + nrp*((nc-1)/2);
813 for (std::size_t j = 0; j < nc/2+1; j++)
821 for (std::size_t i = 0; i < np; i++)
823 for (std::size_t j = 1; j < nr; j++)
824 for (std::size_t k = nc/2+1; k < nc; k++)
825 out[k + (j + i*nr)*nc] =
conj (out[nc - k + ((i+1)*nr - j)*nc]);
827 for (std::size_t j = nc/2+1; j < nc; j++)
828 out[j + i*nr*nc] =
conj (out[(i*nr+1)*nc - j]);
835 std::size_t jstart = dv(0) * dv(1);
836 std::size_t kstep = dv(0);
837 std::size_t nel = dv.
numel ();
839 for (
int inner = 2; inner < dv.
ndims (); inner++)
841 std::size_t jmax = jstart * dv(inner);
842 for (std::size_t i = 0; i < nel; i+=jmax)
843 for (std::size_t j = jstart, jj = jmax-jstart; j < jj;
844 j+=jstart, jj-=jstart)
845 for (std::size_t k = 0; k < jstart; k+= kstep)
846 for (std::size_t l = nc/2+1; l < nc; l++)
848 T tmp = out[i+ j + k + l];
849 out[i + j + k + l] = out[i + jj + k + l];
850 out[i + jj + k + l] = tmp;
863 dist = (dist < 0 ? npts : dist);
867 stride, dist, in, out);
868 fftw_plan m_plan =
reinterpret_cast<fftw_plan
> (vplan);
870 fftw_execute_dft_r2c (m_plan, (
const_cast<double *
> (in)),
871 reinterpret_cast<fftw_complex *
> (out));
875 convert_packcomplex_1d (out, nsamples, npts, stride, dist);
885 dist = (dist < 0 ? npts : dist);
891 fftw_plan m_plan =
reinterpret_cast<fftw_plan
> (vplan);
893 fftw_execute_dft (m_plan,
894 reinterpret_cast<fftw_complex *
> (
const_cast<Complex *
> (in)),
895 reinterpret_cast<fftw_complex *
> (out));
905 dist = (dist < 0 ? npts : dist);
909 stride, dist, in, out);
910 fftw_plan m_plan =
reinterpret_cast<fftw_plan
> (vplan);
912 fftw_execute_dft (m_plan,
913 reinterpret_cast<fftw_complex *
> (
const_cast<Complex *
> (in)),
914 reinterpret_cast<fftw_complex *
> (out));
917 for (std::size_t j = 0; j < nsamples; j++)
918 for (std::size_t i = 0; i < npts; i++)
919 out[i*stride + j*dist] /=
scale;
929 for (
int i = 0; i < rank; i++)
939 fftw_plan m_plan =
reinterpret_cast<fftw_plan
> (vplan);
941 fftw_execute_dft_r2c (m_plan, (
const_cast<double *
> (in)),
942 reinterpret_cast<fftw_complex *
> (out+ offset));
946 convert_packcomplex_Nd (out, dv);
956 for (
int i = 0; i < rank; i++)
960 1, 1, dist, in, out);
961 fftw_plan m_plan =
reinterpret_cast<fftw_plan
> (vplan);
963 fftw_execute_dft (m_plan,
964 reinterpret_cast<fftw_complex *
> (
const_cast<Complex *
> (in)),
965 reinterpret_cast<fftw_complex *
> (out));
975 for (
int i = 0; i < rank; i++)
979 1, 1, dist, in, out);
980 fftw_plan m_plan =
reinterpret_cast<fftw_plan
> (vplan);
982 fftw_execute_dft (m_plan,
983 reinterpret_cast<fftw_complex *
> (
const_cast<Complex *
> (in)),
984 reinterpret_cast<fftw_complex *
> (out));
986 const std::size_t npts = dv.
numel ();
988 for (std::size_t i = 0; i < npts; i++)
999 dist = (dist < 0 ? npts : dist);
1004 fftwf_plan m_plan =
reinterpret_cast<fftwf_plan
> (vplan);
1006 fftwf_execute_dft_r2c (m_plan, (
const_cast<float *
> (in)),
1007 reinterpret_cast<fftwf_complex *
> (out));
1011 convert_packcomplex_1d (out, nsamples, npts, stride, dist);
1021 dist = (dist < 0 ? npts : dist);
1025 nsamples, stride, dist,
1027 fftwf_plan m_plan =
reinterpret_cast<fftwf_plan
> (vplan);
1029 fftwf_execute_dft (m_plan,
1030 reinterpret_cast<fftwf_complex *
> (
const_cast<FloatComplex *
> (in)),
1031 reinterpret_cast<fftwf_complex *
> (out));
1041 dist = (dist < 0 ? npts : dist);
1045 nsamples, stride, dist,
1047 fftwf_plan m_plan =
reinterpret_cast<fftwf_plan
> (vplan);
1049 fftwf_execute_dft (m_plan,
1050 reinterpret_cast<fftwf_complex *
> (
const_cast<FloatComplex *
> (in)),
1051 reinterpret_cast<fftwf_complex *
> (out));
1054 for (std::size_t j = 0; j < nsamples; j++)
1055 for (std::size_t i = 0; i < npts; i++)
1056 out[i*stride + j*dist] /=
scale;
1066 for (
int i = 0; i < rank; i++)
1076 fftwf_plan m_plan =
reinterpret_cast<fftwf_plan
> (vplan);
1078 fftwf_execute_dft_r2c (m_plan, (
const_cast<float *
> (in)),
1079 reinterpret_cast<fftwf_complex *
> (out+ offset));
1083 convert_packcomplex_Nd (out, dv);
1093 for (
int i = 0; i < rank; i++)
1097 1, 1, dist, in, out);
1098 fftwf_plan m_plan =
reinterpret_cast<fftwf_plan
> (vplan);
1100 fftwf_execute_dft (m_plan,
1101 reinterpret_cast<fftwf_complex *
> (
const_cast<FloatComplex *
> (in)),
1102 reinterpret_cast<fftwf_complex *
> (out));
1112 for (
int i = 0; i < rank; i++)
1116 1, 1, dist, in, out);
1117 fftwf_plan m_plan =
reinterpret_cast<fftwf_plan
> (vplan);
1119 fftwf_execute_dft (m_plan,
1120 reinterpret_cast<fftwf_complex *
> (
const_cast<FloatComplex *
> (in)),
1121 reinterpret_cast<fftwf_complex *
> (out));
1123 const std::size_t npts = dv.
numel ();
1125 for (std::size_t i = 0; i < npts; i++)
1136 #if defined (HAVE_FFTW)
1146 #if defined (HAVE_FFTW)
1153 OCTAVE_END_NAMESPACE(
octave)
ComplexColumnVector conj(const ComplexColumnVector &a)
Vector representing the dimensions (size) of an Array.
octave_idx_type numel(int n=0) const
Number of elements that a matrix with this dimensions would have.
octave_idx_type ndims() const
Number of dimensions.
static bool instance_ok()
static void * create_plan(int dir, const int rank, const dim_vector &dims, octave_idx_type howmany, octave_idx_type stride, octave_idx_type dist, const Complex *in, Complex *out)
static int ifftNd(const Complex *, Complex *, const int, const dim_vector &)
static int fft(const double *in, Complex *out, std::size_t npts, std::size_t nsamples=1, octave_idx_type stride=1, octave_idx_type dist=-1)
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 fftNd(const double *, Complex *, const int, const dim_vector &)
static void * create_plan(int dir, const int rank, const dim_vector &dims, octave_idx_type howmany, octave_idx_type stride, octave_idx_type dist, const FloatComplex *in, FloatComplex *out)
static bool instance_ok()
OCTAVE_BEGIN_NAMESPACE(octave) static octave_value daspk_fcn
void scale(Matrix &m, double x, double y, double z)
unsigned long int octave_num_processors_wrapper(enum octave_nproc_query octave_query)
@ OCTAVE_NPROC_CURRENT_OVERRIDABLE
std::complex< double > Complex
std::complex< float > FloatComplex
std::string fftwf_version()
std::string fftw_version()
#define CHECK_SIMD_ALIGNMENT(x)
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
#define OCTAVE_SCOPED_BUFFER_ANCHOR(T, buf)
#define OCTAVE_SCOPED_BUFFER(T, buf, size)