diff --git a/Makefile.PL b/Makefile.PL index 54d9d9d..e7fb714 100644 --- a/Makefile.PL +++ b/Makefile.PL @@ -67,6 +67,7 @@ $descriptor{BUILD_REQUIRES} = {'PDL::PP'=>0}; $descriptor{AUTHOR} = "Dima Kogan , Craig DeForest "; $descriptor{ABSTRACT} = "PDL interface to the Fastest Fourier Transform in the West"; $descriptor{LICENSE} = "gpl_3"; +$descriptor{MIN_PERL_VERSION} = "5.016"; $descriptor{META_MERGE} = { "meta-spec" => { version => 2 }, diff --git a/fftw3.pd b/fftw3.pd index 1c673d5..b986b7c 100644 --- a/fftw3.pd +++ b/fftw3.pd @@ -36,61 +36,39 @@ pp_addhdr( ' #define static_assert_fftw(x) (void)( sizeof( int[ 1 - 2* !(x) ]) ) ' ); - # I want to be able to say $X = fft1($x); rank is required. 'fft()' is ambiguous # about whether threading is desired or if a large fft is desired. Old PDL::FFTW # did one thing, matlab does another, so I do not include this function at all -my $TEMPLATE_REAL = <<'EOF'; -// This is the template used by PP to generate the FFTW routines. - -// This is passed into pp_def() in the 'Code' key. Before this file is passed to -// pp_def, the following strings are replaced: -// -// INVERSE If this is a c2r transform rather than r2c -// RANK The rank of this transform - -{ - // make sure the PDL data type I'm using matches the FFTW data type - static_assert_fftw( sizeof($GENERIC())*2 == sizeof($TFD(fftwf_,fftw_)complex) ); - - $TFD(fftwf_,fftw_)plan plan = INT2PTR( $TFD(fftwf_,fftw_)plan, $COMP(plan)); - - if( !INVERSE ) - $TFD(fftwf_,fftw_)execute_dft_r2c( plan, - ($TFD(float,double)*)$P(real), - ($TFD(fftwf_,fftw_)complex*)$P(complexv) ); - else - { - // FFTW inverse real transforms clobber their input. I thus make a new - // buffer and transform from there - unsigned long nelem = 1; - for( int i=0; i<=RANK; i++ ) - nelem *= $PDL(complexv)->dims[i]; - unsigned long elem_scale = sizeof($GENERIC()) / sizeof( $TFD(float,double) ); /* native complex */ - /* explicit C types here means when changed to $TGC will still be right */ - $TFD(float,double)* input_copy = $TFD(fftwf_,fftw_)alloc_real( nelem * elem_scale ); - memcpy( input_copy, $P(complexv), sizeof($GENERIC()) * nelem ); - - $TFD(fftwf_,fftw_)execute_dft_c2r( plan, - ($TFD(fftwf_,fftw_)complex*)input_copy, - ($TFD(float,double)*)$P(real) ); - - $TFD(fftwf_,fftw_)free( input_copy ); - } -} +my $TEMPLATE_REAL_R2C = <<'EOF'; +// make sure the PDL data type I'm using matches the FFTW data type +static_assert_fftw(sizeof($GENERIC())*2 == sizeof($TFD(fftwf_,fftw_)complex)); +$TFD(fftwf_,fftw_)plan plan = INT2PTR($TFD(fftwf_,fftw_)plan, $COMP(plan)); +$TFD(fftwf_,fftw_)execute_dft_r2c(plan, (void*)$P(real), (void*)$P(complexv)); +EOF +my $TEMPLATE_REAL_C2R = <<'EOF'; +// RANK is replaced with the rank of this transform +// make sure the PDL data type I'm using matches the FFTW data type +static_assert_fftw(sizeof($GENERIC())*2 == sizeof($TFD(fftwf_,fftw_)complex)); +$TFD(fftwf_,fftw_)plan plan = INT2PTR($TFD(fftwf_,fftw_)plan, $COMP(plan)); +// FFTW inverse real transforms clobber their input. I thus make a new +// buffer and transform from there +unsigned long nelem = 1; +for( int i=0; i<=RANK; i++ ) + nelem *= $PDL(complexv)->dims[i]; +unsigned long elem_scale = sizeof($GENERIC()) / sizeof( $TFD(float,double) ); /* native complex */ +/* explicit C types here means when changed to $TGC will still be right */ +$TFD(float,double)* input_copy = $TFD(fftwf_,fftw_)alloc_real( nelem * elem_scale ); +memcpy( input_copy, $P(complexv), sizeof($GENERIC()) * nelem ); +$TFD(fftwf_,fftw_)execute_dft_c2r( plan, (void*)input_copy, (void*)$P(real) ); +$TFD(fftwf_,fftw_)free( input_copy ); EOF my $TEMPLATE_COMPLEX = <<'EOF'; // This is the template used by PP to generate the FFTW routines. -{ - // make sure the PDL data type I'm using matches the FFTW data type - static_assert_fftw( sizeof($GENERIC())*2 == sizeof($TFD(fftwf_,fftw_)complex) ); - - $TFD(fftwf_,fftw_)plan plan = INT2PTR( $TFD(fftwf_,fftw_)plan, $COMP(plan)); - $TFD(fftwf_,fftw_)execute_dft( plan, - ($TFD(fftwf_,fftw_)complex*)$P(in), - ($TFD(fftwf_,fftw_)complex*)$P(out) ); -} +// make sure the PDL data type I'm using matches the FFTW data type +static_assert_fftw(sizeof($GENERIC())*2 == sizeof($TFD(fftwf_,fftw_)complex)); +$TFD(fftwf_,fftw_)plan plan = INT2PTR($TFD(fftwf_,fftw_)plan, $COMP(plan)); +$TFD(fftwf_,fftw_)execute_dft( plan, (void*)$P(in), (void*)$P(out) ); EOF # I define up to rank-10 FFTs. This is annoyingly arbitrary, but hopefully @@ -736,12 +714,7 @@ sub generateDefinitions $dims_complex[1] = 'nhalf'; # first complex dim is real->dim(0)/2+1 my $dims_real_string = join(',', @dims_real); my $dims_complex_string = join(',', @dims_complex); - my $code_real = $TEMPLATE_REAL; - $code_real =~ s/RANK/$rank/ge; - my $code_real_forward = $code_real; - my $code_real_backward = $code_real; - $code_real_forward =~ s/INVERSE/0/g; - $code_real_backward =~ s/INVERSE/1/g; + my $code_real_backward = $TEMPLATE_REAL_C2R =~ s/RANK/$rank/gr; # forward # I have the real dimensions, but not nhalf @@ -750,7 +723,7 @@ if( $PDL(complexv)->ndims <= 1 || $PDL(complexv)->dims[1] <= 0 ) $SIZE(nhalf) = (int)( $PDL(real)->dims[0]/2 ) + 1; EOF $pp_def{Pars} = "real($dims_real_string); [o]complexv($dims_complex_string);"; - $pp_def{Code} = $code_real_forward; + $pp_def{Code} = $TEMPLATE_REAL_R2C; pp_def( "__rfft$rank", %pp_def ); ################################################################################ @@ -772,12 +745,7 @@ EOF $dims_complex[0] = 'nhalf'; # first complex dim is real->dim(0)/2+1 $dims_real_string = join(',', @dims_real); $dims_complex_string = join(',', @dims_complex); - $code_real = $TEMPLATE_REAL; - $code_real =~ s/RANK/$rank-1/ge; - $code_real_forward = $code_real; - $code_real_backward = $code_real; - $code_real_forward =~ s/INVERSE/0/g; - $code_real_backward =~ s/INVERSE/1/g; + $code_real_backward = $TEMPLATE_REAL_C2R =~ s/RANK/$rank-1/gre; # forward $pp_def{RedoDimsCode} = <<'EOF'; @@ -785,7 +753,7 @@ if( $PDL(complexv)->ndims <= 1 || $PDL(complexv)->dims[1] <= 0 ) $SIZE(nhalf) = (int)( $PDL(real)->dims[0]/2 ) + 1; EOF $pp_def{Pars} = "real($dims_real_string); complex [o]complexv($dims_complex_string);"; - $pp_def{Code} = $code_real_forward; + $pp_def{Code} = $TEMPLATE_REAL_R2C; $pp_def{GenericTypes} = [qw(F D)]; pp_def( "__rNfft$rank", %pp_def );