parallel-1d-fast-ft.cpp

#include "fourier-transform-solver/parallel-1d-fast-ft.hpp"
#include "utils/bit-reversal.hpp"
#include <chrono>


void Parallel1DFastFT::compute() {
    /**
     * Parallel Cooley-Tukey Fast Fourier Transform (FFT) Algorithm.
     * The algorithm is intentionally uncommented because it is well commented in the sequential version.
     * Therefore, if you are looking for a well-commented version, please check the sequential version.
     */
    printf("Running the Parallel 1D Fast Fourier Transform solver...\n");
    // 0. Prepare the variables.
    const size_t N = input.size();
    const auto log2N = static_cast<size_t>(std::log2(N));
    constexpr double angle = -2 * M_PI;
    // save time before starting the computation (performance measurement)
    const std::chrono::time_point<std::chrono::system_clock> start_time = std::chrono::high_resolution_clock::now();

    // 1. Bit-Reversal Permutation
    parallelBitReversal(input);

    // 2. Iterative Cooley-Tukey FFT
    // 2.1. Main loop - Stages: iterate log2(N) times.
    for (size_t s = 1; s <= log2N; ++s) {
        const size_t m = 1 << s;  // 2^s
        const size_t m2 = m >> 1; // m/2
        const std::complex<double> wm = std::polar(1.0, angle / static_cast<double>(m));

        /**
         * Loop Parallelization (OpenMP).
         * each iteration of the loop processes a different chunk of size 'm', allowing parallel execution.
         * therefore, each thread processes a different chunk of the input vector.
         */
#pragma omp parallel for
        // 2.2. Inner loop - Iterating over the 'sub-DFTs' in the current stage.
        for (size_t k = 0; k < N; k += m) {
            std::complex<double> w{1.0, 0};
            // 2.3. Butterfly computation.
            for (size_t j = 0; j < m2; ++j) {
                std::complex<double> t = w * input[k + j + m2];
                std::complex<double> u = input[k + j];
                input[k + j] = u + t;
                input[k + j + m2] = u - t;
                w *= wm;
            }
        }
    }

    // save time after the computation is done (performance measurement)
    const std::chrono::time_point<std::chrono::system_clock> end_time = std::chrono::high_resolution_clock::now();

    // save the solution to the output vector
    solution = input;

    printf("Parallel 1D Fast Fourier Transform computed. Elapsed time: %f ms.\n",
        std::chrono::duration<float, std::milli>(end_time - start_time).count());
}