Extract filtering to Filter class (TODO tests). Add cached versions and benchmark.

Author: ProExpertProg

cpp/CMakeLists.txt

@@ -61,6 +61,7 @@ target_link_libraries(adaptive_hermite_refinement argparse mdspan)
 add_library(src-lib OBJECT lib/Naive.cpp lib/Naive.hpp lib/HermiteRunner.cpp lib/HermiteRunner.hpp
         lib/Transformer.hpp lib/Transformer.cpp
         lib/Exporter.hpp lib/Exporter.cpp
+        lib/Filter.hpp lib/Filter.cpp
 )
 target_include_directories(src-lib PUBLIC lib/)
 target_link_libraries(src-lib mdspan fftw-cpp cnpy spdlog::spdlog)

cpp/bench/CMakeLists.txt

@@ -22,6 +22,9 @@ FetchContent_MakeAvailable(google-benchmark)
 add_executable(bench-naive naive.cpp main-no-log.cpp)
 target_link_libraries(bench-naive src-lib benchmark::benchmark)
 
+add_executable(bench-hl-filter hl-filter.cpp)
+target_link_libraries(bench-hl-filter src-lib benchmark::benchmark_main)
+
 if (ENABLE_CILK)
     add_executable(bench-fftw-cilk fftw-cilk.cpp)
     target_link_libraries(bench-fftw-cilk fftw-cpp benchmark::benchmark fftw3_threads)

cpp/bench/hl-filter.cpp

@@ -0,0 +1,37 @@
+#include "Filter.hpp"
+
+#include <benchmark/benchmark.h>
+
+static constexpr auto N_WARMUP_ITERS = 5;
+
+using namespace ahr;
+template <class Filter> static void BM_HouLiFilter(benchmark::State &state) {
+  Dim const X = state.range(0);
+  Dim const Y = state.range(1);
+
+  Grid grid{1, X, Y};
+  Filter filter{grid};
+
+  auto buf = grid.cBufXY();
+
+  // Warm-up
+  for (int i = 0; i < N_WARMUP_ITERS; i++) {
+    filter(buf);
+  }
+
+  for (auto _ : state) {
+    filter(buf);
+  }
+}
+
+BENCHMARK(BM_HouLiFilter<HouLiFilter>)
+    ->ArgsProduct({{2048, 4096, 8192}, {2048, 4096, 8192}})
+    ->Unit(benchmark::kMillisecond);
+
+BENCHMARK(BM_HouLiFilter<HouLiFilterCached>)
+    ->ArgsProduct({{2048, 4096, 8192}, {2048, 4096, 8192}})
+    ->Unit(benchmark::kMillisecond);
+
+BENCHMARK(BM_HouLiFilter<HouLiFilterCached1D>)
+    ->ArgsProduct({{2048, 4096, 8192}, {2048, 4096, 8192}})
+    ->Unit(benchmark::kMillisecond);

cpp/bench/naive.cpp

@@ -5,8 +5,6 @@
 
 using namespace ahr;
 static void BM_Naive(benchmark::State &state) {
-  std::ostringstream oss;
-
   Dim const M = state.range(0);
   Dim const X = state.range(1);
   Dim const N = state.range(2);

cpp/lib/Filter.cpp

@@ -0,0 +1,42 @@
+//
+// Created by Luka on 11/2/2024.
+//
+
+#include "Filter.hpp"
+namespace ahr {
+void HouLiFilter::operator()(Grid::View::C_XY view) {
+  grid.for_each_kxky([&](Dim kx, Dim ky) {
+    view(kx, ky) *=
+        exp(-36.0 * pow(kx_(kx) / grid.KX, 36.0)) * exp(-36.0 * pow(ky_(ky) / grid.KY, 36.0));
+  });
+}
+
+HouLiFilterCached::HouLiFilterCached(Grid const &grid)
+    : HouLiFilter(grid), factors(std::array{grid.KX, grid.KY}) {
+  grid.for_each_kxky([&](Dim kx, Dim ky) {
+    factors(kx, ky) =
+        exp(-36.0 * pow(kx_(kx) / grid.KX, 36.0)) * exp(-36.0 * pow(ky_(ky) / grid.KY, 36.0));
+  });
+}
+
+void HouLiFilterCached::operator()(Grid::View::C_XY view) {
+  grid.for_each_kxky([&](Dim kx, Dim ky) { view(kx, ky) *= factors(kx, ky); });
+}
+
+HouLiFilterCached1D::HouLiFilterCached1D(Grid const &grid)
+    : HouLiFilter(grid), factors_x(grid.KX), factors_y(grid.KY) {
+  for (Dim kx = 0; kx < grid.KX; ++kx) {
+    factors_x.at(kx) = exp(-36.0 * pow(kx_(kx) / grid.KX, 36.0));
+  }
+  for (Dim ky = 0; ky < grid.KY; ++ky) {
+    factors_y.at(ky) = exp(-36.0 * pow(ky_(ky) / grid.KY, 36.0));
+  }
+}
+
+void HouLiFilterCached1D::operator()(Grid::View::C_XY view) {
+  grid.for_each_kxky([&](Dim kx, Dim ky) {
+    // Extra multiplication at runtime for lower memory cost
+    view(kx, ky) *= factors_x[kx] * factors_y[ky];
+  });
+}
+} // namespace ahr

cpp/lib/Filter.hpp

@@ -0,0 +1,42 @@
+#pragma once
+#include "constants.hpp"
+#include "grid.hpp"
+
+namespace ahr {
+
+class HouLiFilter {
+public:
+  explicit HouLiFilter(Grid const &grid) : grid(grid) {}
+
+  void operator()(Grid::View::C_XY view);
+
+protected:
+  Grid const &grid;
+
+  // TODO extract to common utility
+  [[nodiscard]] Real ky_(Dim ky) const {
+    return (ky <= (grid.KY / 2) ? Real(ky) : Real(ky) - Real(grid.KY)) * Real(lx) / Real(ly);
+  }
+  [[nodiscard]] Real kx_(Dim kx) const { return Real(kx); }
+};
+
+class HouLiFilterCached : HouLiFilter {
+public:
+  explicit HouLiFilterCached(Grid const &grid);
+  void operator()(Grid::View::C_XY view);
+
+private:
+  /// Pre-calculated factors for the Hou-Li filter.
+  /// Note that this is a real buffer with dimensions (KX,KY)
+  Grid::Buf::R_XY factors;
+};
+
+class HouLiFilterCached1D : protected HouLiFilter {
+public:
+  explicit HouLiFilterCached1D(Grid const &grid);
+  void operator()(Grid::View::C_XY view);
+
+protected:
+  std::vector<Real> factors_x, factors_y;
+};
+} // namespace ahr

cpp/lib/Naive.cpp

@@ -19,13 +19,6 @@ Naive::Naive(Dim M, Dim X, Dim Y)
   assert((Y & (Y - 1)) == 0);
 }
 
-void Naive::hlFilter(View::C_XY &complexArray) {
-  g.for_each_kxky([&](Dim kx, Dim ky) {
-    complexArray(kx, ky) *=
-        exp(-36.0 * pow(kx_(kx) / g.KX, 36.0)) * exp(-36.0 * pow(ky_(ky) / g.KY, 36.0));
-  });
-}
-
 void Naive::fftHL(View::R_XY in, View::C_XY out) {
   tf.fft(in, out);
   hlFilter(out);

cpp/lib/Naive.hpp

@@ -1,6 +1,7 @@
 #pragma once
 
 #include "Exporter.hpp"
+#include "Filter.hpp"
 #include "HermiteRunner.hpp"
 #include "Transformer.hpp"
 #include "constants.hpp"
@@ -35,6 +36,8 @@ class Naive : public ahr::HermiteRunner {
   Grid g;
   Transformer tf{g};
   Exporter exporter{g, tf};
+  HouLiFilter hlFilter{g};
+
 private:
 
   using View = Grid::View;
@@ -45,7 +48,6 @@ class Naive : public ahr::HermiteRunner {
   Real dt{-1};        ///< timestep
   Real elapsedT{0.0}; ///< total time elapsed
 
-  void hlFilter(View::C_XY &complexArray);
   void fftHL(View::R_XY in, View::C_XY out); ///< FFT with Hou-Li Filter
 
   Real bPerpMax{0};
@@ -127,8 +129,8 @@ class Naive : public ahr::HermiteRunner {
 
   [[nodiscard]] Real ky_(Dim ky) const {
     return (ky <= (g.KY / 2) ? Real(ky) : Real(ky) - Real(g.KY)) * Real(lx) / Real(ly);
-  };
-  [[nodiscard]] Real kx_(Dim kx) const { return Real(kx); };
+  }
+  [[nodiscard]] Real kx_(Dim kx) const { return Real(kx); }
 
   [[nodiscard]] Real kPerp2(Dim kx, Dim ky) const {
     auto dkx = kx_(kx), dky = ky_(ky);