main.cpp

#include "pls/internal/scheduling/scheduler.h"
#include "pls/internal/scheduling/parallel_result.h"
#include "pls/internal/scheduling/scheduler_memory.h"
using namespace pls::internal::scheduling;

#include <iostream>
#include <complex>
#include <vector>
#include <atomic>

static constexpr int CUTOFF = 16;
static constexpr int INPUT_SIZE = 8192;
typedef std::vector<std::complex<double>> complex_vector;

void divide(complex_vector::iterator data, int n) {
  complex_vector tmp_odd_elements(n / 2);
  for (int i = 0; i < n / 2; i++) {
    tmp_odd_elements[i] = data[i * 2 + 1];
  }
  for (int i = 0; i < n / 2; i++) {
    data[i] = data[i * 2];
  }
  for (int i = 0; i < n / 2; i++) {
    data[i + n / 2] = tmp_odd_elements[i];
  }
}

void combine(complex_vector::iterator data, int n) {
  for (int i = 0; i < n / 2; i++) {
    std::complex<double> even = data[i];
    std::complex<double> odd = data[i + n / 2];

    // w is the "twiddle-factor".
    // this could be cached, but we run the same 'data_structures' algorithm parallel/serial,
    // so it won't impact the performance comparison.
    std::complex<double> w = exp(std::complex<double>(0, -2. * M_PI * i / n));

    data[i] = even + w * odd;
    data[i + n / 2] = even - w * odd;
  }
}

void fft_normal(complex_vector::iterator data, int n) {
  if (n < 2) {
    return;
  }

  divide(data, n);
  fft_normal(data, n / 2);
  fft_normal(data + n / 2, n / 2);
  combine(data, n);
}

parallel_result<short> fft(complex_vector::iterator data, int n) {
  if (n < 2) {
    return parallel_result<short>{0};
  }

  divide(data, n);
  if (n <= CUTOFF) {
    fft_normal(data, n / 2);
    fft_normal(data + n / 2, n / 2);
    combine(data, n);
    return parallel_result<short>{0};
  } else {
    return scheduler::par([=]() {
      return fft(data, n / 2);
    }, [=]() {
      return fft(data + n / 2, n / 2);
    }).then([=](int, int) {
      combine(data, n);
      return parallel_result<short>{0};
    });
  }
}

complex_vector prepare_input(int input_size) {
  std::vector<double> known_frequencies{2, 11, 52, 88, 256};
  complex_vector data(input_size);

  // Set our input data to match a time series of the known_frequencies.
  // When applying fft to this time-series we should find these frequencies.
  for (int i = 0; i < input_size; i++) {
    data[i] = std::complex<double>(0.0, 0.0);
    for (auto frequencie : known_frequencies) {
      data[i] += sin(2 * M_PI * frequencie * i / input_size);
    }
  }

  return data;
}

static constexpr int NUM_ITERATIONS = 1000;
constexpr size_t NUM_THREADS = 5;

constexpr size_t NUM_TASKS = 128;

constexpr size_t NUM_CONTS = 128;
constexpr size_t MAX_CONT_SIZE = 512;

int main() {
  complex_vector initial_input = prepare_input(INPUT_SIZE);

  static_scheduler_memory<NUM_THREADS,
                          NUM_TASKS,
                          NUM_CONTS,
                          MAX_CONT_SIZE> static_scheduler_memory;

  scheduler scheduler{static_scheduler_memory, NUM_THREADS};

  auto start = std::chrono::steady_clock::now();
  for (int i = 0; i < NUM_ITERATIONS; i++) {
    complex_vector input_2(initial_input);
    scheduler.perform_work([&]() {
      return scheduler::par([&]() {
        return fft(input_2.begin(), INPUT_SIZE);
      }, []() {
        return parallel_result<int>{0};
      }).then([](int, int) {
        return parallel_result<int>{0};
      });
    });
  }
  auto end = std::chrono::steady_clock::now();
  std::cout << "Framework:  " << std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count()
            << std::endl;

  start = std::chrono::steady_clock::now();
  for (int i = 0; i < NUM_ITERATIONS; i++) {
    complex_vector input_1(initial_input);
    fft_normal(input_1.begin(), INPUT_SIZE);
  }
  end = std::chrono::steady_clock::now();
  std::cout << "Normal:     " << std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count()
            << std::endl;

  return 0;
}