module fft_data_mod

!  use down_process_types
  use fft_types
!  use timers_mod

  implicit none

  ! GOOD For SGI
  integer, parameter ::  FFT_FORW=1
  integer, parameter ::  FFT_BACK=-1

  ! GOOD For LINUX
  !integer, parameter ::  FFT_FORW=-1
  !integer, parameter ::  FFT_BACK=1


  ! Parameters for FFTW package pasted from fftw3.f

      INTEGER FFTW_R2HC
      PARAMETER (FFTW_R2HC=0)
      INTEGER FFTW_HC2R
      PARAMETER (FFTW_HC2R=1)
      INTEGER FFTW_DHT
      PARAMETER (FFTW_DHT=2)
      INTEGER FFTW_REDFT00
      PARAMETER (FFTW_REDFT00=3)
      INTEGER FFTW_REDFT01
      PARAMETER (FFTW_REDFT01=4)
      INTEGER FFTW_REDFT10
      PARAMETER (FFTW_REDFT10=5)
      INTEGER FFTW_REDFT11
      PARAMETER (FFTW_REDFT11=6)
      INTEGER FFTW_RODFT00
      PARAMETER (FFTW_RODFT00=7)
      INTEGER FFTW_RODFT01
      PARAMETER (FFTW_RODFT01=8)
      INTEGER FFTW_RODFT10
      PARAMETER (FFTW_RODFT10=9)
      INTEGER FFTW_RODFT11
      PARAMETER (FFTW_RODFT11=10)
      INTEGER FFTW_FORWARD
      PARAMETER (FFTW_FORWARD=-1)
      INTEGER FFTW_BACKWARD
      PARAMETER (FFTW_BACKWARD=+1)
      INTEGER FFTW_MEASURE
      PARAMETER (FFTW_MEASURE=0)
      INTEGER FFTW_DESTROY_INPUT
      PARAMETER (FFTW_DESTROY_INPUT=1)
      INTEGER FFTW_UNALIGNED
      PARAMETER (FFTW_UNALIGNED=2)
      INTEGER FFTW_CONSERVE_MEMORY
      PARAMETER (FFTW_CONSERVE_MEMORY=4)
      INTEGER FFTW_EXHAUSTIVE
      PARAMETER (FFTW_EXHAUSTIVE=8)
      INTEGER FFTW_PRESERVE_INPUT
      PARAMETER (FFTW_PRESERVE_INPUT=16)
      INTEGER FFTW_PATIENT
      PARAMETER (FFTW_PATIENT=32)
      INTEGER FFTW_ESTIMATE
      PARAMETER (FFTW_ESTIMATE=64)
      INTEGER FFTW_ESTIMATE_PATIENT
      PARAMETER (FFTW_ESTIMATE_PATIENT=128)
      INTEGER FFTW_BELIEVE_PCOST
      PARAMETER (FFTW_BELIEVE_PCOST=256)
      INTEGER FFTW_DFT_R2HC_ICKY
      PARAMETER (FFTW_DFT_R2HC_ICKY=512)
      INTEGER FFTW_NONTHREADED_ICKY
      PARAMETER (FFTW_NONTHREADED_ICKY=1024)
      INTEGER FFTW_NO_BUFFERING
      PARAMETER (FFTW_NO_BUFFERING=2048)
      INTEGER FFTW_NO_INDIRECT_OP
      PARAMETER (FFTW_NO_INDIRECT_OP=4096)
      INTEGER FFTW_ALLOW_LARGE_GENERIC
      PARAMETER (FFTW_ALLOW_LARGE_GENERIC=8192)
      INTEGER FFTW_NO_RANK_SPLITS
      PARAMETER (FFTW_NO_RANK_SPLITS=16384)
      INTEGER FFTW_NO_VRANK_SPLITS
      PARAMETER (FFTW_NO_VRANK_SPLITS=32768)
      INTEGER FFTW_NO_VRECURSE
      PARAMETER (FFTW_NO_VRECURSE=65536)
      INTEGER FFTW_NO_SIMD
      PARAMETER (FFTW_NO_SIMD=131072)

!  integer, external :: omp_get_thread_num

contains

  !----------------------------------------------------------------------
  subroutine  init_fft_1d(fft_id,n1,in,out,fast_init)
    type (fft_type), intent(out) :: fft_id
!  function  init_fft_1d(n1,fast_init) result (fft_id)
!    type (fft_type) fft_id
    integer n1,stat,flag
    complex in(:),out(:)
    logical, optional :: fast_init

    if(present(fast_init)) then
       if(fast_init) then
          flag=FFTW_ESTIMATE
       else
          flag=FFTW_PATIENT
       end if
    else
       flag=FFTW_ESTIMATE
    end if

!    write(0,*)'flag with ior',flag

    fft_id%scale=sqrt(1./(n1))
    fft_id%n_passes=1


!   write(0,*)'fft_id',n1

    call sfftw_plan_dft_1d(fft_id%plan_forw,n1,in,out,FFTW_FORWARD,flag)
!   write(0,*)'fft1d: plan_forw=',fft_id%plan_forw

    call sfftw_plan_dft_1d(fft_id%plan_back,n1,in,out,FFTW_BACKWARD,flag)
!   write(0,*)'fft1d: plan_back=',fft_id%plan_back

    if((fft_id%plan_forw == 0 ).or.(fft_id%plan_back == 0)) then
       fft_id%ierr=1
    else
       fft_id%ierr=0
    end if

    fft_id%ndim=1

!write(0,*)'AFTER PLAN plan_forw=',fft_id%plan_forw
!write(0,*)'AFTER PLAN plan_back=',fft_id%plan_back
!
!write(0,*)'AFTER PLAN scale=',fft_id%scale
!write(0,*)'AFTER PLAN ierr=',fft_id%ierr
!write(0,*)'AFTER PLAN ndim=',fft_id%ndim

!$omp flush (fft_id)

!  end function init_fft_1d
  end subroutine init_fft_1d

  !----------------------------------------------------------------------
  subroutine  init_fft_2d(fft_id,n1,n2,in,out,fast_init) 
    type (fft_type), intent(out) :: fft_id
!  function  init_fft_2d(n1,n2,fast_init) result (fft_id)
!    type (fft_type) fft_id
    integer n1,n2,stat,flag
    complex in(:,:),out(:,:)
    logical, optional :: fast_init

    if(present(fast_init)) then
       if(fast_init) then
          flag=FFTW_ESTIMATE
       else
          flag=FFTW_PATIENT
       end if
    else
       flag=FFTW_ESTIMATE
    end if


    fft_id%scale=sqrt(1./(n1*n2))
    fft_id%n_passes=1

    call sfftw_plan_dft_2d(fft_id%plan_forw,n1,n2,in,out,FFTW_FORWARD,flag)
    !write(0,*)'fft2d: plan_forw=',fft_id%plan_forw

    call sfftw_plan_dft_2d(fft_id%plan_back,n1,n2,in,out,FFTW_BACKWARD,flag)
    !write(0,*)'fft2d: plan_back=',fft_id%plan_back

    if((fft_id%plan_forw ==0 ).or.(fft_id%plan_back == 0)) then
       fft_id%ierr=1
    else
       fft_id%ierr=0
    end if

    fft_id%ndim=2

    !write(0,*)'End of init fft2d'

  end subroutine init_fft_2d
!  end function init_fft_2d

  !----------------------------------------------------------------------
  subroutine init_fft_3d(fft_id,n1,n2,n3,in,out,fast_init) 
    type (fft_type), intent(out) :: fft_id
!  function init_fft_3d(n1,n2,n3,fast_init) result (fft_id)
!    type (fft_type) fft_id
    integer n1,n2,n3,stat,flag
    complex in(:,:,:),out(:,:,:)
    logical, optional :: fast_init
!    character (len=128) :: desc_timer


    if(present(fast_init)) then
       if(fast_init) then
!          desc_timer='Init FFTW 3-D (fast_init)'
!          call SEP_timer_init(20,desc_timer,0.,0)
          flag=FFTW_ESTIMATE
       else
!          desc_timer='Init FFTW 3-D (no fast_init)'
!          call SEP_timer_init(20,desc_timer,0.,0)
          flag=FFTW_PATIENT
       end if
    else
       flag=FFTW_ESTIMATE
    end if

    fft_id%scale=sqrt(1./(n1*n2*n3))
    fft_id%n_passes=1

!    call SEP_timer_start(20)
    call sfftw_plan_dft_3d(fft_id%plan_forw,n1,n2,n3,in,out,FFTW_FORWARD,flag)
    call sfftw_plan_dft_3d(fft_id%plan_back,n1,n2,n3,in,out,FFTW_BACKWARD,flag)
!    call SEP_timer_stop(20)

    if((fft_id%plan_forw ==0 ).or.(fft_id%plan_back == 0)) then
       fft_id%ierr=1
    else
       fft_id%ierr=0
    end if

    fft_id%ndim=3

  end subroutine init_fft_3d

  !----------------------------------------------------------------------
  subroutine init_fft_4d(fft_id,n1,n2,n3,n4,in,out,fast_init) 
    type (fft_type), intent(out) :: fft_id
    integer n1,n2,n3,n4,stat,flag
    complex in(:,:,:,:),out(:,:,:,:)
    logical, optional :: fast_init

    character (len=128) :: desc_timer
    integer rank,dim(4)

    rank=4
    dim(1)=n1
    dim(2)=n2
    dim(3)=n3
    dim(4)=n4


    if(present(fast_init)) then
       if(fast_init) then
!          desc_timer='Init FFTW 4-D (fast_init)'
!          call SEP_timer_init(20,desc_timer,0.,0)
          flag=FFTW_ESTIMATE
       else
!          desc_timer='Init FFTW 4-D (no fast_init)'
!          call SEP_timer_init(20,desc_timer,0.,0)
          flag=FFTW_PATIENT
       end if
    else
       flag=FFTW_ESTIMATE
    end if

    fft_id%scale=sqrt(1./(n1*n2*n3*n4))
    fft_id%n_passes=1

!    call SEP_timer_start(20)
    call sfftw_plan_dft(fft_id%plan_forw,rank,dim,in,out,FFTW_FORWARD,flag)
    call sfftw_plan_dft(fft_id%plan_back,rank,dim,in,out,FFTW_BACKWARD,flag)
!    call SEP_timer_stop(20)

    if((fft_id%plan_forw ==0 ).or.(fft_id%plan_back == 0)) then
       fft_id%ierr=1
    else
       fft_id%ierr=0
    end if

    fft_id%ndim=4

!   desc_timer='execution FFTW 4-D Forward' 
!   call SEP_timer_init(21,desc_timer,0.,0)
!   desc_timer='execution FFTW 4-D Bacward' 
!   call SEP_timer_init(22,desc_timer,0.,0)

  end subroutine init_fft_4d

  !----------------------------------------------------------------------
  subroutine fft_1d_data(fft_id,sign,n1,data_vec)
    type (fft_type) fft_id
    integer :: n1
    integer :: sign
    complex, dimension (:) :: data_vec

!    if(omp_get_thread_num() == 0) write(0,*)'Beg fft 1d'
      data_vec=data_vec*fft_id%scale

    if(sign == FFT_FORW) then
!       write(0,*)'Forw 1d'
       call sfftw_execute(fft_id%plan_forw)
    else
!       if(omp_get_thread_num() == 0) write(0,*)'Back 1d'
       call sfftw_execute(fft_id%plan_back)
    end if


    return	
  end subroutine fft_1d_data

  !----------------------------------------------------------------------
  subroutine fft_2d_data(fft_id,sign,n1,n2,data_slice)
    type (fft_type) fft_id
    integer :: n1,n2
    integer :: sign
    complex, dimension (:,:) :: data_slice

    !write(0,*)'Beg fft 2d'
    data_slice=data_slice*fft_id%scale

    if(sign == FFT_FORW) then
       !write(0,*)'Forw 2d'
       call sfftw_execute(fft_id%plan_forw)
    else
       !write(0,*)'Back 2d'
       call sfftw_execute(fft_id%plan_back)
    end if

    return	
  end subroutine fft_2d_data


  !----------------------------------------------------------------------
  subroutine fft_3d_data(fft_id,sign,n1,n2,n3,data_slice)
    type (fft_type) fft_id
    integer :: n1,n2,n3
    integer :: sign
    complex,  dimension (:,:,:) :: data_slice

    !write(0,*)'Beg fft 3d'

    data_slice=data_slice*fft_id%scale

    ! write(0,*)'fft_3d shape(data_slice)=',shape(data_slice)
!write(0,*)'Inside ffft fft_3d maxval(cabs(data_slice))=',maxval(cabs(data_slice))

!    if(omp_get_thread_num() == 0) then
!!if(omp_get_thread_num() == 0) write(0,*)'#31'
!       call SEP_timer_start(21)
!    end if
    if(sign == FFT_FORW) then
       !write(0,*)'Forw 3d'
!if(omp_get_thread_num() == 0) write(0,*)'#32'
       call sfftw_execute(fft_id%plan_forw)
!if(omp_get_thread_num() == 0) write(0,*)'#33'
    else
       !write(0,*)'Back 3d'
       call sfftw_execute(fft_id%plan_back)
    end if
!    if(omp_get_thread_num() == 0) then
!!if(omp_get_thread_num() == 0) write(0,*)'#34'
!       call SEP_timer_stop(21)
!    end if

    ! write(0,*)'after fft_3d maxval(cabs(data_slice))=',maxval(cabs(data_slice))

    return	
  end subroutine fft_3d_data

  !----------------------------------------------------------------------
  subroutine fft_3d_two_pass_data(fft_1d_id,fft_2d_id,sign,n1,n2,n3,data_slice)
    type (fft_type) fft_1d_id,fft_2d_id
    integer :: n1,n2,n3,i1,i2,i3
    integer :: sign
    real    :: scale
    complex,  dimension (:,:,:) :: data_slice
!    complex,  dimension (n1,n2) :: temp_2d
!    complex,  dimension (n3)    :: temp_1d

    call seperr('fft_3d_two_pass_data not implemented yet in FFTW3')

!    scale=fft_2d_id%scale*fft_1d_id%scale
!
!    !temp_2d(:,:) = data_slice(:,:,1)
!    !call fftwnd_f77_one(fft_2d_id%plan_back, temp_2d(:,:), data_slice(:,:,i3))
!    !write(0,*)'Passed it 2D !'
!
!    !temp_1d(:) = data_slice(1,1,:)
!    !call fftw_f77_one(fft_1d_id%plan_back, temp_1d(:), temp_data(:))
!    !write(0,*)'Passed it 1D !'
!
!    if(sign == FFT_FORW) then
!
!       do i3=1,n3
!          temp_2d = data_slice(:,:,i3)*scale
!          call fftwnd_f77_one(fft_2d_id%plan_forw, temp_2d, data_slice(:,:,i3))
!       end do
!       if(n3 > 1) then
!          do i2=1,n2
!             do i1=1,n1
!                temp_1d = data_slice(i1,i2,:)
!                call fftw_f77_one(fft_1d_id%plan_forw, temp_1d, &
!                data_slice(i1,i2,:))
!             end do
!          end do
!       end if
!
!    else
!
!       do i3=1,n3
!          temp_2d = data_slice(:,:,i3)*scale
!          call fftwnd_f77_one(fft_2d_id%plan_back, temp_2d, data_slice(:,:,i3))
!       end do
!       if(n3 > 1) then
!          do i2=1,n2
!             do i1=1,n1
!                temp_1d = data_slice(i1,i2,:)
!                call fftw_f77_one(fft_1d_id%plan_back, temp_1d, &
!                data_slice(i1,i2,:))
!             end do
!          end do
!       end if
!
!    end if
!
    return	
  end subroutine fft_3d_two_pass_data

  !----------------------------------------------------------------------
  subroutine fft_4d_data(fft_id,sign,n1,n2,n3,n4,data_slice)
    type (fft_type) fft_id
    integer :: n1,n2,n3,n4
    integer :: sign
    complex,  dimension (:,:,:,:) :: data_slice

    !write(0,*)'Beg fft 3d'

    data_slice=data_slice*fft_id%scale

    ! write(0,*)'fft_4d shape(data_slice)=',shape(data_slice)
!write(0,*)'Inside ffft fft_4d maxval(cabs(data_slice))=',maxval(cabs(data_slice))

!    if(omp_get_thread_num() == 0) then
!!if(omp_get_thread_num() == 0) write(0,*)'#31'
!    end if
    if(sign == FFT_FORW) then
       !write(0,*)'Forw 3d'
!if(omp_get_thread_num() == 0) write(0,*)'#32'
!       call SEP_timer_start(21)
       call sfftw_execute(fft_id%plan_forw)
!       call SEP_timer_stop(21)
!if(omp_get_thread_num() == 0) write(0,*)'#33'
    else
       !write(0,*)'Back 3d'
!       call SEP_timer_start(22)
       call sfftw_execute(fft_id%plan_back)
!       call SEP_timer_stop(22)
    end if
!    if(omp_get_thread_num() == 0) then
!!if(omp_get_thread_num() == 0) write(0,*)'#34'
!       call SEP_timer_stop(21)
!    end if

    ! write(0,*)'after fft_4d maxval(cabs(data_slice))=',maxval(cabs(data_slice))

    return	
  end subroutine fft_4d_data


  !----------------------------------------------------------------------
  subroutine fft_4d_two_pass_data(fft_1d_id,fft_3d_id,sign, &
  n1,n2,n3,n4,data_slice)
    type (fft_type), intent(in) :: fft_1d_id,fft_3d_id
    integer, intent(in) :: n1,n2,n3,n4
    integer, intent(in) :: sign
    complex, intent(inout),  dimension (:,:,:,:)  :: data_slice

    integer, automatic :: i1,i2,i3,i4
    real, automatic    :: scale
!    complex, automatic,  dimension (n1,n2,n3) :: temp_3d
!    complex, automatic,  dimension (n4)       :: temp_1d

!write(0,*)'n1=',n1
!write(0,*)'n2=',n2
!write(0,*)'n3=',n3
!write(0,*)'n4=',n4
!write(0,*)'INSIDE'
!write(0,*)'shape(temp_3d)=',shape(temp_3d)
!write(0,*)'shape(temp_1d)=',shape(temp_1d)

    call seperr('fft_4d_two_pass_data not implemented yet in FFTW3')

!    scale=fft_3d_id%scale*fft_1d_id%scale
!!write(0,*)'#31 scale=',scale,'thread#',omp_get_thread_num()
!
!    if(sign == FFT_FORW) then
!
!!if(omp_get_thread_num() == 0) write(0,*)'#31'
!
!       do i4=1,n4
!          temp_3d = data_slice(:,:,:,i4)*scale
!!write(0,*)'#331 scale=',scale,'thread#',omp_get_thread_num()
!          call fftwnd_f77_one(fft_3d_id%plan_forw, temp_3d, data_slice(:,:,:,i4))
!       end do
!
!!if(omp_get_thread_num() == 0) write(0,*)'#32'
!!write(0,*)'#32'
!
!       if(n4 > 1) then
!          do i3=1,n3
!             do i2=1,n2
!                do i1=1,n1
!                   temp_1d = data_slice(i1,i2,i3,:)
!                   call fftw_f77_one(fft_1d_id%plan_forw, temp_1d, &
!                   data_slice(i1,i2,i3,:))
!                end do
!             end do
!          end do
!       end if
!!if(omp_get_thread_num() == 0) write(0,*)'#33'
!!write(0,*)'#33'
!
!    else   
!
!       do i4=1,n4
!          temp_3d = data_slice(:,:,:,i4)*scale
!          call fftwnd_f77_one(fft_3d_id%plan_back, temp_3d, data_slice(:,:,:,i4))
!       end do
!       if(n4 > 1) then
!          do i3=1,n3
!             do i2=1,n2
!                do i1=1,n1
!                   temp_1d = data_slice(i1,i2,i3,:)
!                   call fftw_f77_one(fft_1d_id%plan_back, temp_1d, &
!                   data_slice(i1,i2,i3,:))
!                end do
!             end do end do
!       end if
!
!    end if
!
!    return	

  end subroutine fft_4d_two_pass_data


  !----------------------------------------------------------------------
  subroutine clean_fft(fft_id) 
    type (fft_type) fft_id
    integer ierr

    select case (fft_id%ndim)
    case (0)
       write(0,*)'clean up all FFTW'
       call sfftw_cleanup()
    case (1)
       !write(0,*)'Beg clean fft 1d'
       call sfftw_destroy_plan(fft_id%plan_forw)
       call sfftw_destroy_plan(fft_id%plan_back)
       !write(0,*)'End clean fft 1d'
    case (2)
       !write(0,*)'Beg clean fft 2d'
       call sfftw_destroy_plan(fft_id%plan_forw)
       call sfftw_destroy_plan(fft_id%plan_back)
       !write(0,*)'End clean fft 2d'
    case (3)
       !write(0,*)'Beg clean fft 3d'
       call sfftw_destroy_plan(fft_id%plan_forw)
       call sfftw_destroy_plan(fft_id%plan_back)
       !write(0,*)'End clean fft 3d'
    case default
       write(0,*)'fft_id%ndim=',fft_id%ndim
       call seperr('Not supported fft_id%ndim')
    end select


    return	
  end subroutine clean_fft

end module fft_data_mod
