SUBROUTINE cfft3dsig (f, nx, ny, nz, ldx, ldy, ldz, isign)

!HL This is the driver routine for a 3d complex fft of lengths nx ,ny, nz
!     isign > 0 : f(G) => f(R)    
!     isign < 0 : f(R) => f(G)
!     This is modeled on the cfft3d routine in fft_scalar.f90 and which is used to FFT
!     the potential, and the density. It creates a plan and executes the plan for the 
!     self energy on the smaller real space grid.  
!     This is designed to be used with FFTW3 all the other possible FFT libraries are implement
!     in fft scalar. 

     IMPLICIT NONE
    
     PUBLIC 

     INTEGER, INTENT(IN) :: nx, ny, nz, ldx, ldy, ldz, isign
     COMPLEX (DP) :: f(:)
     INTEGER :: i, k, j, err, idir, ip
     REAL(DP) :: tscale
     INTEGER, SAVE :: icurrent = 1
     INTEGER, SAVE :: dims(3,ndims) = -1

     C_POINTER, SAVE :: fw_plan ( 3, ndims ) = 0
     C_POINTER, SAVE :: bw_plan ( 3, ndims ) = 0

     IF ( nx < 1 ) &
         call errore('cfft3',' nx is less than 1 ', 1)
     IF ( ny < 1 ) &
         call errore('cfft3',' ny is less than 1 ', 1)
     IF ( nz < 1 ) &
         call errore('cfft3',' nz is less than 1 ', 1)

     ip = -1
     DO i = 1, ndims

       !   first check if there is already a table initialized
       !   for this combination of parameters

       IF ( ( nx == dims(1,i) ) .and. &
            ( ny == dims(2,i) ) .and. &
            ( nz == dims(3,i) ) ) THEN
         ip = i
         EXIT
       END IF
     END DO

     IF( ip == -1 ) THEN
       !   no table exist for these parameters
       !   initialize a new one

       IF ( nx /= ldx .or. ny /= ldy .or. nz /= ldz ) &
            call errore('cfft3dsig','not implemented',3)

       IF( fw_plan(icurrent) /= 0 ) CALL dfftw_destroy_plan( fw_plan(icurrent) )
       IF( bw_plan(icurrent) /= 0 ) CALL dfftw_destroy_plan( bw_plan(icurrent) )

       idir = -1

       CALL dfftw_plan_dft_3d ( fw_plan(icurrent), nx, ny, nz, f(1:), &
            f(1:), idir, FFTW_ESTIMATE)

       idir =  1

       CALL dfftw_plan_dft_3d ( bw_plan(icurrent), nx, ny, nz, f(1:), &
            f(1:), idir, FFTW_ESTIMATE)

       IF (err /= 0) CALL errore('cfft3d','FFT init returned an error ', err)

       dims(1,icurrent) = nx; dims(2,icurrent) = ny; dims(3,icurrent) = nz
       ip = icurrent
       icurrent = MOD( icurrent, ndims ) + 1

     END IF

     !
     !   Now perform the 3D FFT using the machine specific driver
     !

    IF( isign < 0 ) THEN
      call dfftw_execute_dft( fw_plan(ip), f(1:), f(1:))
      tscale = 1.0_DP / DBLE( nx * ny * nz )
      call ZDSCAL( nx * ny * nz, tscale, f(1), 1)
    ELSE IF( isign > 0 ) THEN
       call dfftw_execute_dft( bw_plan(ip), f(1:), f(1:))
    END IF

   IF (err /= 0) CALL errore('cfft3d','FFT returned an error ', err)

   DEALLOCATE(cw2)

   RETURN

END SUBROUTINE cfft3dsig