forward and reverse native complex<->complex

Changes

@@ -1,3 +1,5 @@
+- support PDL 2.027+ "native complex" inputs for fftn and ifftn
+
 0.10 2021-03-21
 - fix metadata gremlin
 

FFTW3.pd

@@ -61,7 +61,11 @@ for my $rank (1..$maxrank)
 sub fft$rank { __fft_internal( "fft$rank",\@_ ); }
 *PDL::fft$rank = \\&fft$rank;
 
-sub ifft$rank { my \$a = __fft_internal( "ifft$rank", \@_ ); \$a /= $shapestr; \$a; }
+sub ifft$rank {
+  my \$a = __fft_internal( "ifft$rank", \@_ );
+  \$a /= \$_[0]->type->real ? $shapestr : $rshapestr;
+  \$a;
+}
 *PDL::ifft$rank = \\&ifft$rank;
 
 sub rfft$rank { __fft_internal( "rfft$rank", \@_ ); }
@@ -152,17 +156,14 @@ sub generateDefinitions
   ####### now I generate the definitions for the real-complex and complex-real cases
   my @dims_real    = @dims;
   my @dims_complex = @dims;
-
   shift @dims_real; # get rid of the (real,imag) dimension for the real numbers
   $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');
   $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;
 
@@ -195,6 +196,18 @@ EOF
   $pp_def{Pars} = "complexv($dims_complex_string); [o]real($dims_real_string);";
   $pp_def{Code} = $code_real_backward;
   pp_def( "__irfft$rank", %pp_def );
+
+  ################################################################################
+  ####### now native-complex version
+  shift @dims; # drop the (real,imag) dim
+  $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} =~ s/TFD/TGC/g;
+  $pp_def{Code} =~ s/\*2//g; # the sizeof-doubling
+  $pp_def{GenericTypes} = [qw(G C)];
+  pp_def( "__Nfft$rank", %pp_def );
 }
 
 sub slurp

FFTW3_header_include.pm

@@ -73,9 +73,11 @@ EOF
   my $plan = getPlan( $thisfunction, $rank, $is_real_fft, $do_inverse_fft, $in, $out );
   barf "$thisfunction couldn't make a plan. Giving up\n" unless defined $plan;
 
+  my $is_native = !$in->type->real; # native complex
   # I now have the arguments and the plan. Go!
   my $internal_function = 'PDL::__';
   $internal_function .=
+    $is_native ? 'N' :
     !$is_real_fft ? '' :
     $do_inverse_fft ? 'ir' :
     'r';
@@ -165,14 +167,16 @@ EOF
 sub validateArgumentDimensions_complex
 {
   my ( $rank, $thisfunction, $arg ) = @_;
+  my $is_native = !$arg->type->real;
 
   # complex FFT. Identically-sized inputs/outputs
-  barf <<EOF if $arg->dim(0) != 2;
-$thisfunction must have dim(0) == 2 for the inputs and outputs.
+  barf <<EOF if !$is_native and $arg->dim(0) != 2;
+$thisfunction must have dim(0) == 2 for non-native complex inputs and outputs.
 This is the (real,imag) dimension. Giving up.
 EOF
 
-  barf <<EOF if $arg->ndims-1 < $rank;
+  my $dims_cmp = $arg->ndims - ($is_native ? 0 : 1);
+  barf <<EOF if $dims_cmp < $rank;
 Tried to compute a $rank-dimensional FFT, but an array has fewer than $rank dimensions.
 Giving up.
 EOF
@@ -230,7 +234,6 @@ EOF
 
 sub matchDimensions_real {
   my ($thisfunction, $rank, $do_inverse_fft, $in, $out) = @_;
-
   my ($varname1, $varname2, $var1, $var2);
   if ( !$do_inverse_fft ) {
     # Forward FFT. The input is real, the output is complex. $output->dim(0)
@@ -266,28 +269,19 @@ sub processTypes
   # Input and output types must match, and I can only really deal with float and
   # double. If given an output, I refuse to tweak the type of the output,
   # otherwise, I upgrade to float and then to double
-  if( $$out->isnull )
-  {
-    if( $$in->type < float )
-    {
+  if( $$out->isnull ) {
+    if( $$in->type < float ) {
       forceType( $in, (float) );
     }
-  }
-  else
-  {
+  } else {
     # I'm given an output. Make sure this is of a type I can work with,
     # otherwise give up
-    my $targetType;
-
     my $out_type = $$out->type;
-
     barf <<EOF if $out_type < float;
 $thisfunction can only generate 'float' or 'double' output. You gave an output
 of type '$out_type'. I can't change this so I give up
 EOF
-
-    $targetType = ( $out_type < float ) ? (float) : $out_type;
-
+    my $targetType = ( $out_type < float ) ? (float) : $out_type;
     forceType( $in,  $targetType );
     forceType( $out, $targetType );
   }
@@ -309,9 +303,9 @@ sub getPlan
   my @dims; # the dimensionality of the FFT
   if( !$is_real_fft )
   {
-    # complex FFT - ignore first dimension which is (real, imag)
+    # complex FFT
     @dims = $in->dims;
-    shift @dims;
+    shift @dims if $in->type->real; # ignore first dimension which is (real, imag)
   }
   elsif( !$do_inverse_fft )
   {

Makefile.PL

@@ -59,9 +59,9 @@ $descriptor{depend} = { 'FFTW3.pm' => join(' ', qw(template_complex.c template_r
                                                    compute_plan_template.xs
                                                    FFTW3_header_include.pm)) };
 
-# This Alien::FFTW3 requirement is actually optional. pkg-config is used if this
-# isn't available
-$descriptor{PREREQ_PM} = {'PDL::Core'=>2.001};
+$descriptor{PREREQ_PM} = {
+  'PDL' => '2.030', # working $type->real method
+};
 $descriptor{CONFIGURE_REQUIRES} = {
   'PDL::Core::Dev' =>0,
   'IPC::Run'       =>0,

README.pod

@@ -146,6 +146,10 @@ output will still currently be a similarly-shaped "real" L<PDL>
 object with the initial dimension of 2. This is intended to be changed
 so the output type is the same as the input.
 
+As of version 0.11, you can also pass in piddles with the new "native
+complex" types (C<cfloat>, C<cdouble>), without the initial dimension of
+2. Outputs will also be native complex.
+
 Generally, the sizes of the input and the output must match. This is completely
 true for the complex <-> complex transforms: the output will have the same size
 and the input, and an error will result if this isn't possible for some reason.
@@ -207,8 +211,10 @@ identical. The output parameter is optional and, if present, must be
 the last argument. If the output piddle is passed in, the user I<must>
 make sure the dimensions match.
 
-The 0 dim of the input PDL must have size 2 and run over (real,imaginary)
-components. The transform is carried out over dims 1 through X.
+If PDL 2.027+ "native complex" data is the input, the dimensions are as
+you'd expect. Otherwise, the 0 dim of the input PDL must have size 2 and
+run over (real,imaginary) components. The transform is carried out over
+the remaining dims.
 
 The fftn form takes a minimum of two arguments: the PDL to transform,
 and the number of dimensions to transform as a separate argument.

t/fftw.t

@@ -31,6 +31,14 @@ use constant approx_eps_single => 1e-3;
 
 my $Nplans = 0;
 
+sub other2native {
+  my ($other) = @_;
+  my @other_dims = $other->dims;
+  shift @other_dims; # drop initial 2
+  $other = cplx $other if !UNIVERSAL::isa($other, 'PDL::Complex');
+  zeroes(cdouble, @other_dims) + $other->re + $other->im * ci();
+}
+
 # 1D basic test
 {
   my $x = sequence(10)->cat(sequence(10)**2)->mv(-1,0);
@@ -61,8 +69,15 @@ my $Nplans = 0;
   ok_should_make_plan( all( approx( $cplx_result, $Xref, approx_eps_double) ),
                       "Basic 1D complex FFT with PDL::Complex" );
 
-  ok_should_make_plan( all( approx( ifft1(fft1($x)), $x , approx_eps_double) ),
+  ok_should_make_plan( all( approx( ifft1(fft1($x)), $x, approx_eps_double) ),
                       "Basic 1D complex FFT - inverse(forward) should be the same (normalized)" );
+
+  my $x_nat = other2native($x_cplx);
+  ok_should_make_plan( all( approx( fft1($x_nat), other2native($Xref), approx_eps_double) ),
+                      "Basic 1D native complex FFT" );
+
+  ok_should_make_plan( all( approx( ifft1(fft1($x_nat)), $x_nat, approx_eps_double) ),
+                      "Basic 1D native complex FFT - inverse" );
 }
 
 # 2D basic test
@@ -83,6 +98,10 @@ my $Nplans = 0;
   ok_should_make_plan( all( approx( fft2($x), $Xref, approx_eps_double) ),
      "Basic 2D complex FFT - double precision" );
 
+  my $x_nat = other2native($x);
+  ok_should_make_plan( all( approx( fft2($x_nat), other2native($Xref), approx_eps_double) ),
+     "Basic 2D native complex FFT - double precision" );
+
   ok_should_make_plan( all( approx( fft2(float $x), float($Xref), approx_eps_single) ),
      "Basic 2D complex FFT - single precision" );
 
@@ -145,6 +164,11 @@ my $Nplans = 0;
 
   ok_should_reuse_plan( all( approx( $f, $Xref, approx_eps_double) ),
      "1D FFTs threaded inside a 3D piddle" );
+
+  my $x_nat = other2native($x);
+  my $f_nat = fft1($x_nat);
+  ok_should_reuse_plan( all( approx( $f_nat, other2native($Xref), approx_eps_double) ),
+     "1D native complex FFTs threaded inside a 3D piddle" );
 }
 
 # try out some different ways of calling the module, make sure the argument