    program FFT_2D_openmp
    
    !This program is intedended as an excercise to verify how the MKL FFT routines can 
    !be used within a cycle made parallel by using the OPENMP
    
    !The code is essentially made by a cycle. Within the cycle tha matrix to be transformed is 
    !inizialized and then transformed. The sum of the elements of the transformed matrix is used
    !to verify the correct implementation with/without OpenMP
        
    Use MKL_DFTI
    
    implicit none
    
    real pi,Nx,Ny,M
    parameter (pi=3.14159265,Nx=300,Ny=300,M=100)
    
    integer                           ::shape(2),xx
    complex                           ::data_fft(Nx*Ny),sum_vect(M),sum_serial,sum_parallel
    complex                           ::imag
!   Declaration of the FFT MKL parameters   
    type(DFTI_DESCRIPTOR), POINTER :: Desc_Handle
    integer   status
    real      Scale
    integer   L(2)
    integer   strides_in(3)
    !  Timing variables
	real(8)						:: d1,d2,dclock
	real(8)						:: time_serial,time_parallel,speedup,efficiency
    !  Openmp variables
    integer						::id_thread,nthread,OMP_GET_THREAD_NUM,OMP_GET_NUM_THREADS,num_it_schedule    
    integer                     :: MKL_FFT
!------------------------------------------------------------------
    imag=(0.,1.)
    ! Set size of the FFT
    L(1) = Nx
    L(2) = Ny
    
    
        !$omp parallel
        nthread=OMP_GET_NUM_THREADS()   !number of threads
        !$omp end parallel
        
    ! Create Descriptor
    Status=DftiCreateDescriptor( Desc_Handle, DFTI_SINGLE,DFTI_COMPLEX, 2,L)
    ! Set numthread for FFT
     Status=DftiSetValue(Desc_Handle,DFTI_NUMBER_OF_USER_THREADS,nthread)
    ! Commit Descriptor
    Status=DftiCommitDescriptor(Desc_Handle)
    


!   SERIAL RUN: Serial time evaluation
 
       d1=dclock() !time start evaluation
    
        !Start Cycle
        do xx=1,M
    
            !initalize data to be transformed
            data_fft=xx+imag*xx*2.
        
            ! Perform FFT
            Status=DftiComputeForward(Desc_Handle,data_fft)
        
            !perform sum of the elements
            sum_vect(xx)=sum(data_fft)
        
        end do
       d2=dclock() !time stop evaluation
       
       time_serial=(d2-d1)
       sum_serial=sum(sum_vect(1:M))

!   END SERIAL RUN    

!-------------------------------------


!-------------------------------------
!   Parallel RUN: Parallel time evaluation
    
        num_it_schedule=M/nthread
        d1=dclock() !time start evaluation
        
        !Start Cycle
        
       !$omp parallel
       !$omp do private(data_fft,Status,Desc_Handle) schedule(static,num_it_schedule)
        do xx=1,M
             
            !initalize data to be transformed
            data_fft=xx+imag*xx*2.
        
            ! Perform FFT
            Status=DftiComputeForward(Desc_Handle,data_fft)
            
            !perform sum of the elements
            sum_vect(xx)=sum(data_fft)
        
        end do
        !$omp end do
        !$omp end parallel
        
        
       d2=dclock() !time stop evaluation
       
       time_parallel=(d2-d1)
       sum_parallel=sum(sum_vect(1:M))

!   END PARALLEL RUN    

!   Speedup evaluation    
    speedup=time_serial/time_parallel

!   Efficiency evaluation 
    efficiency=time_serial/time_parallel/nthread

    
    
        
      !   Print execution Data
    print *,'---------- S   t   a   r   t'
    Print *,'+  Matrix size Nx,Ny =',Nx,NY
    Print *,'+  Cycle over M  =',M
    print *,'--> SERIAL'
    Print *,'+  Serial execution time =',time_serial
    Print *,'+  Serial result (sum) =',sum_serial
    print *,'--> PARALLEL'
    print *,'+  Number of threads =',nthread
    Print *,'+  Parallel execution time =',time_parallel
    Print *,'+  Parallel result (sum) =',sum_parallel
    print *,'--> SPEEDUP (ideal = nthread) =',speedup
    print *,'--> EFFICIENCY (ideal =1) =',efficiency
    print *,'---------- S   t   o  p'
    
    end program