Skip to content

Commit

Permalink
split TEMPLATE_REAL into R2C and C2R parts
Browse files Browse the repository at this point in the history
  • Loading branch information
@mohawk2
mohawk2 committed Nov 28, 2024
1 parent b540263 commit 55d746a
Show file tree
Hide file tree
Showing 2 changed files with 31 additions and 62 deletions.
1 change: 1 addition & 0 deletions Makefile.PL
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,6 +67,7 @@ $descriptor{BUILD_REQUIRES} = {'PDL::PP'=>0};
$descriptor{AUTHOR} = "Dima Kogan <dima\@secretsauce.net>, Craig DeForest <deforest\@boulder.swri.edu>";
$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 },
Expand Down
92 changes: 30 additions & 62 deletions fftw3.pd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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 );

################################################################################
Expand All @@ -772,20 +745,15 @@ 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';
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 );

Expand Down

0 comments on commit 55d746a

Please sign in to comment.