fft_interface.F90 1.99 KB
Newer Older
Daniel Wortmann's avatar
Daniel Wortmann committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
MODULE m_fft_interface

 use m_cfft
 use m_juDFT
#ifdef CPP_FFT_MKL        
 use mkl_dfti
#endif

IMPLICIT NONE
CONTAINS

    subroutine fft_interface(dimen,length,dat,forw)
    ! provides interfaces to ftt subroutines

    integer, intent(in)    :: dimen         !dimension of fft transformation
    integer, intent(in)    :: length(dimen) !length of data in each direction
    complex, intent(inout) :: dat(:)        !data to be transformed, size(dat) should be sum(length)
    logical, intent(in)    :: forw          !.true. for the forward transformation, .false. for the backward one

#ifdef CPP_FFT_MKL        
    type(dfti_descriptor),pointer :: dfti_handle
    integer :: dfti_status
#else
    real,allocatable :: afft(:),bfft(:)
    integer :: isn
#endif
    integer :: size_dat,i
    
       size_dat = 1
       do i = 1,dimen
          size_dat = size_dat * length(i)
       enddo
       if (size(dat) .ne. size_dat) call juDFT_error('array bounds are inconsistent',calledby ='fft_interface')
       if (dimen .ne. 3 )  call juDFT_error('sorry, not implemented yet for this value of dimen',calledby ='fft_interface')
 
#ifdef CPP_FFT_MKL        
       !using MKL library
       dfti_status = DftiCreateDescriptor(dfti_handle,dfti_double,dfti_complex,3,length)
       dfti_status = DftiCommitDescriptor(dfti_handle)
       if (forw) then
           dfti_status = DftiComputeForward(dfti_handle,dat)
       else
           dfti_status = DftiComputeBackward(dfti_handle,dat)
       end if
       dfti_status = DftiFreeDescriptor(dfti_handle)

#else
       allocate(afft(size_dat),bfft(size_dat))
       afft = real(dat)
       bfft = aimag(dat)
       if (forw) then 
          isn = -1
       else
          isn = 1
       end if
       CALL cfft(afft,bfft,size_dat,length(1),length(1),isn)
       CALL cfft(afft,bfft,size_dat,length(2),length(1)*length(2),isn)
       CALL cfft(afft,bfft,size_dat,length(3),size_dat,isn)

       dat = cmplx(afft,bfft)
#endif

    end subroutine fft_interface

END MODULE m_fft_interface