23 #if defined (HAVE_CONFIG_H) 27 #if defined (HAVE_FFTW3_H) 37 #if defined (HAVE_FFTW3_THREADS) || defined (HAVE_FFTW3F_THREADS) 43 #if defined (HAVE_FFTW) 66 : meth (ESTIMATE), rplan (nullptr), rd (0), rs (0), rr (0), rh (0), rn (),
67 rsimd_align (
false), nthreads (1)
70 d[0] =
d[1] =
s[0] =
s[1] =
r[0] =
r[1] =
h[0] =
h[1] = 0;
75 #if defined (HAVE_FFTW3_THREADS) 76 int init_ret = fftw_init_threads ();
78 (*current_liboctave_error_handler) (
"Error initializing FFTW threads");
87 fftw_import_system_wisdom ();
94 plan_p =
reinterpret_cast<fftw_plan *
> (&
rplan);
96 fftw_destroy_plan (*plan_p);
98 plan_p =
reinterpret_cast<fftw_plan *
> (&
plan[0]);
100 fftw_destroy_plan (*plan_p);
102 plan_p =
reinterpret_cast<fftw_plan *
> (&
plan[1]);
104 fftw_destroy_plan (*plan_p);
121 (*current_liboctave_error_handler)
122 (
"unable to create fftw_planner object!");
130 #if defined (HAVE_FFTW3_THREADS) 134 fftw_plan_with_nthreads (nt);
139 (*current_liboctave_warning_handler)
140 (
"unable to change number of threads without FFTW thread support");
144 #define CHECK_SIMD_ALIGNMENT(x) \ 145 (((reinterpret_cast<ptrdiff_t> (x)) & 0xF) == 0) 155 int which = (dir == FFTW_FORWARD) ? 0 : 1;
156 fftw_plan *cur_plan_p =
reinterpret_cast<fftw_plan *
> (&
plan[which]);
157 bool create_new_plan =
false;
159 bool ioinplace = (in == out);
165 if (
plan[which] ==
nullptr ||
d[which] != dist ||
s[which] != stride
166 ||
r[which] != rank ||
h[which] != howmany
168 || ((ioalign !=
simd_align[which]) ? ! ioalign :
false))
169 create_new_plan =
true;
174 for (
int i = 0;
i < rank;
i++)
177 create_new_plan =
true;
196 for (
int i = 0, j = rank-1;
i < rank;
i++, j--)
203 bool plan_destroys_in =
true;
209 plan_flags |= FFTW_ESTIMATE;
210 plan_destroys_in =
false;
213 plan_flags |= FFTW_MEASURE;
216 plan_flags |= FFTW_PATIENT;
219 plan_flags |= FFTW_EXHAUSTIVE;
223 plan_flags |= FFTW_MEASURE;
226 plan_flags |= FFTW_ESTIMATE;
227 plan_destroys_in =
false;
233 plan_flags &= ~FFTW_UNALIGNED;
235 plan_flags |= FFTW_UNALIGNED;
238 fftw_destroy_plan (*cur_plan_p);
240 if (plan_destroys_in)
244 itmp =
reinterpret_cast<Complex *
> 245 (((
reinterpret_cast<ptrdiff_t
>(itmp) + 15) & ~ 0xF) +
246 ((reinterpret_cast<ptrdiff_t> (in)) & 0xF));
249 fftw_plan_many_dft (rank,
tmp, howmany,
250 reinterpret_cast<fftw_complex *> (itmp),
251 nullptr, stride, dist,
252 reinterpret_cast<fftw_complex *> (out),
253 nullptr, stride, dist, dir, plan_flags);
258 fftw_plan_many_dft (rank,
tmp, howmany,
259 reinterpret_cast<fftw_complex *> (const_cast<Complex *> (in)),
260 nullptr, stride, dist,
261 reinterpret_cast<fftw_complex *> (out),
262 nullptr, stride, dist, dir, plan_flags);
265 if (*cur_plan_p ==
nullptr)
266 (*current_liboctave_error_handler) (
"Error creating fftw plan");
277 const double *in,
Complex *out)
279 fftw_plan *cur_plan_p =
reinterpret_cast<fftw_plan *
> (&
rplan);
280 bool create_new_plan =
false;
287 if (
rplan ==
nullptr ||
rd != dist ||
rs != stride ||
rr != rank
288 ||
rh != howmany || ((ioalign !=
rsimd_align) ? ! ioalign :
false))
289 create_new_plan =
true;
294 for (
int i = 0;
i < rank;
i++)
297 create_new_plan =
true;
315 for (
int i = 0, j = rank-1;
i < rank;
i++, j--)
322 bool plan_destroys_in =
true;
328 plan_flags |= FFTW_ESTIMATE;
329 plan_destroys_in =
false;
332 plan_flags |= FFTW_MEASURE;
335 plan_flags |= FFTW_PATIENT;
338 plan_flags |= FFTW_EXHAUSTIVE;
342 plan_flags |= FFTW_MEASURE;
345 plan_flags |= FFTW_ESTIMATE;
346 plan_destroys_in =
false;
352 plan_flags &= ~FFTW_UNALIGNED;
354 plan_flags |= FFTW_UNALIGNED;
357 fftw_destroy_plan (*cur_plan_p);
359 if (plan_destroys_in)
363 itmp =
reinterpret_cast<double *
> 364 (((
reinterpret_cast<ptrdiff_t
>(itmp) + 15) & ~ 0xF) +
365 ((reinterpret_cast<ptrdiff_t> (in)) & 0xF));
368 fftw_plan_many_dft_r2c (rank,
tmp, howmany, itmp,
369 nullptr, stride, dist,
370 reinterpret_cast<fftw_complex *> (out),
371 nullptr, stride, dist, plan_flags);
376 fftw_plan_many_dft_r2c (rank,
tmp, howmany,
377 (const_cast<double *> (in)),
378 nullptr, stride, dist,
379 reinterpret_cast<fftw_complex *> (out),
380 nullptr, stride, dist, plan_flags);
383 if (*cur_plan_p ==
nullptr)
384 (*current_liboctave_error_handler) (
"Error creating fftw plan");
408 fftw_destroy_plan (reinterpret_cast<fftw_plan> (
rplan));
410 fftw_destroy_plan (reinterpret_cast<fftw_plan> (
plan[0]));
412 fftw_destroy_plan (reinterpret_cast<fftw_plan> (
plan[1]));
424 : meth (ESTIMATE), rplan (nullptr), rd (0), rs (0), rr (0), rh (0), rn (),
425 rsimd_align (
false), nthreads (1)
428 d[0] =
d[1] =
s[0] =
s[1] =
r[0] =
r[1] =
h[0] =
h[1] = 0;
433 #if defined (HAVE_FFTW3F_THREADS) 434 int init_ret = fftwf_init_threads ();
436 (*current_liboctave_error_handler) (
"Error initializing FFTW3F threads");
441 fftwf_plan_with_nthreads (
nthreads);
445 fftwf_import_system_wisdom ();
452 plan_p =
reinterpret_cast<fftwf_plan *
> (&
rplan);
454 fftwf_destroy_plan (*plan_p);
456 plan_p =
reinterpret_cast<fftwf_plan *
> (&
plan[0]);
458 fftwf_destroy_plan (*plan_p);
460 plan_p =
reinterpret_cast<fftwf_plan *
> (&
plan[1]);
462 fftwf_destroy_plan (*plan_p);
479 (*current_liboctave_error_handler)
480 (
"unable to create fftw_planner object!");
488 #if defined (HAVE_FFTW3F_THREADS) 492 fftwf_plan_with_nthreads (nt);
497 (*current_liboctave_warning_handler)
498 (
"unable to change number of threads without FFTW thread support");
511 int which = (dir == FFTW_FORWARD) ? 0 : 1;
512 fftwf_plan *cur_plan_p =
reinterpret_cast<fftwf_plan *
> (&
plan[which]);
513 bool create_new_plan =
false;
515 bool ioinplace = (in == out);
521 if (
plan[which] ==
nullptr ||
d[which] != dist ||
s[which] != stride
522 ||
r[which] != rank ||
h[which] != howmany
524 || ((ioalign !=
simd_align[which]) ? ! ioalign :
false))
525 create_new_plan =
true;
530 for (
int i = 0;
i < rank;
i++)
533 create_new_plan =
true;
552 for (
int i = 0, j = rank-1;
i < rank;
i++, j--)
559 bool plan_destroys_in =
true;
565 plan_flags |= FFTW_ESTIMATE;
566 plan_destroys_in =
false;
569 plan_flags |= FFTW_MEASURE;
572 plan_flags |= FFTW_PATIENT;
575 plan_flags |= FFTW_EXHAUSTIVE;
579 plan_flags |= FFTW_MEASURE;
582 plan_flags |= FFTW_ESTIMATE;
583 plan_destroys_in =
false;
589 plan_flags &= ~FFTW_UNALIGNED;
591 plan_flags |= FFTW_UNALIGNED;
594 fftwf_destroy_plan (*cur_plan_p);
596 if (plan_destroys_in)
601 (((
reinterpret_cast<ptrdiff_t
>(itmp) + 15) & ~ 0xF) +
602 ((reinterpret_cast<ptrdiff_t> (in)) & 0xF));
605 fftwf_plan_many_dft (rank,
tmp, howmany,
606 reinterpret_cast<fftwf_complex *> (itmp),
607 nullptr, stride, dist,
608 reinterpret_cast<fftwf_complex *> (out),
609 nullptr, stride, dist, dir, plan_flags);
614 fftwf_plan_many_dft (rank,
tmp, howmany,
615 reinterpret_cast<fftwf_complex *> (const_cast<FloatComplex *> (in)),
616 nullptr, stride, dist,
617 reinterpret_cast<fftwf_complex *> (out),
618 nullptr, stride, dist, dir, plan_flags);
621 if (*cur_plan_p ==
nullptr)
622 (*current_liboctave_error_handler) (
"Error creating fftw plan");
635 fftwf_plan *cur_plan_p =
reinterpret_cast<fftwf_plan *
> (&
rplan);
636 bool create_new_plan =
false;
643 if (
rplan ==
nullptr ||
rd != dist ||
rs != stride ||
rr != rank
644 ||
rh != howmany || ((ioalign !=
rsimd_align) ? ! ioalign :
false))
645 create_new_plan =
true;
650 for (
int i = 0;
i < rank;
i++)
653 create_new_plan =
true;
671 for (
int i = 0, j = rank-1;
i < rank;
i++, j--)
678 bool plan_destroys_in =
true;
684 plan_flags |= FFTW_ESTIMATE;
685 plan_destroys_in =
false;
688 plan_flags |= FFTW_MEASURE;
691 plan_flags |= FFTW_PATIENT;
694 plan_flags |= FFTW_EXHAUSTIVE;
698 plan_flags |= FFTW_MEASURE;
701 plan_flags |= FFTW_ESTIMATE;
702 plan_destroys_in =
false;
708 plan_flags &= ~FFTW_UNALIGNED;
710 plan_flags |= FFTW_UNALIGNED;
713 fftwf_destroy_plan (*cur_plan_p);
715 if (plan_destroys_in)
719 itmp =
reinterpret_cast<float *
> 720 (((
reinterpret_cast<ptrdiff_t
>(itmp) + 15) & ~ 0xF) +
721 ((reinterpret_cast<ptrdiff_t> (in)) & 0xF));
724 fftwf_plan_many_dft_r2c (rank,
tmp, howmany, itmp,
725 nullptr, stride, dist,
726 reinterpret_cast<fftwf_complex *> (out),
727 nullptr, stride, dist, plan_flags);
732 fftwf_plan_many_dft_r2c (rank,
tmp, howmany,
733 (const_cast<float *> (in)),
734 nullptr, stride, dist,
735 reinterpret_cast<fftwf_complex *> (out),
736 nullptr, stride, dist, plan_flags);
739 if (*cur_plan_p ==
nullptr)
740 (*current_liboctave_error_handler) (
"Error creating fftw plan");
764 fftwf_destroy_plan (reinterpret_cast<fftwf_plan> (
rplan));
766 fftwf_destroy_plan (reinterpret_cast<fftwf_plan> (
plan[0]));
768 fftwf_destroy_plan (reinterpret_cast<fftwf_plan> (
plan[1]));
777 template <
typename T>
786 for (
size_t i = 0;
i < nr;
i++)
787 for (
size_t j = nc/2+1; j < nc; j++)
788 out[j*stride +
i*dist] =
conj (out[(nc - j)*stride +
i*dist]);
793 template <
typename T>
800 size_t nrp = nr * np;
807 for (
size_t i = 0;
i < nrp;
i++)
809 ptr1 = out +
i * (nc/2 + 1) + nrp*((nc-1)/2);
811 for (
size_t j = 0; j < nc/2+1; j++)
819 for (
size_t i = 0;
i < np;
i++)
821 for (
size_t j = 1; j < nr; j++)
822 for (
size_t k = nc/2+1;
k < nc;
k++)
823 out[
k + (j +
i*nr)*nc] =
conj (out[nc -
k + ((
i+1)*nr - j)*nc]);
825 for (
size_t j = nc/2+1; j < nc; j++)
826 out[j +
i*nr*nc] =
conj (out[(
i*nr+1)*nc - j]);
833 size_t jstart =
dv(0) *
dv(1);
834 size_t kstep =
dv(0);
837 for (
int inner = 2; inner <
dv.
ndims (); inner++)
839 size_t jmax = jstart *
dv(inner);
840 for (
size_t i = 0;
i < nel;
i+=jmax)
841 for (
size_t j = jstart, jj = jmax-jstart; j < jj;
842 j+=jstart, jj-=jstart)
843 for (
size_t k = 0;
k < jstart;
k+= kstep)
844 for (
size_t l = nc/2+1; l < nc; l++)
846 T
tmp = out[
i+ j +
k + l];
847 out[
i + j +
k + l] = out[
i + jj +
k + l];
848 out[
i + jj +
k + l] =
tmp;
860 dist = (dist < 0 ? npts : dist);
864 stride, dist, in, out);
865 fftw_plan plan =
reinterpret_cast<fftw_plan
> (vplan);
867 fftw_execute_dft_r2c (plan, (const_cast<double *>(in)),
868 reinterpret_cast<fftw_complex *> (out));
881 dist = (dist < 0 ? npts : dist);
887 fftw_plan plan =
reinterpret_cast<fftw_plan
> (vplan);
889 fftw_execute_dft (plan,
890 reinterpret_cast<fftw_complex *> (const_cast<Complex *>(in)),
891 reinterpret_cast<fftw_complex *> (out));
901 dist = (dist < 0 ? npts : dist);
905 stride, dist, in, out);
906 fftw_plan plan =
reinterpret_cast<fftw_plan
> (vplan);
908 fftw_execute_dft (plan,
909 reinterpret_cast<fftw_complex *> (const_cast<Complex *>(in)),
910 reinterpret_cast<fftw_complex *> (out));
913 for (
size_t j = 0; j < nsamples; j++)
914 for (
size_t i = 0;
i < npts;
i++)
915 out[
i*stride + j*dist] /=
scale;
925 for (
int i = 0;
i < rank;
i++)
935 fftw_plan plan =
reinterpret_cast<fftw_plan
> (vplan);
937 fftw_execute_dft_r2c (plan, (const_cast<double *>(in)),
938 reinterpret_cast<fftw_complex *> (out+ offset));
952 for (
int i = 0;
i < rank;
i++)
956 1, 1, dist, in, out);
957 fftw_plan plan =
reinterpret_cast<fftw_plan
> (vplan);
959 fftw_execute_dft (plan,
960 reinterpret_cast<fftw_complex *> (const_cast<Complex *>(in)),
961 reinterpret_cast<fftw_complex *> (out));
971 for (
int i = 0;
i < rank;
i++)
975 1, 1, dist, in, out);
976 fftw_plan plan =
reinterpret_cast<fftw_plan
> (vplan);
978 fftw_execute_dft (plan,
979 reinterpret_cast<fftw_complex *> (const_cast<Complex *>(in)),
980 reinterpret_cast<fftw_complex *> (out));
982 const size_t npts =
dv.
numel ();
984 for (
size_t i = 0;
i < npts;
i++)
994 dist = (dist < 0 ? npts : dist);
999 fftwf_plan plan =
reinterpret_cast<fftwf_plan
> (vplan);
1001 fftwf_execute_dft_r2c (plan, (const_cast<float *>(in)),
1002 reinterpret_cast<fftwf_complex *> (out));
1015 dist = (dist < 0 ? npts : dist);
1019 nsamples, stride, dist,
1021 fftwf_plan plan =
reinterpret_cast<fftwf_plan
> (vplan);
1023 fftwf_execute_dft (plan,
1024 reinterpret_cast<fftwf_complex *> (const_cast<FloatComplex *>(in)),
1025 reinterpret_cast<fftwf_complex *> (out));
1035 dist = (dist < 0 ? npts : dist);
1039 nsamples, stride, dist,
1041 fftwf_plan plan =
reinterpret_cast<fftwf_plan
> (vplan);
1043 fftwf_execute_dft (plan,
1044 reinterpret_cast<fftwf_complex *> (const_cast<FloatComplex *>(in)),
1045 reinterpret_cast<fftwf_complex *> (out));
1048 for (
size_t j = 0; j < nsamples; j++)
1049 for (
size_t i = 0;
i < npts;
i++)
1050 out[
i*stride + j*dist] /=
scale;
1060 for (
int i = 0;
i < rank;
i++)
1070 fftwf_plan plan =
reinterpret_cast<fftwf_plan
> (vplan);
1072 fftwf_execute_dft_r2c (plan, (const_cast<float *>(in)),
1073 reinterpret_cast<fftwf_complex *> (out+ offset));
1087 for (
int i = 0;
i < rank;
i++)
1091 1, 1, dist, in, out);
1092 fftwf_plan plan =
reinterpret_cast<fftwf_plan
> (vplan);
1094 fftwf_execute_dft (plan,
1095 reinterpret_cast<fftwf_complex *> (const_cast<FloatComplex *>(in)),
1096 reinterpret_cast<fftwf_complex *> (out));
1106 for (
int i = 0;
i < rank;
i++)
1110 1, 1, dist, in, out);
1111 fftwf_plan plan =
reinterpret_cast<fftwf_plan
> (vplan);
1113 fftwf_execute_dft (plan,
1114 reinterpret_cast<fftwf_complex *> (const_cast<FloatComplex *>(in)),
1115 reinterpret_cast<fftwf_complex *> (out));
1117 const size_t npts =
dv.
numel ();
1119 for (
size_t i = 0;
i < npts;
i++)
1130 #if defined (HAVE_FFTW) 1140 #if defined (HAVE_FFTW)
FftwMethod do_method(void)
static void convert_packcomplex_Nd(T *out, const dim_vector &dv)
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 void convert_packcomplex_1d(T *out, size_t nr, size_t nc, octave_idx_type stride, octave_idx_type dist)
static bool instance_ok(void)
static int ifft(const Complex *in, Complex *out, size_t npts, 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 &)
~float_fftw_planner(void)
void * do_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)
ComplexColumnVector conj(const ComplexColumnVector &a)
static int fft(const double *in, Complex *out, size_t npts, size_t nsamples=1, octave_idx_type stride=1, octave_idx_type dist=-1)
static void cleanup_instance(void)
std::string fftwf_version(void)
unsigned long int octave_num_processors_wrapper(enum octave_nproc_query octave_query)
std::string fftw_version(void)
static fftw_planner * instance
static float_fftw_planner * instance
void * do_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)
the exceeded dimensions are set to if fewer subscripts than dimensions are the exceeding dimensions are merged into the final requested dimension For consider the following dims
static bool instance_ok(void)
octave_idx_type numel(int n=0) const
Number of elements that a matrix with this dimensions would have.
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)
FftwMethod do_method(void)
static int ifftNd(const Complex *, Complex *, const int, const dim_vector &)
void scale(Matrix &m, double x, double y, double z)
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
#define CHECK_SIMD_ALIGNMENT(x)
std::complex< float > FloatComplex
static void cleanup_instance(void)
octave_idx_type ndims(void) const
Number of dimensions.
std::complex< double > Complex
Vector representing the dimensions (size) of an Array.
If this string is the system will ring the terminal sometimes it is useful to be able to print the original representation of the string