program test_mkl_fftw

    use omp_lib
    implicit none

    include 'fftw3.f'

    integer, parameter :: nfft=32768, nthreads=16, niter=20000
    complex(4), dimension(nfft,nthreads) :: work
    real, dimension(nfft) :: r
    integer(8) :: plan_forward, plan_backward, i, j, k, id

    save work

    write (*,*) nfft, omp_get_num_threads()

    call sfftw_plan_dft_1d(plan_forward, nfft, work(1,1), work(1,1), &
                           FFTW_FORWARD, FFTW_MEASURE)

    call sfftw_plan_dft_1d(plan_backward, nfft, work(1,1), work(1,1), &
                           FFTW_BACKWARD, FFTW_MEASURE)

    call random_seed()
    call random_number(r)

    call omp_set_num_threads(nthreads)

!$omp parallel do private(i,j,k,id)
    do i=1,niter
        id = 1 + omp_get_thread_num()
        if (mod(i,100) .eq. 0) then
            write (*,*) i, id
        endif
        do j=1,nfft
            k = mod(i+j+id,nfft)
            work(j,id) = r(k)
        enddo
        call sfftw_execute_dft(plan_forward, work(1,id), work(1,id))
        call sfftw_execute_dft(plan_backward, work(1,id), work(1,id))
    enddo

end program
