stand-alone implementations of NAG DFTs C06FPF and C06FRF

my_nag.f90

module my_nag

  use stel_kinds, only: dp

  implicit none

  real(dp), parameter :: twopi = 6.283185307179586_dp

contains

! one-dimensional r2c DFT
SUBROUTINE MY_C06FPF (M, N, X, INIT, TRIG, WORK, IFAIL)

  INTEGER,      intent(in)    :: M     !< number of sequences to transform
  INTEGER,      intent(in)    :: N     !< length of each sequence
  INTEGER,      intent(inout) :: IFAIL
  REAL(dp),     intent(inout) :: X(M*N) ! on input: data to transform; on exit: DFT result
  REAL(dp),     intent(inout) :: TRIG(2*N)
  REAL(dp),     intent(inout) :: WORK(M*N) ! workspace array for in-place transform
  CHARACTER(1), intent(in)    :: INIT

  integer :: i, j, k, src_idx, tgt_idx
  real(dp) :: sqrt_n

  sqrt_n = sqrt(real(N, kind=dp))

  ! loop over all sequences to transform
  do i=1,M

    ! loop over output Fourier coefficients
    do j=0, N/2
      tgt_idx = j*M + i
      work(tgt_idx) = 0.0_dp

      ! loop over input data points
      do k=0, N-1
        src_idx = k*M + i

        work(tgt_idx) = work(tgt_idx) &
          + x(src_idx) * cos(twopi * j*k / real(N, kind=dp))

      end do ! k
      work(tgt_idx) = work(tgt_idx) / sqrt_n
    end do ! j

    do j=1, N/2-1
      tgt_idx = (N-j)*M + i
      work(tgt_idx) = 0.0_dp

      ! loop over input data points
      do k=0, N-1
        src_idx = k*M + i

        work(tgt_idx) = work(tgt_idx) &
          - x(src_idx) * sin(twopi * j*k / real(N, kind=dp))

      end do ! k
      work(tgt_idx) = work(tgt_idx) / sqrt_n
    end do ! j

  end do ! i

  ! copy result in workspace to output array
  x = work

end subroutine

! one-dimensional c2c DFT
SUBROUTINE MY_C06FRF (M, N, X, Y, INIT, TRIG, WORK, IFAIL)

  INTEGER,      intent(in)    :: M
  INTEGER,      intent(in)    :: N
  INTEGER,      intent(inout) :: IFAIL
  REAL(dp),     intent(inout) :: X(M*N)
  REAL(dp),     intent(inout) :: Y(M*N)
  REAL(dp),     intent(inout) :: TRIG(2*N)
  REAL(dp),     intent(inout) :: WORK(2*M*N)
  CHARACTER(1), intent(in)    :: INIT

  integer :: i, j, k, src_idx, tgt_idx_r, tgt_idx_i
  real(dp) :: sqrt_n

  sqrt_n = sqrt(real(N, kind=dp))

  ! loop over all sequences to transform
  do i=1,M

    ! loop over output Fourier coefficients
    do j=0, N-1

      tgt_idx_r = j*M + i
      tgt_idx_i = j*M + i + M*N

      work(tgt_idx_r) = 0.0_dp
      work(tgt_idx_i) = 0.0_dp

      ! loop over input data points
      do k=0, N-1
        src_idx = k*M + i

        work(tgt_idx_r) = work(tgt_idx_r) &
          + x(src_idx) * cos(twopi * j*k / real(N, kind=dp)) &
          + y(src_idx) * sin(twopi * j*k / real(N, kind=dp))

        work(tgt_idx_i) = work(tgt_idx_i) &
          + y(src_idx) * cos(twopi * j*k / real(N, kind=dp)) &
          - x(src_idx) * sin(twopi * j*k / real(N, kind=dp))

      end do ! k

      ! apply normalization
      work(tgt_idx_r) = work(tgt_idx_r) / sqrt_n
      work(tgt_idx_i) = work(tgt_idx_i) / sqrt_n
    end do ! j

  end do ! i

  ! copy result in workspace to output array
  x = work(    1:  M*N)
  y = work(M*N+1:2*M*N)

end subroutine

end module

These are stand-alone implementations of the NAG discrete Fourier transforms described here:
C06FPF
C06FRF