contrib/benchmarks/StructuredMatrices/bench_structured_hankel.cpp

// Benchmarks for the fast Hankel matrix-vector product: the structured operator
// (O(n log n) via a cached DFT symbol, with direct O(n) dispatch for single-row
// and single-column shapes) against the equivalent dense GEMV. The *Dense
// variants evaluate the same product through a materialized dense matrix.
// SPDX-FileCopyrightText: The Eigen Authors
// SPDX-License-Identifier: MPL-2.0

#include <benchmark/benchmark.h>
#include <Eigen/Core>
#include <contrib/Eigen/StructuredMatrices>

using namespace Eigen;

typedef Matrix<double, Dynamic, 1> Vec;
typedef Matrix<double, Dynamic, Dynamic> Mat;

// --- Square n x n product: FFT-based operator vs dense GEMV ---
static void BM_HankelProduct(benchmark::State& state) {
  const Index n = state.range(0);
  Vec h = Vec::Random(2 * n - 1);
  Hankel<double> H(h.head(n), h.tail(n));  // the FFT symbol is paid once, here
  Vec x = Vec::Random(n), y(n);
  for (auto _ : state) {
    y.noalias() = H * x;
    benchmark::DoNotOptimize(y.data());
  }
}
BENCHMARK(BM_HankelProduct)->Arg(256)->Arg(1024)->Arg(4096)->Arg(16384);

static void BM_HankelProductDense(benchmark::State& state) {
  const Index n = state.range(0);
  Vec h = Vec::Random(2 * n - 1);
  Mat dense = Hankel<double>(h.head(n), h.tail(n));
  Vec x = Vec::Random(n), y(n);
  for (auto _ : state) {
    y.noalias() = dense * x;
    benchmark::DoNotOptimize(y.data());
  }
}
BENCHMARK(BM_HankelProductDense)->Arg(256)->Arg(1024)->Arg(4096)->Arg(16384);

// --- Skinny 1 x n: dispatched to a direct O(n) dot-product kernel, not the FFT
// (whose O(n log n) cost plus setup is asymptotically worse than direct here) ---
static void BM_HankelProductSkinnyRow(benchmark::State& state) {
  const Index n = state.range(0);
  Vec h = Vec::Random(n);  // generating sequence of the 1 x n operator
  Hankel<double> H(h.head(1), h.tail(n));
  Vec x = Vec::Random(n);
  Matrix<double, 1, 1> y;
  for (auto _ : state) {
    y.noalias() = H * x;
    benchmark::DoNotOptimize(y.data());
  }
}
BENCHMARK(BM_HankelProductSkinnyRow)->Arg(4096)->Arg(65536)->Arg(1048576);

static void BM_HankelProductSkinnyRowDense(benchmark::State& state) {
  const Index n = state.range(0);
  Vec h = Vec::Random(n);
  Mat dense = Hankel<double>(h.head(1), h.tail(n));
  Vec x = Vec::Random(n);
  Matrix<double, 1, 1> y;
  for (auto _ : state) {
    y.noalias() = dense * x;
    benchmark::DoNotOptimize(y.data());
  }
}
BENCHMARK(BM_HankelProductSkinnyRowDense)->Arg(4096)->Arg(65536)->Arg(1048576);

// --- Skinny n x 1: a single axpy of the generating sequence ---
static void BM_HankelProductSkinnyCol(benchmark::State& state) {
  const Index n = state.range(0);
  Vec h = Vec::Random(n);  // generating sequence of the n x 1 operator
  Hankel<double> H(h.head(n), h.tail(1));
  Matrix<double, 1, 1> x = Matrix<double, 1, 1>::Random();
  Vec y(n);
  for (auto _ : state) {
    y.noalias() = H * x;
    benchmark::DoNotOptimize(y.data());
  }
}
BENCHMARK(BM_HankelProductSkinnyCol)->Arg(4096)->Arg(65536)->Arg(1048576);

static void BM_HankelProductSkinnyColDense(benchmark::State& state) {
  const Index n = state.range(0);
  Vec h = Vec::Random(n);
  Mat dense = Hankel<double>(h.head(n), h.tail(1));
  Matrix<double, 1, 1> x = Matrix<double, 1, 1>::Random();
  Vec y(n);
  for (auto _ : state) {
    y.noalias() = dense * x;
    benchmark::DoNotOptimize(y.data());
  }
}
BENCHMARK(BM_HankelProductSkinnyColDense)->Arg(4096)->Arg(65536)->Arg(1048576);

template <typename Scalar>
static void BM_HankelTranspose(benchmark::State& state) {
  using Real = typename NumTraits<Scalar>::Real;
  using Vector = Matrix<Scalar, Dynamic, 1>;
  const Index m = state.range(0), n = state.range(1);
  const Vector h = Vector::Random(m + n - 1);
  const Hankel<Scalar> op(h.head(m), h.tail(n));
  const auto symbol = op.symbol();
  const auto roundTrip = op.transpose().transpose().symbol();
  const Real error = (roundTrip - symbol).norm() / symbol.norm();
  if (!(error <= Real(16) * NumTraits<Real>::epsilon())) {
    state.SkipWithError("incorrect transpose symbol");
    return;
  }
  state.counters["relative_error"] = double(error);
  for (auto _ : state) {
    auto transposed = op.transpose();
    benchmark::DoNotOptimize(transposed);
  }
}
BENCHMARK_TEMPLATE(BM_HankelTranspose, float)->Args({48, 33})->Args({256, 512})->Args({2048, 1024});
BENCHMARK_TEMPLATE(BM_HankelTranspose, double)->Args({48, 33})->Args({256, 512})->Args({2048, 1024});
BENCHMARK_TEMPLATE(BM_HankelTranspose, std::complex<double>)->Args({48, 33})->Args({256, 512})->Args({2048, 1024});