program main
    implicit none
    integer,parameter :: FTSIZE=1048576
    integer :: tic_start,tic_finish,tic_rate
    real*8 :: t
    real*8,parameter :: M_PI=3.141592653589793238460D0
    integer :: i
    complex*16,dimension(:),allocatable:: X,B
    allocate(X(0:FTSIZE-1))
    allocate(B(0:FTSIZE-1))
    do i=0,FTSIZE-1
        X(i)=cmplx(1D0/(i+1),1D0/(FTSIZE-i),kind(X(i)));
    end do
    print '("N=",I0)',FTSIZE
    call tic
    call fft(X,B)
    t=toc()
    print '("B(0)=(",G0," ",G0,")")',B(0)
    print '("fft took ",G0," seconds.")',t
!$omp parallel
!$omp single
    call tic
    call pfft(X,B)
    t=toc()
!$omp end single
!$omp end parallel
    print '("B(0)=(",G0," ",G0,")")',B(0)
    print '("parallel fft took ",G0," seconds.")',t
    return

contains 

subroutine tic
    call system_clock(count_rate=tic_rate)
    call system_clock(tic_start)
end subroutine

real*8 function toc()
    call system_clock(tic_finish)
    toc=float(tic_finish-tic_start)/tic_rate
end function

recursive subroutine fftwork(N,s,o,x,p,b)
    integer :: k,N2,N,s,o,p
    complex*16,dimension(0:) :: x,b
    complex*16 :: even,odd,w
    if(N.eq.1) then
        b(p)=x(o)
        return
    end if
    if(mod(N,2).ne.0) then
        print '("N not divisible by 2!")'
        stop
    end if
    N2=N/2
    call fftwork(N2,2*s,o,x,p,b)
    call fftwork(N2,2*s,o+s,x,p+N2,b)
    do k=0,N2-1
        even=b(p+k)
        odd=b(p+k+N2)
        w=exp(cmplx(0D0,-2*M_PI*k/N,kind(w)))
        b(p+k)=even+w*odd
        b(p+k+N2)=even-w*odd
    end do
end subroutine

subroutine fft(x,b)
    complex*16,dimension(0:) :: x,b
    integer :: N
    N=size(x)
    call fftwork(N,1,0,x,0,b)
end subroutine

recursive subroutine pfftwork(N,s,o,x,p,b)
    integer :: k,N2,N,s,o,p
    complex*16,dimension(0:) :: x,b
    complex*16 :: even,odd,w
    if(N.eq.1) then
        b(p)=x(o)
        return
    end if
    if(mod(N,2).ne.0) then
        print '("N not divisible by 2!")'
        stop
    end if
    N2=N/2
!$omp task default(shared) final(s<=32)
    call pfftwork(N2,2*s,o,x,p,b)
!$omp end task
    call pfftwork(N2,2*s,o+s,x,p+N2,b)
!$omp taskwait
    do k=0,N2-1
        even=b(p+k)
        odd=b(p+k+N2)
        w=exp(cmplx(0D0,-2*M_PI*k/N,kind(w)))
        b(p+k)=even+w*odd
        b(p+k+N2)=even-w*odd
    end do
end subroutine

subroutine pfft(x,b)
    complex*16,dimension(0:) :: x,b
    integer :: N
    N=size(x)
    call pfftwork(N,1,0,x,0,b)
end subroutine

subroutine cvecprint(x)
    complex*16,dimension(0:) :: x
    integer :: N,i
    N=size(x)
    do i=0,N-1
        print '("(",G0," ",G0,")")',x(i)
    end do
end subroutine

end program
