Add new files

cxx/heat/heat_ghosts_mini.cxx

@@ -0,0 +1,167 @@
+//  Copyright (c) 2023 AUTHORS
+//
+//  SPDX-License-Identifier: BSL-1.0
+//  Distributed under the Boost Software License, Version 1.0. (See accompanying
+//  file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#include <iostream>
+#include <thread>
+#include <functional>
+#include <algorithm>
+#include <chrono>
+#include <filesystem>
+#include <fstream>
+#include <assert.h>
+#include <atomic>
+#include <mutex>
+#include <condition_variable>
+
+using std::size_t;
+
+const size_t ghosts = 1;
+size_t nx = 20;//100000;
+const double k = 0.4;
+size_t nt = 1;
+const double dt = 1.0;
+const double dx = 1.0;
+size_t threads = 1;
+
+void pr(const std::vector<double>& total) {
+  std::cout << "[";
+  for(size_t i=0;i<total.size();i++) {
+    std::cout << " " << total[i];
+  }
+  std::cout << " ]" << std::endl;
+}
+
+class Queue {
+  static const size_t sz = 20;
+  double data[sz];
+  size_t head=0, tail=0;
+  std::mutex m;
+  std::condition_variable cv;
+
+public:
+
+  void push(double d) {
+    while(tail - head >= sz) {
+      // this should never happen
+      std::this_thread::yield();
+    }
+    size_t new_tail = tail + 1;
+    data[tail % sz] = d;
+    std::atomic_thread_fence(std::memory_order_seq_cst);
+    tail = new_tail;
+    std::atomic_thread_fence(std::memory_order_seq_cst);
+    if(head+1 == tail) {
+        std::unique_lock lk(m);
+        cv.notify_one();
+    }
+  }
+
+  double pop() {
+    while(head == tail) {
+      std::unique_lock lk(m);
+      cv.wait(lk,[this]()->bool { return this->head < this->tail; });
+      std::this_thread::yield();
+    }
+    double result = data[head % sz];
+    std::atomic_thread_fence(std::memory_order_seq_cst);
+    ++head;
+    return result;
+  }
+};
+
+class Worker : public std::thread {
+public:
+  size_t num;
+  size_t lo, hi, sz;
+  std::vector<double> data, data2;
+  Queue left, right;
+  Worker *leftThread = nullptr;
+  Worker *rightThread = nullptr;
+
+  void start() {
+    std::thread t(std::bind(&Worker::run, this));
+    this->swap(t);
+  }
+
+  Worker(Worker&& w) : lo(w.lo), hi(w.hi), sz(w.sz), data(w.data), data2(w.data2), leftThread(w.leftThread), rightThread(w.rightThread) {}
+
+  Worker(size_t num_, size_t tx) : num(num_) {
+    lo = tx * num;
+    hi = tx * (num + 1);
+    if (hi > nx) hi = nx;
+    lo -= ghosts;
+    hi += ghosts;
+    sz = hi - lo;
+    data.resize(sz);
+    data2.resize(sz);
+
+    size_t off=1;
+    for(size_t n=0; n < sz; n++) {
+      data.at(n) = n + lo + off;
+      data2.at(n) = 0.0;
+    }
+  }
+  //Worker(Worker&&) noexcept = default;
+  ~Worker() {}
+
+  void recv_ghosts() {
+    data[0] = left.pop();
+    data[sz-1] = right.pop();
+  }
+
+  void update() {
+    recv_ghosts();
+
+    for(size_t n=1;n < sz-1;n++)
+      data2[n] = data[n] + k*dt/(dx*dx) * (data[n+1] + data[n-1] - 2*data[n]);
+    data.swap(data2);
+
+    send_ghosts();
+  }
+
+  void send_ghosts() {
+    leftThread->right.push(data[1]);
+    rightThread->left.push(data[sz-2]);
+  }
+
+  void run() {
+    send_ghosts();
+    for(size_t t=0; t < nt; t++)
+      update();
+    recv_ghosts();
+  }
+};
+
+int main(int argc, char **argv) {
+  assert(argc == 4);
+  threads = std::stoi(argv[1]);
+  nt = std::stoi(argv[2]);
+  nx = std::stoi(argv[3]);
+  std::vector<Worker> workers;
+  size_t tx = (2*ghosts + nx)/threads;
+  for(size_t th=0;th < threads;th++) {
+    workers.emplace_back(Worker(th, tx));
+  }
+  for(size_t th=0;th < threads;th++) {
+    size_t next = (th + 1) % threads;
+    size_t prev = (threads + th - 1) % threads;
+    Worker& w = workers[th];
+    w.rightThread = &workers[next];
+    w.leftThread = &workers[prev];
+  }
+  auto t1 = std::chrono::high_resolution_clock::now();
+  for(Worker& w : workers) {
+    w.start();
+  }
+  for(Worker& w : workers) {
+    w.join();
+  }
+  auto t2 = std::chrono::high_resolution_clock::now();
+  double elapsed = std::chrono::duration_cast<std::chrono::nanoseconds>(t2 -t1).count()*1e-9;
+  std::cout << "elapsed: " << elapsed << std::endl;
+
+  return 0;
+}

hpx/heat/CMakeLists.txt

@@ -21,3 +21,4 @@ set(EXECUTABLE_OUTPUT_PATH "${PROJECT_BINARY_DIR}/bin")
 find_package(HPX REQUIRED)
 
 add_hpx_executable(heat_ghosts SOURCES heat_ghosts.cxx)
+add_hpx_executable(heat_ghosts_mini SOURCES heat_ghosts_mini.cxx)

hpx/heat/heat_ghosts_mini.cxx

@@ -0,0 +1,133 @@
+
+//  Copyright (c) 2023 AUTHORS
+//
+//  SPDX-License-Identifier: BSL-1.0
+//  Distributed under the Boost Software License, Version 1.0. (See accompanying
+//  file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#include <hpx/hpx.hpp>
+#include <hpx/hpx_main.hpp>
+#include <functional>
+#include <algorithm>
+#include <hpx/channel.hpp>
+#include <chrono>
+#include <filesystem>
+#include <fstream>
+
+using std::size_t;
+
+const size_t ghosts = 1;
+size_t nx = 20;//100000;
+const double k = 0.4;
+size_t nt = 1;
+const double dt = 1.0;
+const double dx = 1.0;
+size_t threads = 1;
+
+void pr(const std::vector<double>& total) {
+  std::cout << "[";
+  for(size_t i=0;i<total.size();i++) {
+    std::cout << " " << total[i];
+  }
+  std::cout << " ]" << std::endl;
+}
+
+class Worker : public hpx::thread {
+public:
+  size_t num;
+  size_t lo, hi, sz;
+  std::vector<double> data, data2;
+  static const size_t chsz = 20;
+  hpx::lcos::local::channel_spsc<double> left, right;
+  Worker *leftThread = nullptr;
+  Worker *rightThread = nullptr;
+
+  void start() {
+    hpx::thread t(std::bind(&Worker::run, this));
+    this->swap(t);
+  }
+
+  Worker(size_t num_, size_t tx) : num(num_), left(chsz), right(chsz) {
+    lo = tx * num;
+    hi = tx * (num + 1);
+    if (hi > nx) hi = nx;
+    lo -= ghosts;
+    hi += ghosts;
+    sz = hi - lo;
+    data.resize(sz);
+    data2.resize(sz);
+
+    size_t off=1;
+    for(size_t n=0; n < sz; n++) {
+      data.at(n) = n + lo + off;
+      data2.at(n) = 0.0;
+    }
+  }
+  Worker(Worker&&) noexcept = default;
+  ~Worker() {}
+
+  void recv_ghosts() {
+    //data[0] = left.get().get();
+    //data[sz-1] = right.get().get();
+    while(!left.get(&data[0]))
+      hpx::this_thread::yield();
+    while(!right.get(&data[sz-1]))
+      hpx::this_thread::yield();
+  }
+
+  void update() {
+    recv_ghosts();
+
+    for(size_t n=1;n < sz-1;n++)
+      data2[n] = data[n] + k*dt/(dx*dx) * (data[n+1] + data[n-1] - 2*data[n]);
+    data.swap(data2);
+
+    send_ghosts();
+  }
+
+  void send_ghosts() {
+    double d1 = data[1];
+    double d2 = data[sz-2];
+    leftThread->right.set(std::move(d1));
+    rightThread->left.set(std::move(d2));
+  }
+
+  void run() {
+    send_ghosts();
+    for(size_t t=0; t < nt; t++)
+      update();
+    recv_ghosts();
+  }
+};
+
+int main(int argc, char **argv) {
+  assert(argc == 4);
+  threads = std::stoi(argv[1]);
+  nt = std::stoi(argv[2]);
+  nx = std::stoi(argv[3]);
+  std::vector<Worker> workers;
+  size_t tx = (2*ghosts + nx)/threads;
+  for(size_t th=0;th < threads;th++) {
+    Worker w(th, tx);
+    workers.emplace_back(std::move(w));
+  }
+  for(size_t th=0;th < threads;th++) {
+    size_t next = (th + 1) % threads;
+    size_t prev = (threads + th - 1) % threads;
+    Worker& w = workers[th];
+    w.rightThread = &workers[next];
+    w.leftThread = &workers[prev];
+  }
+  auto t1 = std::chrono::high_resolution_clock::now();
+  for(Worker& w : workers) {
+    w.start();
+  }
+  for(Worker& w : workers) {
+    w.join();
+  }
+  auto t2 = std::chrono::high_resolution_clock::now();
+  double elapsed = std::chrono::duration_cast<std::chrono::nanoseconds>(t2 -t1).count()*1e-9;
+  std::cout << "elapsed: " << elapsed << std::endl;
+
+  return 0;
+}