incorporate template_*.c into .pd

MANIFEST

@@ -6,6 +6,4 @@ MANIFEST			This list of files
 README.pod
 t/fftw.t
 t/threads.t
-template_complex.c
-template_real.c
 TODO

Makefile.PL

@@ -55,8 +55,6 @@ push @{$descriptor{LIBS}  },  $libs;
 $descriptor{INC} = '' unless defined $descriptor{INC};
 $descriptor{INC} .= " $cflags";
 
-$descriptor{depend} = { 'FFTW3.pm' => join(' ', qw(template_complex.c template_real.c)) };
-
 $descriptor{PREREQ_PM} = {
   'PDL' => '2.049', # as_native
 };

fftw3.pd

@@ -31,15 +31,67 @@ EOF
 
 pp_addhdr( '
 #include <fftw3.h>
+
+/* the Linux kernel does something similar to assert at compile time */
+#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 );
+  }
+}
+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) );
+}
+EOF
 
 # I define up to rank-10 FFTs. This is annoyingly arbitrary, but hopefully
 # should be sufficient
@@ -653,7 +705,6 @@ pp_add_exported( map "${_}fftn", '','i','r','ir' );
 
 pp_done();
 
-
 sub generateDefinitions
 {
   my $rank = shift;
@@ -676,7 +727,7 @@ sub generateDefinitions
   unshift @dims, 'n0=2';
   my $dims_string = join(',', @dims);
   $pp_def{Pars} = "in($dims_string); [o]out($dims_string);";
-  $pp_def{Code} = slurp('template_complex.c');
+  $pp_def{Code} = $TEMPLATE_COMPLEX;
   pp_def( "__fft$rank", %pp_def );
 
   ##################################################################################
@@ -687,7 +738,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          =  slurp('template_real.c');
+  my $code_real          =  $TEMPLATE_REAL;
   $code_real             =~ s/RANK/$rank/ge;
   my $code_real_forward  =  $code_real;
   my $code_real_backward =  $code_real;
@@ -730,7 +781,7 @@ EOF
   $dims_string = join(',', @dims);
   delete $pp_def{RedoDimsCode};
   $pp_def{Pars} = "in($dims_string); [o]out($dims_string);";
-  $pp_def{Code} = slurp('template_complex.c');
+  $pp_def{Code} = $TEMPLATE_COMPLEX;
   $pp_def{Code} =~ s/TFD/TGC/g;
   $pp_def{Code} =~ s/\*2//g; # the sizeof-doubling
   $pp_def{GenericTypes} = [qw(G C)];
@@ -743,7 +794,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          =  slurp('template_real.c');
+  $code_real          =  $TEMPLATE_REAL;
   $code_real             =~ s/RANK/$rank-1/ge;
   $code_real_forward  =  $code_real;
   $code_real_backward =  $code_real;
@@ -773,13 +824,9 @@ EOF
   pp_def( "__irNfft$rank", %pp_def );
 }
 
-sub slurp
-{
+sub slurp {
   my $filename = shift;
-  open FD, '<', $filename or die "Couldn't open '$filename' for rading";
-
+  open my $fh, '<', $filename or die "open '$filename' for reading: $!";
   local $/ = undef;
-  my $contents = <FD>;
-  close FD;
-  return qq{\n#line 0 "$filename"\n\n} . $contents;
+  return qq{\n#line 0 "$filename"\n\n} . <$fh>;
 }