program ffttest
      implicit none
      real(8) pi
      parameter(pi=3.141592653589793D0)

      integer fftw_real_to_complex,fftw_complex_to_real
      integer fftw_measure, fftw_in_place
      parameter (fftw_real_to_complex=-1,fftw_complex_to_real=1)
      parameter (fftw_measure=1,fftw_in_place=8)
      common /fftwinit/ fplan,bplan
      integer*8 fplan,bplan

      integer n1,n2
      parameter(n1=32,n2=32)
      real(8) scale
      parameter(scale=n1*n2)
      real(8) x(0:n1+1,0:n2-1)
      double complex y(0:n1/2,0:n2-1)
      equivalence (x(0,0),y(0,0))
      real(8) dum

      integer i,j

      write(6,*) 'create plan, this may take a minute'
      call rfftw2d_f77_create_plan(fplan,n1,n2,fftw_real_to_complex,
     &                             fftw_measure + fftw_in_place)
      call rfftw2d_f77_create_plan(bplan,n1,n2,fftw_complex_to_real,
     &                             fftw_measure + fftw_in_place)

      write(6,*) 'compute x(i,j)'
      do j=0,n2-1
        do i=0,n1-1
          x(i,j)=dcos(2.D0*pi/dble(n1)*dble(i))
     &      *dsin(2.D0*pi/dble(n2)*dble(j))
        enddo
      enddo
      write(6,*) 'check a few elements'
      write(6,*) 'x(1,1)=',x(1,1)
      write(6,*) 'x(2,2)=',x(2,2)

c     forward fft
      write(6,*) 'compute forward fft'
      call rfftwnd_f77_one_real_to_complex(fplan,y,dum)
      write(6,*) 'check a few elements'
      write(6,*) 'y(1,1)=',y(1,1)
      write(6,*) 'y(1,n2-1)=',y(1,n2-1)
c     backward fft
      write(6,*) 'compute backward fft'
      call rfftwnd_f77_one_complex_to_real(bplan,x,dum)
c     need to divide by n1*n2 to get original x(i,j),back
      do j=0,n2-1
        do i=0,n1-1
          x(i,j)=x(i,j)/scale
        enddo
      enddo
      write(6,*) 'check elements, should get original x(i,j) back'
      write(6,*) 'x(1,1)=',x(1,1)
      write(6,*) 'x(2,2)=',x(2,2)

      call rfftwnd_f77_destroy_plan(fplan)
      call rfftwnd_f77_destroy_plan(bplan)

      end