Why is this test FFT program giving incorrect results?

Question

I am writing a Fortran code which requires an FFT. I am using the double precision version of four1 from Numerical Recipes in Fortran 77 (page 501). Below is a program to test the FFT. Below that is the output I get from the FFT along with the output I expect to get. The real part of the transform is correct within rounding errors, but the imaginary part is not. However, for my purposes, I do not even need the imaginary part, so should I be able to just grab the real part and continue with that? It still bothers me that the output is not correct though and makes me think I don't understand something about the implimentation of this subroutine.

The way I understand it, the values given to the array "data" when it is constructed are all real (1,1,1,1,0,0,0,0). And the imaginary values are all zero. Is that correct? I am concerned that I am not giving the input array to the FFT in the form that it needs. When I use this subroutine in my actual program, what comes out of the FFT is somewhat nonsensical and the whole thing goes to NaNs after several timesteps.

program fftTest
    implicit none   
    complex(kind=8), dimension(8) :: data = (/1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0/)
    integer :: i

    do i=1,8
        write(*,'("(", F20.15, ",", F20.15, "i )")') data(i)
    end do
print *, '-----'

call dfour1(data,8,1)

    do i=1,8
        write(*,'("(", F20.15, ",", F20.15, "i )")') data(i)
    end do
print *, '-----'

call dfour1(data,8,-1)
data = data/8

do i=1,8
        write(*,'("(", F20.15, ",", F20.15, "i )")') data(i)
    end do

end program fftTest

    SUBROUTINE dfour1(data,nn,isign)
        INTEGER isign,nn
        DOUBLE PRECISION data(2*nn)
        INTEGER i,istep,j,m,mmax,n
        DOUBLE PRECISION tempi,tempr
        DOUBLE PRECISION theta,wi,wpi,wpr,wr,wtemp 
        n=2*nn
        j=1
        do 11 i=1,n,2 !This is the bit reversal section of the routine.
            if(j.gt.i)then
            tempr=data(j) !Exchange the two complex numbers.
            tempi=data(j+1)
            data(j)=data(i)
            data(j+1)=data(i+1)
            data(i)=tempr
            data(i+1)=tempi
        endif
        m=n/2
1       if ((m.ge.2).and.(j.gt.m)) then
            j=j-m
            m=m/2
        goto 1
        endif
        j=j+m
11    continue
      mmax=2 !Here begins the Danielson-Lanczos section of the routine.
2     if (n.gt.mmax) then                   
          istep=2*mmax
          theta=6.28318530717959d0/(isign*mmax)
          wpr=-2.d0*sin(0.5d0*theta)**2
          wpi=sin(theta)
          wr=1.d0
          wi=0.d0
          do 13 m=1,mmax,2  !Here are the two nested inner loops.
              do 12 i=m,n,istep
                 j=i+mmax   !This is the Danielson-Lanczos formula:
                 tempr=wr*data(j)-wi*data(j+1)
                 tempi=wr*data(j+1)+wi*data(j)
                 data(j)=data(i)-tempr
                 data(j+1)=data(i+1)-tempi
                 data(i)=data(i)+tempr
                 data(i+1)=data(i+1)+tempi
12        continue
          wtemp=wr  !Trigonometric recurrence
          wr=wr*wpr-wi*wpi+wr
          wi=wi*wpr+wtemp*wpi+wi
13      continue
        mmax=istep
        goto 2  !Not yet done.
        endif   !All done.
        return 
    END

Expected (correct) output:
( 4.000000000000000, 0.000000000000000i )
( 1.000000000000000, -2.414213562373095i )
( 0.000000000000000, 0.000000000000000i )
( 1.000000000000000, -0.414213562373095i )
( 0.000000000000000, 0.000000000000000i )
( 1.000000000000000, 0.414213562373095i )
( 0.000000000000000, 0.000000000000000i )
( 1.000000000000000, 2.414213562373095i )

Actual output from test program:
( 1.000000000000000, 0.000000000000000i )
( 1.000000000000000, 0.000000000000000i )
( 1.000000000000000, 0.000000000000000i )
( 1.000000000000000, 0.000000000000000i )
( 0.000000000000000, 0.000000000000000i )
( 0.000000000000000, 0.000000000000000i )
( 0.000000000000000, 0.000000000000000i )
( 0.000000000000000, 0.000000000000000i )
-----
( 4.000000000000000, 0.000000000000000i )
( 0.999999999999998, 2.414213562373094i )
( 0.000000000000000, 0.000000000000000i )
( 0.999999999999999, 0.414213562373096i )
( 0.000000000000000, 0.000000000000000i )
( 1.000000000000000, -0.414213562373094i )
( 0.000000000000000, 0.000000000000000i )
( 1.000000000000003, -2.414213562373096i )
-----
( 1.000000000000000, 0.000000000000000i )
( 1.000000000000000, -0.000000000000000i )
( 1.000000000000000, 0.000000000000000i )
( 1.000000000000000, 0.000000000000000i )
( 0.000000000000000, 0.000000000000000i )
( 0.000000000000000, 0.000000000000000i )
( 0.000000000000000, -0.000000000000000i )
( 0.000000000000000, -0.000000000000000i )

Answer 1

program fftTest
    use iso_fortran_env
    implicit none
    integer, parameter :: nn = 8
    real(real64), dimension(2*nn) :: data = [ 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0,   &
         0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 ]
    integer :: i

    do i=1,2*nn,2
        write(*,'("(", F20.15, ",", F20.15, "i )")') data(i),data(i+1)
    end do
print *, '-----'

call dfour1(data,nn,-1)

    do i=1,2*nn,2
        write(*,'("(", F20.15, ",", F20.15, "i )")') data(i),data(i+1)
    end do
print *, '-----'

call dfour1(data,nn,+1)
data = data/8

do i=1,2*nn,2
        write(*,'("(", F20.15, ",", F20.15, "i )")') data(i),data(i+1)
    end do

end program fftTest

and a demo of using the GNU Scientific Library, which is Open Source, via Forran's ISO C Binding:

program fftTest
    use iso_c_binding
    implicit none
    integer (c_size_t), parameter :: nn = 8
    real(C_DOUBLE), dimension(2*nn) :: data = [ 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0,   &
         0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 ]
    integer :: i
    integer :: result

    interface
    !  Adapted from Fortran GSL Library, FGSL,
    !  http://www.lrz.de/services/software/mathematik/gsl/fortran/
      function gsl_fft_complex_radix2_backward(data, stride, n) bind(c)
         import :: c_size_t, c_double, c_int
         integer(c_size_t), value :: stride, n
         real (c_double), dimension(*), intent (inout)  :: data
         integer(c_int) :: gsl_fft_complex_radix2_backward
      end function gsl_fft_complex_radix2_backward
      function gsl_fft_complex_radix2_forward(data, stride, n) bind(c)
         import :: c_size_t, c_double, c_int
         integer(c_size_t), value :: stride, n
         real (c_double), dimension(*), intent (inout)  :: data
         integer(c_int) :: gsl_fft_complex_radix2_forward
      end function gsl_fft_complex_radix2_forward
    end interface

    do i=1,2*nn,2
        write(*,'("(", F20.15, ",", F20.15, "i )")') data(i),data(i+1)
    end do
   print *, '-----'

   result = gsl_fft_complex_radix2_forward (data, 1_c_size_t, nn)
   write (*, '("GSL return code: ", I0)' )  result

    do i=1,2*nn,2
        write(*,'("(", F20.15, ",", F20.15, "i )")') data(i),data(i+1)
    end do
   print *, '-----'

   result = gsl_fft_complex_radix2_backward (data, 1_c_size_t, nn)
   write (*, '("GSL return code: ", I0)' )  result
   data = data/8

   do i=1,2*nn,2
      write(*,'("(", F20.15, ",", F20.15, "i )")') data(i),data(i+1)
   end do

end program fftTest