23 #if defined (HAVE_CONFIG_H)
31 #if defined (HAVE_FFTW3_H)
41 #if defined (HAVE_FFTW3_THREADS) || defined (HAVE_FFTW3F_THREADS)
45 #if defined (HAVE_FFTW)
68 : meth (ESTIMATE), rplan (0), rd (0), rs (0), rr (0), rh (0), rn (),
69 rsimd_align (
false), nthreads (1)
72 d[0] =
d[1] =
s[0] =
s[1] =
r[0] =
r[1] =
h[0] =
h[1] = 0;
77 #if defined (HAVE_FFTW3_THREADS)
78 int init_ret = fftw_init_threads ();
80 (*current_liboctave_error_handler) (
"Error initializing FFTW threads");
85 fftw_plan_with_nthreads (nthreads);
89 fftw_import_system_wisdom ();
96 plan_p =
reinterpret_cast<fftw_plan *
> (&
rplan);
98 fftw_destroy_plan (*plan_p);
100 plan_p =
reinterpret_cast<fftw_plan *
> (&
plan[0]);
102 fftw_destroy_plan (*plan_p);
104 plan_p =
reinterpret_cast<fftw_plan *
> (&
plan[1]);
106 fftw_destroy_plan (*plan_p);
123 (*current_liboctave_error_handler)
124 (
"unable to create octave_fftw_planner object!");
132 #if defined (HAVE_FFTW3_THREADS)
136 fftw_plan_with_nthreads (nt);
141 (*current_liboctave_warning_handler)
142 (
"unable to change number of threads without FFTW thread support");
146 #define CHECK_SIMD_ALIGNMENT(x) \
147 (((reinterpret_cast<ptrdiff_t> (x)) & 0xF) == 0)
157 int which = (dir == FFTW_FORWARD) ? 0 : 1;
158 fftw_plan *cur_plan_p =
reinterpret_cast<fftw_plan *
> (&
plan[which]);
159 bool create_new_plan =
false;
161 bool ioinplace = (in == out);
167 if (
plan[which] == 0 ||
d[which] != dist ||
s[which] != stride
168 ||
r[which] != rank ||
h[which] != howmany
170 || ((ioalign !=
simd_align[which]) ? ! ioalign :
false))
171 create_new_plan =
true;
176 for (
int i = 0;
i < rank;
i++)
179 create_new_plan =
true;
198 for (
int i = 0, j = rank-1;
i < rank;
i++, j--)
205 bool plan_destroys_in =
true;
211 plan_flags |= FFTW_ESTIMATE;
212 plan_destroys_in =
false;
215 plan_flags |= FFTW_MEASURE;
218 plan_flags |= FFTW_PATIENT;
221 plan_flags |= FFTW_EXHAUSTIVE;
225 plan_flags |= FFTW_MEASURE;
228 plan_flags |= FFTW_ESTIMATE;
229 plan_destroys_in =
false;
235 plan_flags &= ~FFTW_UNALIGNED;
237 plan_flags |= FFTW_UNALIGNED;
240 fftw_destroy_plan (*cur_plan_p);
242 if (plan_destroys_in)
246 itmp =
reinterpret_cast<Complex *
>
247 (((
reinterpret_cast<ptrdiff_t
>(itmp) + 15) & ~ 0xF) +
248 ((reinterpret_cast<ptrdiff_t> (in)) & 0xF));
251 fftw_plan_many_dft (rank,
tmp, howmany,
252 reinterpret_cast<fftw_complex *> (itmp),
254 reinterpret_cast<fftw_complex *> (out),
255 0, stride, dist, dir, plan_flags);
260 fftw_plan_many_dft (rank,
tmp, howmany,
261 reinterpret_cast<fftw_complex *> (const_cast<Complex *> (in)),
263 reinterpret_cast<fftw_complex *> (out),
264 0, stride, dist, dir, plan_flags);
267 if (*cur_plan_p == 0)
268 (*current_liboctave_error_handler) (
"Error creating fftw plan");
279 const double *in,
Complex *out)
281 fftw_plan *cur_plan_p =
reinterpret_cast<fftw_plan *
> (&
rplan);
282 bool create_new_plan =
false;
289 if (
rplan == 0 ||
rd != dist ||
rs != stride ||
rr != rank
290 ||
rh != howmany || ((ioalign !=
rsimd_align) ? ! ioalign :
false))
291 create_new_plan =
true;
296 for (
int i = 0;
i < rank;
i++)
299 create_new_plan =
true;
317 for (
int i = 0, j = rank-1;
i < rank;
i++, j--)
324 bool plan_destroys_in =
true;
330 plan_flags |= FFTW_ESTIMATE;
331 plan_destroys_in =
false;
334 plan_flags |= FFTW_MEASURE;
337 plan_flags |= FFTW_PATIENT;
340 plan_flags |= FFTW_EXHAUSTIVE;
344 plan_flags |= FFTW_MEASURE;
347 plan_flags |= FFTW_ESTIMATE;
348 plan_destroys_in =
false;
354 plan_flags &= ~FFTW_UNALIGNED;
356 plan_flags |= FFTW_UNALIGNED;
359 fftw_destroy_plan (*cur_plan_p);
361 if (plan_destroys_in)
365 itmp =
reinterpret_cast<double *
>
366 (((
reinterpret_cast<ptrdiff_t
>(itmp) + 15) & ~ 0xF) +
367 ((reinterpret_cast<ptrdiff_t> (in)) & 0xF));
370 fftw_plan_many_dft_r2c (rank,
tmp, howmany, itmp,
372 reinterpret_cast<fftw_complex *> (out),
373 0, stride, dist, plan_flags);
378 fftw_plan_many_dft_r2c (rank,
tmp, howmany,
379 (const_cast<double *> (in)),
381 reinterpret_cast<fftw_complex *> (out),
382 0, stride, dist, plan_flags);
385 if (*cur_plan_p == 0)
386 (*current_liboctave_error_handler) (
"Error creating fftw plan");
410 fftw_destroy_plan (reinterpret_cast<fftw_plan> (
rplan));
412 fftw_destroy_plan (reinterpret_cast<fftw_plan> (
plan[0]));
414 fftw_destroy_plan (reinterpret_cast<fftw_plan> (
plan[1]));
426 : meth (ESTIMATE), rplan (0), rd (0), rs (0), rr (0), rh (0), rn (),
427 rsimd_align (
false), nthreads (1)
430 d[0] =
d[1] =
s[0] =
s[1] =
r[0] =
r[1] =
h[0] =
h[1] = 0;
435 #if defined (HAVE_FFTW3F_THREADS)
436 int init_ret = fftwf_init_threads ();
438 (*current_liboctave_error_handler) (
"Error initializing FFTW3F threads");
443 fftwf_plan_with_nthreads (nthreads);
447 fftwf_import_system_wisdom ();
454 plan_p =
reinterpret_cast<fftwf_plan *
> (&
rplan);
456 fftwf_destroy_plan (*plan_p);
458 plan_p =
reinterpret_cast<fftwf_plan *
> (&
plan[0]);
460 fftwf_destroy_plan (*plan_p);
462 plan_p =
reinterpret_cast<fftwf_plan *
> (&
plan[1]);
464 fftwf_destroy_plan (*plan_p);
481 (*current_liboctave_error_handler)
482 (
"unable to create octave_fftw_planner object!");
490 #if defined (HAVE_FFTW3F_THREADS)
494 fftwf_plan_with_nthreads (nt);
499 (*current_liboctave_warning_handler)
500 (
"unable to change number of threads without FFTW thread support");
513 int which = (dir == FFTW_FORWARD) ? 0 : 1;
514 fftwf_plan *cur_plan_p =
reinterpret_cast<fftwf_plan *
> (&
plan[which]);
515 bool create_new_plan =
false;
517 bool ioinplace = (in == out);
523 if (
plan[which] == 0 ||
d[which] != dist ||
s[which] != stride
524 ||
r[which] != rank ||
h[which] != howmany
526 || ((ioalign !=
simd_align[which]) ? ! ioalign :
false))
527 create_new_plan =
true;
532 for (
int i = 0;
i < rank;
i++)
535 create_new_plan =
true;
554 for (
int i = 0, j = rank-1;
i < rank;
i++, j--)
561 bool plan_destroys_in =
true;
567 plan_flags |= FFTW_ESTIMATE;
568 plan_destroys_in =
false;
571 plan_flags |= FFTW_MEASURE;
574 plan_flags |= FFTW_PATIENT;
577 plan_flags |= FFTW_EXHAUSTIVE;
581 plan_flags |= FFTW_MEASURE;
584 plan_flags |= FFTW_ESTIMATE;
585 plan_destroys_in =
false;
591 plan_flags &= ~FFTW_UNALIGNED;
593 plan_flags |= FFTW_UNALIGNED;
596 fftwf_destroy_plan (*cur_plan_p);
598 if (plan_destroys_in)
603 (((
reinterpret_cast<ptrdiff_t
>(itmp) + 15) & ~ 0xF) +
604 ((reinterpret_cast<ptrdiff_t> (in)) & 0xF));
607 fftwf_plan_many_dft (rank,
tmp, howmany,
608 reinterpret_cast<fftwf_complex *> (itmp),
610 reinterpret_cast<fftwf_complex *> (out),
611 0, stride, dist, dir, plan_flags);
616 fftwf_plan_many_dft (rank,
tmp, howmany,
617 reinterpret_cast<fftwf_complex *> (const_cast<FloatComplex *> (in)),
619 reinterpret_cast<fftwf_complex *> (out),
620 0, stride, dist, dir, plan_flags);
623 if (*cur_plan_p == 0)
624 (*current_liboctave_error_handler) (
"Error creating fftw plan");
638 fftwf_plan *cur_plan_p =
reinterpret_cast<fftwf_plan *
> (&
rplan);
639 bool create_new_plan =
false;
646 if (
rplan == 0 ||
rd != dist ||
rs != stride ||
rr != rank
647 ||
rh != howmany || ((ioalign !=
rsimd_align) ? ! ioalign :
false))
648 create_new_plan =
true;
653 for (
int i = 0;
i < rank;
i++)
656 create_new_plan =
true;
674 for (
int i = 0, j = rank-1;
i < rank;
i++, j--)
681 bool plan_destroys_in =
true;
687 plan_flags |= FFTW_ESTIMATE;
688 plan_destroys_in =
false;
691 plan_flags |= FFTW_MEASURE;
694 plan_flags |= FFTW_PATIENT;
697 plan_flags |= FFTW_EXHAUSTIVE;
701 plan_flags |= FFTW_MEASURE;
704 plan_flags |= FFTW_ESTIMATE;
705 plan_destroys_in =
false;
711 plan_flags &= ~FFTW_UNALIGNED;
713 plan_flags |= FFTW_UNALIGNED;
716 fftwf_destroy_plan (*cur_plan_p);
718 if (plan_destroys_in)
722 itmp =
reinterpret_cast<float *
>
723 (((
reinterpret_cast<ptrdiff_t
>(itmp) + 15) & ~ 0xF) +
724 ((reinterpret_cast<ptrdiff_t> (in)) & 0xF));
727 fftwf_plan_many_dft_r2c (rank,
tmp, howmany, itmp,
729 reinterpret_cast<fftwf_complex *> (out),
730 0, stride, dist, plan_flags);
735 fftwf_plan_many_dft_r2c (rank,
tmp, howmany,
736 (const_cast<float *> (in)),
738 reinterpret_cast<fftwf_complex *> (out),
739 0, stride, dist, plan_flags);
742 if (*cur_plan_p == 0)
743 (*current_liboctave_error_handler) (
"Error creating fftw plan");
767 fftwf_destroy_plan (reinterpret_cast<fftwf_plan> (
rplan));
769 fftwf_destroy_plan (reinterpret_cast<fftwf_plan> (
plan[0]));
771 fftwf_destroy_plan (reinterpret_cast<fftwf_plan> (
plan[1]));
780 template <
typename T>
789 for (
size_t i = 0;
i < nr;
i++)
790 for (
size_t j = nc/2+1; j < nc; j++)
791 out[j*stride +
i*dist] =
conj (out[(nc - j)*stride +
i*dist]);
796 template <
typename T>
802 size_t np = (dv.
ndims () > 2 ? dv.
numel () / nc / nr : 1);
803 size_t nrp = nr * np;
810 for (
size_t i = 0;
i < nrp;
i++)
812 ptr1 = out +
i * (nc/2 + 1) + nrp*((nc-1)/2);
814 for (
size_t j = 0; j < nc/2+1; j++)
822 for (
size_t i = 0;
i < np;
i++)
824 for (
size_t j = 1; j < nr; j++)
825 for (
size_t k = nc/2+1;
k < nc;
k++)
826 out[
k + (j +
i*nr)*nc] =
conj (out[nc -
k + ((
i+1)*nr - j)*nc]);
828 for (
size_t j = nc/2+1; j < nc; j++)
829 out[j +
i*nr*nc] =
conj (out[(
i*nr+1)*nc - j]);
836 size_t jstart =
dv(0) *
dv(1);
837 size_t kstep =
dv(0);
838 size_t nel = dv.
numel ();
840 for (
int inner = 2; inner < dv.
ndims (); inner++)
842 size_t jmax = jstart *
dv(inner);
843 for (
size_t i = 0;
i < nel;
i+=jmax)
844 for (
size_t j = jstart, jj = jmax-jstart; j < jj;
845 j+=jstart, jj-=jstart)
846 for (
size_t k = 0;
k < jstart;
k+= kstep)
847 for (
size_t l = nc/2+1; l < nc; l++)
849 T
tmp = out[
i+ j +
k + l];
850 out[
i + j +
k + l] = out[
i + jj +
k + l];
851 out[
i + jj +
k + l] =
tmp;
863 dist = (dist < 0 ? npts : dist);
867 stride, dist, in, out);
868 fftw_plan plan =
reinterpret_cast<fftw_plan
> (vplan);
870 fftw_execute_dft_r2c (plan, (const_cast<double *>(in)),
871 reinterpret_cast<fftw_complex *> (out));
884 dist = (dist < 0 ? npts : dist);
890 fftw_plan plan =
reinterpret_cast<fftw_plan
> (vplan);
892 fftw_execute_dft (plan,
893 reinterpret_cast<fftw_complex *> (const_cast<Complex *>(in)),
894 reinterpret_cast<fftw_complex *> (out));
904 dist = (dist < 0 ? npts : dist);
910 fftw_plan plan =
reinterpret_cast<fftw_plan
> (vplan);
912 fftw_execute_dft (plan,
913 reinterpret_cast<fftw_complex *> (const_cast<Complex *>(in)),
914 reinterpret_cast<fftw_complex *> (out));
917 for (
size_t j = 0; j < nsamples; j++)
918 for (
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 plan =
reinterpret_cast<fftw_plan
> (vplan);
941 fftw_execute_dft_r2c (plan, (const_cast<double *>(in)),
942 reinterpret_cast<fftw_complex *> (out+ offset));
956 for (
int i = 0;
i < rank;
i++)
960 dv, 1, 1, dist, in, out);
961 fftw_plan plan =
reinterpret_cast<fftw_plan
> (vplan);
963 fftw_execute_dft (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 dv, 1, 1, dist, in, out);
980 fftw_plan plan =
reinterpret_cast<fftw_plan
> (vplan);
982 fftw_execute_dft (plan,
983 reinterpret_cast<fftw_complex *> (const_cast<Complex *>(in)),
984 reinterpret_cast<fftw_complex *> (out));
986 const size_t npts = dv.
numel ();
988 for (
size_t i = 0;
i < npts;
i++)
998 dist = (dist < 0 ? npts : dist);
1004 fftwf_plan plan =
reinterpret_cast<fftwf_plan
> (vplan);
1006 fftwf_execute_dft_r2c (plan, (const_cast<float *>(in)),
1007 reinterpret_cast<fftwf_complex *> (out));
1020 dist = (dist < 0 ? npts : dist);
1027 fftwf_plan plan =
reinterpret_cast<fftwf_plan
> (vplan);
1029 fftwf_execute_dft (plan,
1030 reinterpret_cast<fftwf_complex *> (const_cast<FloatComplex *>(in)),
1031 reinterpret_cast<fftwf_complex *> (out));
1041 dist = (dist < 0 ? npts : dist);
1048 fftwf_plan plan =
reinterpret_cast<fftwf_plan
> (vplan);
1050 fftwf_execute_dft (plan,
1051 reinterpret_cast<fftwf_complex *> (const_cast<FloatComplex *>(in)),
1052 reinterpret_cast<fftwf_complex *> (out));
1055 for (
size_t j = 0; j < nsamples; j++)
1056 for (
size_t i = 0;
i < npts;
i++)
1057 out[
i*stride + j*dist] /= scale;
1067 for (
int i = 0;
i < rank;
i++)
1078 fftwf_plan plan =
reinterpret_cast<fftwf_plan
> (vplan);
1080 fftwf_execute_dft_r2c (plan, (const_cast<float *>(in)),
1081 reinterpret_cast<fftwf_complex *> (out+ offset));
1095 for (
int i = 0;
i < rank;
i++)
1101 fftwf_plan plan =
reinterpret_cast<fftwf_plan
> (vplan);
1103 fftwf_execute_dft (plan,
1104 reinterpret_cast<fftwf_complex *> (const_cast<FloatComplex *>(in)),
1105 reinterpret_cast<fftwf_complex *> (out));
1115 for (
int i = 0;
i < rank;
i++)
1121 fftwf_plan plan =
reinterpret_cast<fftwf_plan
> (vplan);
1123 fftwf_execute_dft (plan,
1124 reinterpret_cast<fftwf_complex *> (const_cast<FloatComplex *>(in)),
1125 reinterpret_cast<fftwf_complex *> (out));
1127 const size_t npts = dv.
numel ();
1129 for (
size_t i = 0;
i < npts;
i++)
1140 #if defined (HAVE_FFTW)
1141 return fftw_version;
1150 #if defined (HAVE_FFTW)
1151 return fftwf_version;
static octave_fftw_planner * instance
static bool instance_ok(void)
~octave_float_fftw_planner(void)
static void convert_packcomplex_1d(T *out, size_t nr, size_t nc, octave_idx_type stride, octave_idx_type dist)
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 octave_idx_type nn
octave_fftw_planner(void)
static octave_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 Complex *in, Complex *out)
octave_float_fftw_planner(void)
octave_idx_type numel(int n=0) const
Number of elements that a matrix with this dimensions would have.
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 Complex *in, Complex *out)
ComplexColumnVector conj(const ComplexColumnVector &a)
FftwMethod do_method(void)
std::string octave_fftwf_version(void)
static void convert_packcomplex_Nd(T *out, const dim_vector &dv)
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 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)
unsigned long int octave_num_processors_wrapper(enum octave_nproc_query octave_query)
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 int ifftNd(const Complex *, Complex *, const int, const dim_vector &)
std::string octave_fftw_version(void)
=val(i)}if ode{val(i)}occurs in table i
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)
static bool instance_ok(void)
octave_idx_type ndims(void) const
Number of dimensions.
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
std::complex< double > Complex
~octave_fftw_planner(void)
FftwMethod do_method(void)
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
static void cleanup_instance(void)
static void cleanup_instance(void)