  PROGRAM JACOBI_MATRIX
 
  Use MKL_DFTI
  IMPLICIT NONE
  INCLUDE 'MKL_RCI.fi'
 
! EXTERNAL FCN
  INTEGER N,M,k
  PARAMETER(N=10,M=1)
  double precision Fun,jac_eps,gnew(1,N)
  real g(N),res
  	
  !   FFT MKL DECLARATIONS
  double complex, dimension(:),allocatable      :: data_in
  type(DFTI_DESCRIPTOR), POINTER :: Desc_Handle_FFT
  integer   Status_FFT,length_FFT(2),dim1,dim2
  
  !	FUNCTION YN
  double precision::X(N)
  
  ! OMP VARIABLES	
  integer							  :: nprocs,OMP_GET_NUM_PROCS,id_thread,nthread,OMP_GET_THREAD_NUM,OMP_GET_NUM_THREADS
  
  !Timing variables
  real(8)                             ::d1,d2,dclock
!-------------------------------------------------------------------------

!-------------------------------------------------------------------------
    nthread=0
    !$omp parallel shared(nthread)
	id_thread=OMP_GET_THREAD_NUM()
	
	if (id_thread.eq.0) then	
	nthread=OMP_GET_NUM_THREADS()
	print*, 'NUMBER OF ACTIVE THREADS ---->', nthread
	endif

	!$omp end parallel
    
!	CREATE FFT DESCRIPTOR 
    allocate(data_in(N))
    length_FFT(1)=N;length_FFT(2)=1
    ! Create Descriptor
    Status_FFT=DftiCreateDescriptor( Desc_Handle_FFT, DFTI_DOUBLE,DFTI_COMPLEX,2,length_FFT)
    ! Commit Descriptor
    Status_FFT=DftiCommitDescriptor(Desc_Handle_FFT)
	
!	INITIALIZE X
	X(1:3)=0.05*70;X(4:7)=0.1*70;X(8:10)=0.15*70
	
!	EVALUATE FCN
	call FCN(M,N,X,Fun)
	print *,'Starting Fun =',Fun
	
	d1=dclock()		  
!   Evaluate Gradient	
    gnew=0;jac_eps=1D-4;
    do k=1,1000
    gnew=0.
    res=djacobi(FCN,N,1,gnew(1,1:N),X(1:N),jac_eps)
    end do
!   Print gradient  
    print *,gnew(1,1:N)
!	Print gradient evaluation time
    d2=dclock()		  
    print *,'GRADIENT EVALUATION TIME=',d2-d1  

 
    CONTAINS

! ***********************************************************************    
        subroutine EVAL_FFT(X_IN,N_X,OUT_VETT,N_OUT)
!       Perform the FFT of (exp(j*X_IN))
    
        implicit none
        integer N_X,N_OUT
        double precision X_IN(N_X)
        complex imag
        parameter (imag=(0.,1.))
        real dummy(N_X)
        complex OUT_VETT(N_OUT)
!-----------------------------------    
	    ! Data_in=exp(j*X_IN)
	    dummy(1:N_X)=X(1:N_X)
        data_in=cexp(imag*dummy)
    
 	    ! OUT_VETT=FFT(data_in)     
        Status_FFT=DftiComputeForward(Desc_Handle_FFT,data_in)
        OUT_VETT=data_in
    
        end subroutine EVAL_FFT
        
! ***********************************************************************        
	    subroutine FCN(M,N,X_IN,F)
    
        implicit none
        integer M,N,M_out
	    external sdot
	    double precision X_IN(N),F
	    real m1(N),ps1,sdot
	    complex out(N)
!-----------------------------------    
       	M_out=N
		! perform the FFT of exp(j*X_IN)
        call EVAL_FFT(X_IN,N,out,M_out)

		! normalize the output
        out=out/maxval(cabs(out))
		
		! perform the sum of the element of abs(out)
		F=sum(cabs(out))


    	end subroutine FCN
! ***********************************************************************

 
    end program JACOBI_MATRIX
 
