module wei_fft_mod
  use wei_fft_fftw_mod
  use wei_fft_mkl_mod
  use wei_step_mod
  
  implicit none
  character(len=128), save , private :: method
  real, private,save :: fract



contains

  subroutine wei_fft_init_ser(meth,wfld)
    character(len=*)   :: meth
    complex            :: wfld(:,:,:,:,:)

#ifdef DEBUG 
    write(0,*) "IN FFT INIT SER"
#endif
    method=meth
    fract=sqrt(real(size(wfld,5))/real(size(wfld)))
    if(method(1:3)=="MKL") then
#ifdef MKL
      call weifft_mkl_init_ser(wfld)
#endif
    else if(method(1:4)=="FFTW") then
#ifdef FFTW
      call fftw_init_ser(wfld)
#endif
    else
      call seperr("Unknown FFT method")
    end if


#ifdef DEBUG 
    write(0,*) "OUT FFT INIT SER"
#endif
  end subroutine

  subroutine wei_fft_init_par(step)
    type(westep)       :: step

#ifdef DEBUG 
    write(0,*) "IN FFT INIT PAR"
#endif


    if(method(1:3)=="MKL") then
#ifdef MKL
      call weifft_mkl_init_par(step)
#endif
    else if(method(1:4)=="FFTW") then
#ifdef FFTW
      call fftw_init_par(step)
#endif
    else
      call seperr("Unknown FFT method")
    end if

#ifdef DEBUG 
    write(0,*) "OUT FFT INIT PAR"
#endif
  end subroutine


  subroutine wei_ffft(wfld,step)
    complex      :: wfld(:,:,:,:,:)
    type(westep)  :: step 
    integer        :: i1,i2,i3,i4

#ifdef DEBUG 
    write(0,*) "IN WEI FORWARD FFT"
#endif

    if(method(1:3)=="MKL") then
#ifdef MKL
      call wei_ffft_mkl(wfld,step)
#endif
    else if(method(1:4)=="FFTW") then
#ifdef FFTW
      call wei_ffft_fftw(wfld,step)
#endif
    else
      call seperr("Unknown FFT method")
    end if
    do i4=1,size(wfld,4)
      do i3=1,size(wfld,3)
        do i2=1,size(wfld,2)
          do i1=1,size(wfld,1)
            wfld(i1,i2,i3,i4,step%ithread)=wfld(i1,i2,i3,i4,step%ithread)*fract
          end do
        end do
      end do
    end do

#ifdef DEBUG 
    write(0,*) "OUT WEI FORWARD FFT"
#endif

  end subroutine


    
  subroutine wei_ifft(wfld,step)
    complex      :: wfld(:,:,:,:,:)
    type(westep)  :: step 
    integer       :: i1,i2,i3,i4

#ifdef DEBUG 
    write(0,*) "IN WEI INVERSE FFT"
#endif

    if(method(1:3)=="MKL") then
#ifdef MKL
      call wei_ifft_mkl(wfld,step)
#endif
    else if(method(1:4)=="FFTW") then
#ifdef FFTW
      call wei_ifft_fftw(wfld,step)
#endif
    else
      call seperr("Unknown FFT method")
    end if
    do i4=1,size(wfld,4)
      do i3=1,size(wfld,3)
        do i2=1,size(wfld,2)
          do i1=1,size(wfld,1)
            wfld(i1,i2,i3,i4,step%ithread)=wfld(i1,i2,i3,i4,step%ithread)*fract
          end do
        end do
      end do
    end do


#ifdef DEBUG 
    write(0,*) "OUT WEI INVERSE FFT"
#endif

  end subroutine
end module
