r2c/c2r FFTW advanced interface question/example

I'd like to use FFTW's "advanced" interface to compute the Fourier
transform of the "fast" and "slow" axis of a 2D array individually.
The input data are real.  An in-place transform is required.

I will attach a simple Fortran 90 test program which illustrates a
problem that I'm having.

The program initializes a FFTW "padded" real array to a constant value
of 1.0.  It runs the r2c/c2r FFTW advanced planner, accouting (or
attempting to account) for the requisite doubling of the output stride
and halving of the output distance for the r2c transform.  It computes
the forward FFT and shows the first three transformed vectors in the
output.  As the only nonzero Fourier component should be the DC, the
program sums the DC components of each transformed vector and sums the
remaining components.  Then the program computes the inverse FFT and
sums the output, checking consistency with input.

!-----------------------------------------------------------------------
program FFTW_2D_R2C_Decomp

  implicit none

  include 'fftw3.f'

  integer   :: n1, n1p, n2, s1, d1, h1
  integer*8 :: pf1, pb1, pf2, pb2
  real, dimension(:,:), allocatable :: data

  ! n1-> "fast" axis, n2-> "slow" axis, n1p-> padding on "fast" axis
  n1=300
  n2=500
  n1p=2*(n1/2+1)

  allocate( data(n1p,n2) )

  ! Input data = constant value of 1.0 for all x,y
  data(1:n1,:) = 1.0

  ! "howmany", "dist", and "stride" for FFTW advanced interface
  h1 = n2;  d1 = n1p;  s1 = 1

  !---------------------------------------------------------------------
  ! First run forward/inverse 1-D FFT along "fast" axis of 2-D data
  !---------------------------------------------------------------------
  call sfftw_plan_many_dft_r2c(pf1, 1, n1, h1, data, n1, s1, d1, data,
&
                               n1, s1*2, d1/2, FFTW_ESTIMATE)
  call sfftw_plan_many_dft_c2r(pb1, 1, n1, h1, data, n1, s1*2, d1/2,
&
                               data, n1, s1, d1, FFTW_ESTIMATE)
  call sfftw_execute(pf1)
  write(0,*) '--------------------------------------------'
  write(0,*) 'Running 1-D Forward FFT of 1-axis of 2D data'
  write(0,*) '--------------------------------------------'
  write(0,'(a,10i3)') ' 5 (real,imag) pairs for 3 1-axis FFTd vectors;
each should be: ',n1,0,0,0,0,0,0,0,0,0
  write(0,'(" ",10f6.2)') data(1:10,1)
  write(0,'(" ",10f6.2)') data(1:10,2)
  write(0,'(" ",10f6.2)') data(1:10,3)
  write(0,'(" DC sum (should be",i7,"):",f10.2)') n1*n2,sum(data(1,:))
  write(0,'(" non-DC sum (should be 0):",f10.2)') sum(data(2:,:))
  call sfftw_execute(pb1)
  data = data/n1
  write(0,*) '--------------------------------------------'
  write(0,*) 'Running 1-D Inverse FFT of 1-axis of 2D data'
  write(0,*) '--------------------------------------------'
  write(0,'(a,5i3)') ' 5 real values for 3 1-axis reconstructed
vectors; each should be: ',1,1,1,1,1
  write(0,'(" ",5f6.2)') data(1:5,1)
  write(0,'(" ",5f6.2)') data(1:5,2)
  write(0,'(" ",5f6.2)') data(1:5,3)
  write(0,'(" sum (should be",i7,"):",f10.2)') n1*n2,sum(data)
  call sfftw_destroy_plan(pf1)
  call sfftw_destroy_plan(pb1)

  deallocate(data)

  call exit(0)

end program FFTW_2D_R2C_Decomp
!-------------------------------------------------------------------------------

If you save this program as FFTW_2D_R2C_Decomp.f90, you can compile it
with the following ifort command, provided you define "/usr/local/
include" and "/usr/local/lib" to point to your FFTW include and lib
directories, respectively.

ifort  -O3 -vms -w -WB -I/usr/local/include FFTW_2D_R2C_Decomp.f90 -L/
usr/local/lib -lfftw3f -lfftw3 -lfftw3f_threads -lpthread -lm -o

When I run ./FFTW_2D_R2C_Decomp, I obtain the following:

--------------------------------------------
Running 1-D Forward FFT of 1-axis of 2D data
--------------------------------------------
5 (real,imag) pairs for 3 1-axis FFTd vectors; each should be: 300  0
0  0  0  0  0  0  0  0
300.00  0.00  1.00  1.00  0.00  0.00  1.00  1.00  0.00  0.00
300.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
300.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
DC sum (should be 150000): 150000.00
non-DC sum (should be 0):    150.00
--------------------------------------------
Running 1-D Inverse FFT of 1-axis of 2D data
--------------------------------------------
5 real values for 3 1-axis reconstructed vectors; each should be:   1
1  1  1  1
  1.00  1.00  1.00  1.00  1.00
  1.00  1.00  1.00  1.00  1.00
  1.00  1.00  1.00  1.00  1.00
sum (should be 150000): 149999.98
--------------------------------------------

In other words, the transformed vectors look fine, except for the
first vector, which appears to have something resembling the FT of a
cosine around Nyquist.  The DC components have the correct magnitude.
Oddly, the inverse FFT perfectly reconstructs the input.

To pre-empt some questions: yes, I have RTFM'd ;-). I have run the
same test program using complex-to-complex FFTs and obtain the correct
result.  I was even able to use the "guru" planner to compute the c2c
FFT of the 2-axis of a 3D array!

Although I have read the salient sections of the FFTW manual, from my
results, I suspect that I've still got an "off by 1" bug in one of my
dimensions.  I'd appreciate it if someone who's successfully
implemented this transform pair could help me out.

Morgan Brown