Reorganize benchmarks into subdirectories and clean up Eigen sources

libeigen/eigen!2176

Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>

benchmarks/Core/CMakeLists.txt

@@ -0,0 +1,7 @@
+eigen_add_benchmark(bench_gemm bench_gemm.cpp)
+eigen_add_benchmark(bench_gemm_double bench_gemm.cpp DEFINITIONS SCALAR=double)
+eigen_add_benchmark(bench_gemv bench_gemv.cpp)
+eigen_add_benchmark(bench_vecadd bench_vecadd.cpp)
+eigen_add_benchmark(bench_trsm bench_trsm.cpp)
+eigen_add_benchmark(bench_reverse bench_reverse.cpp)
+eigen_add_benchmark(bench_move_semantics bench_move_semantics.cpp)

benchmarks/Core/bench_gemm.cpp

@@ -0,0 +1,81 @@
+#include <benchmark/benchmark.h>
+#include <Eigen/Core>
+
+using namespace Eigen;
+
+#ifndef SCALAR
+#define SCALAR float
+#endif
+
+typedef SCALAR Scalar;
+typedef Matrix<Scalar, Dynamic, Dynamic> Mat;
+
+template <typename A, typename B, typename C>
+EIGEN_DONT_INLINE void gemm(const A& a, const B& b, C& c) {
+  c.noalias() += a * b;
+}
+
+static void BM_EigenGemm(benchmark::State& state) {
+  int m = state.range(0);
+  int n = state.range(1);
+  int p = state.range(2);
+  Mat a(m, p);
+  a.setRandom();
+  Mat b(p, n);
+  b.setRandom();
+  Mat c = Mat::Zero(m, n);
+  for (auto _ : state) {
+    c.setZero();
+    gemm(a, b, c);
+    benchmark::DoNotOptimize(c.data());
+    benchmark::ClobberMemory();
+  }
+  state.counters["GFLOPS"] =
+      benchmark::Counter(2.0 * m * n * p, benchmark::Counter::kIsIterationInvariantRate, benchmark::Counter::kIs1000);
+}
+
+static void GemmSizes(::benchmark::Benchmark* b) {
+  for (int size : {8, 16, 32, 64, 96, 128, 160, 192, 224, 256, 288, 320, 384, 448, 512, 768, 1024, 1536, 2048}) {
+    b->Args({size, size, size});
+  }
+  // Non-square sizes
+  b->Args({64, 64, 1024});
+  b->Args({1024, 64, 64});
+  b->Args({64, 1024, 64});
+  b->Args({256, 256, 1024});
+  b->Args({1024, 256, 256});
+}
+
+BENCHMARK(BM_EigenGemm)->Apply(GemmSizes);
+
+#ifdef HAVE_BLAS
+extern "C" {
+#include <Eigen/src/misc/blas.h>
+}
+
+static void BM_BlasGemm(benchmark::State& state) {
+  int m = state.range(0);
+  int n = state.range(1);
+  int p = state.range(2);
+  Mat a(m, p);
+  a.setRandom();
+  Mat b(p, n);
+  b.setRandom();
+  Mat c = Mat::Zero(m, n);
+  char notrans = 'N';
+  Scalar one = 1, zero = 0;
+  for (auto _ : state) {
+    c.setZero();
+    if constexpr (std::is_same_v<Scalar, float>) {
+      sgemm_(&notrans, &notrans, &m, &n, &p, &one, a.data(), &m, b.data(), &p, &one, c.data(), &m);
+    } else {
+      dgemm_(&notrans, &notrans, &m, &n, &p, &one, a.data(), &m, b.data(), &p, &one, c.data(), &m);
+    }
+    benchmark::DoNotOptimize(c.data());
+    benchmark::ClobberMemory();
+  }
+  state.counters["GFLOPS"] =
+      benchmark::Counter(2.0 * m * n * p, benchmark::Counter::kIsIterationInvariantRate, benchmark::Counter::kIs1000);
+}
+BENCHMARK(BM_BlasGemm)->Apply(GemmSizes);
+#endif

benchmarks/Core/bench_gemv.cpp

@@ -0,0 +1,157 @@
+// Benchmark for dense general matrix-vector multiplication (GEMV).
+//
+// Tests performance of y += op(A) * x for various matrix sizes, aspect ratios,
+// scalar types, and operation variants (transpose, conjugate, adjoint).
+//
+// The Eigen GEMV kernel (Eigen/src/Core/products/GeneralMatrixVector.h) has
+// two main specializations:
+//   - ColMajor kernel: used for y += A * x with column-major A.
+//     Processes vertical panels, vectorizes along rows.
+//   - RowMajor kernel: used for y += A^T * x with column-major A.
+//     Processes groups of rows, vectorizes the dot product along columns.
+//
+// For complex scalars, conjugation flags (ConjugateLhs, ConjugateRhs) select
+// additional code paths within each kernel via conj_helper.
+//
+// Operation mapping (for column-major stored A):
+//   Gemv       y += A * x           -> ColMajor kernel, no conjugation
+//   GemvTrans  y += A^T * x         -> RowMajor kernel, no conjugation
+//   GemvConj   y += conj(A) * x     -> ColMajor kernel, ConjugateLhs=true
+//   GemvAdj    y += A^H * x         -> RowMajor kernel, ConjugateLhs=true
+
+#include <benchmark/benchmark.h>
+#include <Eigen/Core>
+
+using namespace Eigen;
+
+// ---------- Benchmark helpers ----------
+
+// GEMV flop count: 2*m*n for real, 8*m*n for complex.
+template <typename Scalar>
+double gemvFlops(Index m, Index n) {
+  return (NumTraits<Scalar>::IsComplex ? 8.0 : 2.0) * m * n;
+}
+
+// ---------- y += A * x  (ColMajor GEMV kernel, no conjugation) ----------
+
+template <typename Scalar>
+static void BM_Gemv(benchmark::State& state) {
+  using Mat = Matrix<Scalar, Dynamic, Dynamic>;
+  using Vec = Matrix<Scalar, Dynamic, 1>;
+  const Index m = state.range(0);
+  const Index n = state.range(1);
+  Mat A = Mat::Random(m, n);
+  Vec x = Vec::Random(n);
+  Vec y = Vec::Random(m);
+  for (auto _ : state) {
+    y.noalias() += A * x;
+    benchmark::DoNotOptimize(y.data());
+    benchmark::ClobberMemory();
+  }
+  state.counters["GFLOPS"] = benchmark::Counter(gemvFlops<Scalar>(m, n), benchmark::Counter::kIsIterationInvariantRate,
+                                                benchmark::Counter::kIs1000);
+}
+
+// ---------- y += A^T * x  (RowMajor GEMV kernel, no conjugation) ----------
+
+template <typename Scalar>
+static void BM_GemvTrans(benchmark::State& state) {
+  using Mat = Matrix<Scalar, Dynamic, Dynamic>;
+  using Vec = Matrix<Scalar, Dynamic, 1>;
+  const Index m = state.range(0);
+  const Index n = state.range(1);
+  Mat A = Mat::Random(m, n);
+  Vec x = Vec::Random(m);
+  Vec y = Vec::Random(n);
+  for (auto _ : state) {
+    y.noalias() += A.transpose() * x;
+    benchmark::DoNotOptimize(y.data());
+    benchmark::ClobberMemory();
+  }
+  state.counters["GFLOPS"] = benchmark::Counter(gemvFlops<Scalar>(m, n), benchmark::Counter::kIsIterationInvariantRate,
+                                                benchmark::Counter::kIs1000);
+}
+
+// ---------- y += conj(A) * x  (ColMajor kernel, ConjugateLhs=true) ----------
+
+template <typename Scalar>
+static void BM_GemvConj(benchmark::State& state) {
+  using Mat = Matrix<Scalar, Dynamic, Dynamic>;
+  using Vec = Matrix<Scalar, Dynamic, 1>;
+  const Index m = state.range(0);
+  const Index n = state.range(1);
+  Mat A = Mat::Random(m, n);
+  Vec x = Vec::Random(n);
+  Vec y = Vec::Random(m);
+  for (auto _ : state) {
+    y.noalias() += A.conjugate() * x;
+    benchmark::DoNotOptimize(y.data());
+    benchmark::ClobberMemory();
+  }
+  state.counters["GFLOPS"] = benchmark::Counter(gemvFlops<Scalar>(m, n), benchmark::Counter::kIsIterationInvariantRate,
+                                                benchmark::Counter::kIs1000);
+}
+
+// ---------- y += A^H * x  (RowMajor kernel, ConjugateLhs=true) ----------
+
+template <typename Scalar>
+static void BM_GemvAdj(benchmark::State& state) {
+  using Mat = Matrix<Scalar, Dynamic, Dynamic>;
+  using Vec = Matrix<Scalar, Dynamic, 1>;
+  const Index m = state.range(0);
+  const Index n = state.range(1);
+  Mat A = Mat::Random(m, n);
+  Vec x = Vec::Random(m);
+  Vec y = Vec::Random(n);
+  for (auto _ : state) {
+    y.noalias() += A.adjoint() * x;
+    benchmark::DoNotOptimize(y.data());
+    benchmark::ClobberMemory();
+  }
+  state.counters["GFLOPS"] = benchmark::Counter(gemvFlops<Scalar>(m, n), benchmark::Counter::kIsIterationInvariantRate,
+                                                benchmark::Counter::kIs1000);
+}
+
+// ---------- Size configurations ----------
+// All sizes refer to the stored matrix A (m rows, n cols).
+
+static void GemvSizes(::benchmark::Benchmark* b) {
+  // Square matrices: exercises balanced kernel behavior.
+  for (int size : {8, 32, 128, 512, 1024}) {
+    b->Args({size, size});
+  }
+  // Tall-thin (m >> n): in ColMajor kernel, the inner vectorized loop over rows
+  // is long while the outer column loop is short. In RowMajor kernel (transpose),
+  // there are many rows to process but short dot products.
+  for (int n : {1, 16}) {
+    for (int m : {256, 1024}) {
+      b->Args({m, n});
+    }
+  }
+  // Short-wide (m << n): in ColMajor kernel, the outer column loop is long but
+  // the inner vectorized loop over rows is short. In RowMajor kernel (transpose),
+  // there are few rows but long dot products.
+  for (int m : {1, 16}) {
+    for (int n : {256, 1024}) {
+      b->Args({m, n});
+    }
+  }
+}
+
+// ---------- Register benchmarks ----------
+
+// Real types: Gemv and GemvTrans exercise the two kernel specializations.
+// Conjugation is a no-op for real scalars.
+
+BENCHMARK(BM_Gemv<float>)->Apply(GemvSizes)->Name("Gemv_float");
+BENCHMARK(BM_Gemv<double>)->Apply(GemvSizes)->Name("Gemv_double");
+BENCHMARK(BM_GemvTrans<float>)->Apply(GemvSizes)->Name("GemvTrans_float");
+BENCHMARK(BM_GemvTrans<double>)->Apply(GemvSizes)->Name("GemvTrans_double");
+
+// Complex types: all four variants exercise distinct kernel code paths.
+// Only cfloat is benchmarked since cdouble exercises the same paths but slower.
+
+BENCHMARK(BM_Gemv<std::complex<float>>)->Apply(GemvSizes)->Name("Gemv_cfloat");
+BENCHMARK(BM_GemvTrans<std::complex<float>>)->Apply(GemvSizes)->Name("GemvTrans_cfloat");
+BENCHMARK(BM_GemvConj<std::complex<float>>)->Apply(GemvSizes)->Name("GemvConj_cfloat");
+BENCHMARK(BM_GemvAdj<std::complex<float>>)->Apply(GemvSizes)->Name("GemvAdj_cfloat");

benchmarks/Core/bench_move_semantics.cpp

@@ -0,0 +1,41 @@
+#include <benchmark/benchmark.h>
+#include <Eigen/Core>
+#include "../../test/MovableScalar.h"
+#include <utility>
+
+template <typename MatrixType>
+void copy_matrix(MatrixType& m) {
+  MatrixType tmp(m);
+  m = tmp;
+}
+
+template <typename MatrixType>
+void move_matrix(MatrixType&& m) {
+  MatrixType tmp(std::move(m));
+  m = std::move(tmp);
+}
+
+template <typename Scalar>
+static void BM_CopySemantics(benchmark::State& state) {
+  using MatrixType = Eigen::Matrix<Eigen::MovableScalar<Scalar>, 1, 10>;
+  MatrixType data = MatrixType::Random().eval();
+  for (auto _ : state) {
+    copy_matrix(data);
+    benchmark::DoNotOptimize(data.data());
+  }
+}
+
+template <typename Scalar>
+static void BM_MoveSemantics(benchmark::State& state) {
+  using MatrixType = Eigen::Matrix<Eigen::MovableScalar<Scalar>, 1, 10>;
+  MatrixType data = MatrixType::Random().eval();
+  for (auto _ : state) {
+    move_matrix(std::move(data));
+    benchmark::DoNotOptimize(data.data());
+  }
+}
+
+BENCHMARK(BM_CopySemantics<float>);
+BENCHMARK(BM_MoveSemantics<float>);
+BENCHMARK(BM_CopySemantics<double>);
+BENCHMARK(BM_MoveSemantics<double>);

benchmarks/Core/bench_reverse.cpp

@@ -0,0 +1,30 @@
+#include <benchmark/benchmark.h>
+#include <Eigen/Core>
+
+using namespace Eigen;
+
+static void BM_MatrixReverse(benchmark::State& state) {
+  int n = state.range(0);
+  typedef Matrix<double, Dynamic, Dynamic> MatrixType;
+  MatrixType a = MatrixType::Random(n, n);
+  MatrixType b(n, n);
+  for (auto _ : state) {
+    b = a.reverse();
+    benchmark::DoNotOptimize(b.data());
+  }
+  state.SetBytesProcessed(state.iterations() * n * n * sizeof(double));
+}
+BENCHMARK(BM_MatrixReverse)->RangeMultiplier(2)->Range(4, 512);
+
+static void BM_VectorReverse(benchmark::State& state) {
+  int n = state.range(0);
+  typedef Matrix<double, Dynamic, 1> VectorType;
+  VectorType a = VectorType::Random(n);
+  VectorType b(n);
+  for (auto _ : state) {
+    b = a.reverse();
+    benchmark::DoNotOptimize(b.data());
+  }
+  state.SetBytesProcessed(state.iterations() * n * sizeof(double));
+}
+BENCHMARK(BM_VectorReverse)->RangeMultiplier(4)->Range(16, 1 << 18);

benchmarks/Core/bench_trsm.cpp

@@ -0,0 +1,94 @@
+#include <benchmark/benchmark.h>
+#include <Eigen/Dense>
+
+using namespace Eigen;
+
+// ---------- TRSV: triangular solve with single RHS vector ----------
+
+template <typename Scalar, unsigned int Mode>
+static void BM_TRSV(benchmark::State& state) {
+  using Mat = Matrix<Scalar, Dynamic, Dynamic>;
+  using Vec = Matrix<Scalar, Dynamic, 1>;
+  const Index n = state.range(0);
+  Mat A = Mat::Random(n, n);
+  // Make diagonally dominant to ensure well-conditioned triangular part.
+  A.diagonal().array() += Scalar(n);
+  Vec x = Vec::Random(n);
+  Vec b = x;
+  for (auto _ : state) {
+    x = b;
+    A.template triangularView<Mode>().solveInPlace(x);
+    benchmark::DoNotOptimize(x.data());
+  }
+  state.SetItemsProcessed(state.iterations() * n * n);
+}
+
+// ---------- TRSM: triangular solve with multiple RHS (OnTheLeft) ----------
+
+template <typename Scalar, unsigned int Mode>
+static void BM_TRSM_Left(benchmark::State& state) {
+  using Mat = Matrix<Scalar, Dynamic, Dynamic>;
+  const Index n = state.range(0);
+  const Index nrhs = state.range(1);
+  Mat A = Mat::Random(n, n);
+  A.diagonal().array() += Scalar(n);
+  Mat X = Mat::Random(n, nrhs);
+  Mat B = X;
+  for (auto _ : state) {
+    X = B;
+    A.template triangularView<Mode>().solveInPlace(X);
+    benchmark::DoNotOptimize(X.data());
+  }
+  state.SetItemsProcessed(state.iterations() * n * n * nrhs);
+}
+
+// ---------- TRSM: triangular solve with multiple RHS (OnTheRight) ----------
+
+template <typename Scalar, unsigned int Mode>
+static void BM_TRSM_Right(benchmark::State& state) {
+  using Mat = Matrix<Scalar, Dynamic, Dynamic>;
+  const Index n = state.range(0);
+  const Index nrhs = state.range(1);
+  Mat A = Mat::Random(n, n);
+  A.diagonal().array() += Scalar(n);
+  Mat X = Mat::Random(nrhs, n);
+  Mat B = X;
+  for (auto _ : state) {
+    X = B;
+    A.template triangularView<Mode>().template solveInPlace<OnTheRight>(X);
+    benchmark::DoNotOptimize(X.data());
+  }
+  state.SetItemsProcessed(state.iterations() * n * n * nrhs);
+}
+
+// ---------- Size configurations ----------
+
+static void TrsvSizes(::benchmark::Benchmark* b) {
+  for (int n : {32, 128, 512}) {
+    b->Args({n});
+  }
+}
+
+static void TrsmSizes(::benchmark::Benchmark* b) {
+  for (int n : {64, 256, 512}) {
+    for (int nrhs : {1, 16, 64}) {
+      b->Args({n, nrhs});
+    }
+  }
+}
+
+// ---------- TRSV benchmarks ----------
+// Only Lower is benchmarked; Upper exercises the same kernel via transposed storage.
+
+BENCHMARK(BM_TRSV<float, Lower>)->Apply(TrsvSizes)->Name("TRSV_float_Lower");
+BENCHMARK(BM_TRSV<double, Lower>)->Apply(TrsvSizes)->Name("TRSV_double_Lower");
+
+// ---------- TRSM Left benchmarks ----------
+
+BENCHMARK(BM_TRSM_Left<float, Lower>)->Apply(TrsmSizes)->Name("TRSM_Left_float_Lower");
+BENCHMARK(BM_TRSM_Left<double, Lower>)->Apply(TrsmSizes)->Name("TRSM_Left_double_Lower");
+
+// ---------- TRSM Right benchmarks ----------
+
+BENCHMARK(BM_TRSM_Right<float, Lower>)->Apply(TrsmSizes)->Name("TRSM_Right_float_Lower");
+BENCHMARK(BM_TRSM_Right<double, Lower>)->Apply(TrsmSizes)->Name("TRSM_Right_double_Lower");

benchmarks/Core/bench_vecadd.cpp

@@ -0,0 +1,28 @@
+#include <benchmark/benchmark.h>
+#include <Eigen/Core>
+
+using namespace Eigen;
+
+static void BM_VecAdd(benchmark::State& state) {
+  int size = state.range(0);
+  VectorXf a = VectorXf::Random(size);
+  VectorXf b = VectorXf::Random(size);
+  for (auto _ : state) {
+    a = a + b;
+    benchmark::DoNotOptimize(a.data());
+  }
+  state.SetBytesProcessed(state.iterations() * size * sizeof(float) * 3);
+}
+BENCHMARK(BM_VecAdd)->RangeMultiplier(4)->Range(64, 1 << 20);
+
+static void BM_MatAdd(benchmark::State& state) {
+  int n = state.range(0);
+  MatrixXf a = MatrixXf::Random(n, n);
+  MatrixXf b = MatrixXf::Random(n, n);
+  for (auto _ : state) {
+    a = a + b;
+    benchmark::DoNotOptimize(a.data());
+  }
+  state.SetBytesProcessed(state.iterations() * n * n * sizeof(float) * 3);
+}
+BENCHMARK(BM_MatAdd)->RangeMultiplier(2)->Range(8, 512);