673 lines
25 KiB
C++
673 lines
25 KiB
C++
// SPDX-FileCopyrightText: The Eigen Authors
|
||
// SPDX-License-Identifier: MPL-2.0
|
||
|
||
// Benchmark: SparseMatrix vs BlockSparseMatrix, real and complex scalars.
|
||
//
|
||
// Naming convention:
|
||
// BM_Sm_* — SparseMatrix (standard Eigen sparse)
|
||
// BM_BSM_* — BlockSparseMatrix
|
||
// *_SpMV — sparse × dense vector (GEMV)
|
||
// *_TriMV — triangular-view × dense vector
|
||
// *_SymmMV — selfadjoint-view × dense vector
|
||
// *_TriSolve — triangular solve in-place
|
||
// BM_Sm_Sm_* — SparseMatrix × SparseMatrix
|
||
// BM_BSM_BSM_* — BlockSparseMatrix × BlockSparseMatrix
|
||
// DiagT — DiagIsTriangular=true (diagonal blocks are actually triangular)
|
||
// DiagNSA — DiagIsSelfAdjoint=false with Hermitian diagonal blocks
|
||
// DiagSA — DiagIsSelfAdjoint=true with Hermitian diagonal blocks
|
||
|
||
#include <benchmark/benchmark.h>
|
||
#include <Eigen/Sparse>
|
||
|
||
#include <complex>
|
||
#include <random>
|
||
#include <set>
|
||
|
||
using namespace Eigen;
|
||
using cd = std::complex<double>;
|
||
using cf = std::complex<float>;
|
||
|
||
// ---------------------------------------------------------------------------
|
||
// Decode benchmark args: range(0)=nB, range(1)=sparsity% → nnzPerCol.
|
||
// ---------------------------------------------------------------------------
|
||
static void parseArgs(const benchmark::State& state, int& nB, int& nnz) {
|
||
nB = state.range(0);
|
||
nnz = std::max(1, (int)(state.range(1) * nB / 100));
|
||
}
|
||
|
||
// ---------------------------------------------------------------------------
|
||
// Random-value helper — works for float, double, complex<float>, complex<double>.
|
||
// ---------------------------------------------------------------------------
|
||
template <typename Scalar>
|
||
static typename std::enable_if<std::is_floating_point<Scalar>::value, Scalar>::type randVal(
|
||
std::mt19937& rng, std::normal_distribution<double>& d) {
|
||
return Scalar(d(rng));
|
||
}
|
||
|
||
template <typename Scalar>
|
||
static typename std::enable_if<!std::is_floating_point<Scalar>::value, Scalar>::type randVal(
|
||
std::mt19937& rng, std::normal_distribution<double>& d) {
|
||
using R = typename Scalar::value_type;
|
||
return Scalar{R(d(rng)), R(d(rng))};
|
||
}
|
||
|
||
// ---------------------------------------------------------------------------
|
||
// Build a general (full) block-sparse pair.
|
||
// ---------------------------------------------------------------------------
|
||
template <typename Scalar, int B>
|
||
static void buildPair(int nB, int nnzPerCol, unsigned seed, BlockSparseMatrix<Scalar, ColMajor, B, B>& bsm,
|
||
SparseMatrix<Scalar>& sm) {
|
||
using BSM = BlockSparseMatrix<Scalar, ColMajor, B, B>;
|
||
using BT = typename BSM::BlockType;
|
||
using Triplet = typename BSM::TripletType;
|
||
|
||
std::mt19937 rng(seed);
|
||
std::uniform_int_distribution<int> rowDist(0, nB - 1);
|
||
std::normal_distribution<double> vd;
|
||
|
||
std::vector<Triplet> triplets;
|
||
triplets.reserve(nB * nnzPerCol);
|
||
for (int j = 0; j < nB; ++j) {
|
||
std::set<int> rows;
|
||
rows.insert(j % nB);
|
||
while ((int)rows.size() < std::min(nnzPerCol, nB)) rows.insert(rowDist(rng));
|
||
for (int bi : rows) {
|
||
BT blk;
|
||
for (int r = 0; r < B; ++r)
|
||
for (int c = 0; c < B; ++c) blk(r, c) = randVal<Scalar>(rng, vd);
|
||
triplets.emplace_back(bi, j, blk);
|
||
}
|
||
}
|
||
bsm = BSM(nB, nB);
|
||
bsm.setFromTriplets(triplets.begin(), triplets.end());
|
||
sm = bsm.toSparse();
|
||
}
|
||
|
||
// ---------------------------------------------------------------------------
|
||
// Build an upper-triangular block-sparse pair (diagonal blocks are general).
|
||
// forSolve=true → diagonal blocks diagonally dominant (well-conditioned).
|
||
// ---------------------------------------------------------------------------
|
||
template <typename Scalar, int B>
|
||
static void buildUpperTriPair(int nB, int nnzPerCol, unsigned seed, BlockSparseMatrix<Scalar, ColMajor, B, B>& bsm,
|
||
SparseMatrix<Scalar>& sm, bool forSolve = false) {
|
||
using BSM = BlockSparseMatrix<Scalar, ColMajor, B, B>;
|
||
using BT = typename BSM::BlockType;
|
||
using Triplet = typename BSM::TripletType;
|
||
|
||
std::mt19937 rng(seed);
|
||
std::uniform_int_distribution<int> rowDist(0, nB - 1);
|
||
std::normal_distribution<double> vd;
|
||
|
||
std::vector<Triplet> triplets;
|
||
triplets.reserve(nB * nnzPerCol);
|
||
for (int j = 0; j < nB; ++j) {
|
||
BT diag;
|
||
for (int r = 0; r < B; ++r)
|
||
for (int c = 0; c < B; ++c) diag(r, c) = randVal<Scalar>(rng, vd);
|
||
if (forSolve) {
|
||
diag *= Scalar(0.1);
|
||
for (int d = 0; d < B; ++d) diag(d, d) += Scalar(double(B));
|
||
}
|
||
triplets.emplace_back(j, j, diag);
|
||
|
||
std::set<int> rows;
|
||
while ((int)rows.size() < std::min(nnzPerCol - 1, j)) rows.insert(rowDist(rng) % j);
|
||
for (int bi : rows) {
|
||
BT blk;
|
||
for (int r = 0; r < B; ++r)
|
||
for (int c = 0; c < B; ++c) blk(r, c) = randVal<Scalar>(rng, vd);
|
||
triplets.emplace_back(bi, j, blk);
|
||
}
|
||
}
|
||
bsm = BSM(nB, nB);
|
||
bsm.setFromTriplets(triplets.begin(), triplets.end());
|
||
sm = bsm.toSparse();
|
||
}
|
||
|
||
// ---------------------------------------------------------------------------
|
||
// Build an upper-triangular pair with actually-triangular diagonal blocks
|
||
// (strict lower triangle zeroed). Valid for DiagIsTriangular=true.
|
||
// ---------------------------------------------------------------------------
|
||
template <typename Scalar, int B>
|
||
static void buildActuallyTriPair(int nB, int nnzPerCol, unsigned seed, BlockSparseMatrix<Scalar, ColMajor, B, B>& bsm,
|
||
SparseMatrix<Scalar>& sm) {
|
||
using BSM = BlockSparseMatrix<Scalar, ColMajor, B, B>;
|
||
using BT = typename BSM::BlockType;
|
||
using Triplet = typename BSM::TripletType;
|
||
|
||
std::mt19937 rng(seed);
|
||
std::uniform_int_distribution<int> rowDist(0, nB - 1);
|
||
std::normal_distribution<double> vd;
|
||
|
||
std::vector<Triplet> triplets;
|
||
triplets.reserve(nB * nnzPerCol);
|
||
for (int j = 0; j < nB; ++j) {
|
||
BT diag = BT::Zero();
|
||
for (int r = 0; r < B; ++r)
|
||
for (int c = r; c < B; ++c) // upper triangle only
|
||
diag(r, c) = randVal<Scalar>(rng, vd);
|
||
triplets.emplace_back(j, j, diag);
|
||
|
||
std::set<int> rows;
|
||
while ((int)rows.size() < std::min(nnzPerCol - 1, j)) rows.insert(rowDist(rng) % j);
|
||
for (int bi : rows) {
|
||
BT blk;
|
||
for (int r = 0; r < B; ++r)
|
||
for (int c = 0; c < B; ++c) blk(r, c) = randVal<Scalar>(rng, vd);
|
||
triplets.emplace_back(bi, j, blk);
|
||
}
|
||
}
|
||
bsm = BSM(nB, nB);
|
||
bsm.setFromTriplets(triplets.begin(), triplets.end());
|
||
sm = bsm.toSparse();
|
||
}
|
||
|
||
// ---------------------------------------------------------------------------
|
||
// Build an upper-triangular pair with Hermitian diagonal blocks.
|
||
// Valid for DiagIsSelfAdjoint=true.
|
||
// ---------------------------------------------------------------------------
|
||
template <typename Scalar, int B>
|
||
static void buildHermDiagUpperTriPair(int nB, int nnzPerCol, unsigned seed,
|
||
BlockSparseMatrix<Scalar, ColMajor, B, B>& bsm, SparseMatrix<Scalar>& sm) {
|
||
using BSM = BlockSparseMatrix<Scalar, ColMajor, B, B>;
|
||
using BT = typename BSM::BlockType;
|
||
using Triplet = typename BSM::TripletType;
|
||
|
||
std::mt19937 rng(seed);
|
||
std::uniform_int_distribution<int> rowDist(0, nB - 1);
|
||
std::normal_distribution<double> vd;
|
||
|
||
std::vector<Triplet> triplets;
|
||
triplets.reserve(nB * nnzPerCol);
|
||
for (int j = 0; j < nB; ++j) {
|
||
BT raw;
|
||
for (int r = 0; r < B; ++r)
|
||
for (int c = 0; c < B; ++c) raw(r, c) = randVal<Scalar>(rng, vd);
|
||
BT diag = (raw + raw.adjoint()) / Scalar(2);
|
||
triplets.emplace_back(j, j, diag);
|
||
|
||
std::set<int> rows;
|
||
while ((int)rows.size() < std::min(nnzPerCol - 1, j)) rows.insert(rowDist(rng) % j);
|
||
for (int bi : rows) {
|
||
BT blk;
|
||
for (int r = 0; r < B; ++r)
|
||
for (int c = 0; c < B; ++c) blk(r, c) = randVal<Scalar>(rng, vd);
|
||
triplets.emplace_back(bi, j, blk);
|
||
}
|
||
}
|
||
bsm = BSM(nB, nB);
|
||
bsm.setFromTriplets(triplets.begin(), triplets.end());
|
||
sm = bsm.toSparse();
|
||
}
|
||
|
||
// ---------------------------------------------------------------------------
|
||
// Sparse×Sparse: Addition
|
||
// ---------------------------------------------------------------------------
|
||
template <typename Scalar, int B>
|
||
static void BM_Sm_Sm_Add(benchmark::State& state) {
|
||
int nB, nnz;
|
||
parseArgs(state, nB, nnz);
|
||
SparseMatrix<Scalar> smA, smB, smC;
|
||
BlockSparseMatrix<Scalar, ColMajor, B, B> tmp;
|
||
buildPair<Scalar, B>(nB, nnz, 1, tmp, smA);
|
||
buildPair<Scalar, B>(nB, nnz, 2, tmp, smB);
|
||
for (auto _ : state) {
|
||
smC = smA + smB;
|
||
benchmark::DoNotOptimize(smC.valuePtr());
|
||
}
|
||
state.counters["n"] = smA.rows();
|
||
}
|
||
|
||
template <typename Scalar, int B>
|
||
static void BM_BSM_BSM_Add(benchmark::State& state) {
|
||
int nB, nnz;
|
||
parseArgs(state, nB, nnz);
|
||
using BSM = BlockSparseMatrix<Scalar, ColMajor, B, B>;
|
||
BSM bsmA, bsmB, bsmC;
|
||
SparseMatrix<Scalar> smTmp;
|
||
buildPair<Scalar, B>(nB, nnz, 1, bsmA, smTmp);
|
||
buildPair<Scalar, B>(nB, nnz, 2, bsmB, smTmp);
|
||
for (auto _ : state) {
|
||
bsmC = bsmA + bsmB;
|
||
benchmark::DoNotOptimize(bsmC.valuePtr());
|
||
}
|
||
state.counters["n"] = bsmA.rows();
|
||
}
|
||
|
||
// ---------------------------------------------------------------------------
|
||
// Sparse×Dense: GEMV
|
||
// ---------------------------------------------------------------------------
|
||
template <typename Scalar, int B>
|
||
static void BM_Sm_SpMV(benchmark::State& state) {
|
||
int nB, nnz;
|
||
parseArgs(state, nB, nnz);
|
||
SparseMatrix<Scalar> sm;
|
||
BlockSparseMatrix<Scalar, ColMajor, B, B> tmp;
|
||
buildPair<Scalar, B>(nB, nnz, 1, tmp, sm);
|
||
Matrix<Scalar, Dynamic, 1> x = Matrix<Scalar, Dynamic, 1>::Random(sm.cols());
|
||
Matrix<Scalar, Dynamic, 1> y(sm.rows());
|
||
for (auto _ : state) {
|
||
y.noalias() = sm * x;
|
||
benchmark::DoNotOptimize(y.data());
|
||
}
|
||
state.counters["n"] = sm.rows();
|
||
}
|
||
|
||
template <typename Scalar, int B>
|
||
static void BM_BSM_SpMV(benchmark::State& state) {
|
||
int nB, nnz;
|
||
parseArgs(state, nB, nnz);
|
||
using BSM = BlockSparseMatrix<Scalar, ColMajor, B, B>;
|
||
BSM bsm;
|
||
SparseMatrix<Scalar> smTmp;
|
||
buildPair<Scalar, B>(nB, nnz, 1, bsm, smTmp);
|
||
Matrix<Scalar, Dynamic, 1> x = Matrix<Scalar, Dynamic, 1>::Random(bsm.cols());
|
||
Matrix<Scalar, Dynamic, 1> y(bsm.rows());
|
||
for (auto _ : state) {
|
||
y.noalias() = bsm * x;
|
||
benchmark::DoNotOptimize(y.data());
|
||
}
|
||
state.counters["n"] = bsm.rows();
|
||
}
|
||
|
||
// Accumulating GEMV: y += A * x.
|
||
template <typename Scalar, int B>
|
||
static void BM_Sm_SpMV_Acc(benchmark::State& state) {
|
||
int nB, nnz;
|
||
parseArgs(state, nB, nnz);
|
||
SparseMatrix<Scalar> sm;
|
||
BlockSparseMatrix<Scalar, ColMajor, B, B> tmp;
|
||
buildPair<Scalar, B>(nB, nnz, 1, tmp, sm);
|
||
Matrix<Scalar, Dynamic, 1> x = Matrix<Scalar, Dynamic, 1>::Random(sm.cols());
|
||
Matrix<Scalar, Dynamic, 1> y = Matrix<Scalar, Dynamic, 1>::Random(sm.rows());
|
||
for (auto _ : state) {
|
||
y.noalias() += sm * x;
|
||
benchmark::DoNotOptimize(y.data());
|
||
}
|
||
state.counters["n"] = sm.rows();
|
||
}
|
||
|
||
template <typename Scalar, int B>
|
||
static void BM_BSM_SpMV_Acc(benchmark::State& state) {
|
||
int nB, nnz;
|
||
parseArgs(state, nB, nnz);
|
||
using BSM = BlockSparseMatrix<Scalar, ColMajor, B, B>;
|
||
BSM bsm;
|
||
SparseMatrix<Scalar> smTmp;
|
||
buildPair<Scalar, B>(nB, nnz, 1, bsm, smTmp);
|
||
Matrix<Scalar, Dynamic, 1> x = Matrix<Scalar, Dynamic, 1>::Random(bsm.cols());
|
||
Matrix<Scalar, Dynamic, 1> y = Matrix<Scalar, Dynamic, 1>::Random(bsm.rows());
|
||
for (auto _ : state) {
|
||
y.noalias() += bsm * x;
|
||
benchmark::DoNotOptimize(y.data());
|
||
}
|
||
state.counters["n"] = bsm.rows();
|
||
}
|
||
|
||
// ---------------------------------------------------------------------------
|
||
// Sparse×Dense: Triangular MV
|
||
// ---------------------------------------------------------------------------
|
||
template <typename Scalar, int B>
|
||
static void BM_Sm_TriMV(benchmark::State& state) {
|
||
int nB, nnz;
|
||
parseArgs(state, nB, nnz);
|
||
SparseMatrix<Scalar> sm;
|
||
BlockSparseMatrix<Scalar, ColMajor, B, B> tmp;
|
||
buildUpperTriPair<Scalar, B>(nB, nnz, 1, tmp, sm);
|
||
Matrix<Scalar, Dynamic, 1> x = Matrix<Scalar, Dynamic, 1>::Random(sm.cols());
|
||
Matrix<Scalar, Dynamic, 1> y(sm.rows());
|
||
for (auto _ : state) {
|
||
y.noalias() = sm.template triangularView<Upper>() * x;
|
||
benchmark::DoNotOptimize(y.data());
|
||
}
|
||
state.counters["n"] = sm.rows();
|
||
}
|
||
|
||
// BSM triangular view, DiagIsTriangular=false: diagonal blocks treated as triangular via triangularView<>.
|
||
template <typename Scalar, int B>
|
||
static void BM_BSM_TriMV(benchmark::State& state) {
|
||
int nB, nnz;
|
||
parseArgs(state, nB, nnz);
|
||
using BSM = BlockSparseMatrix<Scalar, ColMajor, B, B>;
|
||
BSM bsm;
|
||
SparseMatrix<Scalar> smTmp;
|
||
buildUpperTriPair<Scalar, B>(nB, nnz, 1, bsm, smTmp);
|
||
Matrix<Scalar, Dynamic, 1> x = Matrix<Scalar, Dynamic, 1>::Random(bsm.cols());
|
||
Matrix<Scalar, Dynamic, 1> y(bsm.rows());
|
||
for (auto _ : state) {
|
||
y.noalias() = bsm.template triangularView<Upper, false>() * x;
|
||
benchmark::DoNotOptimize(y.data());
|
||
}
|
||
state.counters["n"] = bsm.rows();
|
||
}
|
||
|
||
// BSM triangular view, DiagIsTriangular=true: diagonal blocks are actually triangular, uses full GEMV.
|
||
template <typename Scalar, int B>
|
||
static void BM_BSM_TriMV_DiagT(benchmark::State& state) {
|
||
int nB, nnz;
|
||
parseArgs(state, nB, nnz);
|
||
using BSM = BlockSparseMatrix<Scalar, ColMajor, B, B>;
|
||
BSM bsm;
|
||
SparseMatrix<Scalar> smTmp;
|
||
buildActuallyTriPair<Scalar, B>(nB, nnz, 1, bsm, smTmp);
|
||
Matrix<Scalar, Dynamic, 1> x = Matrix<Scalar, Dynamic, 1>::Random(bsm.cols());
|
||
Matrix<Scalar, Dynamic, 1> y(bsm.rows());
|
||
for (auto _ : state) {
|
||
y.noalias() = bsm.template triangularView<Upper, true>() * x;
|
||
benchmark::DoNotOptimize(y.data());
|
||
}
|
||
state.counters["n"] = bsm.rows();
|
||
}
|
||
|
||
// Accumulating triangular MV: y += A * x.
|
||
template <typename Scalar, int B>
|
||
static void BM_Sm_TriMV_Acc(benchmark::State& state) {
|
||
int nB, nnz;
|
||
parseArgs(state, nB, nnz);
|
||
SparseMatrix<Scalar> sm;
|
||
BlockSparseMatrix<Scalar, ColMajor, B, B> tmp;
|
||
buildUpperTriPair<Scalar, B>(nB, nnz, 1, tmp, sm);
|
||
Matrix<Scalar, Dynamic, 1> x = Matrix<Scalar, Dynamic, 1>::Random(sm.cols());
|
||
Matrix<Scalar, Dynamic, 1> y = Matrix<Scalar, Dynamic, 1>::Random(sm.rows());
|
||
for (auto _ : state) {
|
||
y.noalias() += sm.template triangularView<Upper>() * x;
|
||
benchmark::DoNotOptimize(y.data());
|
||
}
|
||
state.counters["n"] = sm.rows();
|
||
}
|
||
|
||
template <typename Scalar, int B>
|
||
static void BM_BSM_TriMV_Acc(benchmark::State& state) {
|
||
int nB, nnz;
|
||
parseArgs(state, nB, nnz);
|
||
using BSM = BlockSparseMatrix<Scalar, ColMajor, B, B>;
|
||
BSM bsm;
|
||
SparseMatrix<Scalar> smTmp;
|
||
buildUpperTriPair<Scalar, B>(nB, nnz, 1, bsm, smTmp);
|
||
Matrix<Scalar, Dynamic, 1> x = Matrix<Scalar, Dynamic, 1>::Random(bsm.cols());
|
||
Matrix<Scalar, Dynamic, 1> y = Matrix<Scalar, Dynamic, 1>::Random(bsm.rows());
|
||
for (auto _ : state) {
|
||
y.noalias() += bsm.template triangularView<Upper, false>() * x;
|
||
benchmark::DoNotOptimize(y.data());
|
||
}
|
||
state.counters["n"] = bsm.rows();
|
||
}
|
||
|
||
template <typename Scalar, int B>
|
||
static void BM_BSM_TriMV_DiagT_Acc(benchmark::State& state) {
|
||
int nB, nnz;
|
||
parseArgs(state, nB, nnz);
|
||
using BSM = BlockSparseMatrix<Scalar, ColMajor, B, B>;
|
||
BSM bsm;
|
||
SparseMatrix<Scalar> smTmp;
|
||
buildActuallyTriPair<Scalar, B>(nB, nnz, 1, bsm, smTmp);
|
||
Matrix<Scalar, Dynamic, 1> x = Matrix<Scalar, Dynamic, 1>::Random(bsm.cols());
|
||
Matrix<Scalar, Dynamic, 1> y = Matrix<Scalar, Dynamic, 1>::Random(bsm.rows());
|
||
for (auto _ : state) {
|
||
y.noalias() += bsm.template triangularView<Upper, true>() * x;
|
||
benchmark::DoNotOptimize(y.data());
|
||
}
|
||
state.counters["n"] = bsm.rows();
|
||
}
|
||
|
||
// ---------------------------------------------------------------------------
|
||
// Sparse×Dense: Selfadjoint MV
|
||
// ---------------------------------------------------------------------------
|
||
template <typename Scalar, int B>
|
||
static void BM_Sm_SymmMV(benchmark::State& state) {
|
||
int nB, nnz;
|
||
parseArgs(state, nB, nnz);
|
||
SparseMatrix<Scalar> sm;
|
||
BlockSparseMatrix<Scalar, ColMajor, B, B> tmp;
|
||
buildUpperTriPair<Scalar, B>(nB, nnz, 1, tmp, sm);
|
||
Matrix<Scalar, Dynamic, 1> x = Matrix<Scalar, Dynamic, 1>::Random(sm.cols());
|
||
Matrix<Scalar, Dynamic, 1> y(sm.rows());
|
||
for (auto _ : state) {
|
||
y.noalias() = sm.template selfadjointView<Upper>() * x;
|
||
benchmark::DoNotOptimize(y.data());
|
||
}
|
||
state.counters["n"] = sm.rows();
|
||
}
|
||
|
||
// BSM selfadjoint view, general (non-Hermitian) diagonal blocks, DiagIsSA=false.
|
||
template <typename Scalar, int B>
|
||
static void BM_BSM_SymmMV(benchmark::State& state) {
|
||
int nB, nnz;
|
||
parseArgs(state, nB, nnz);
|
||
using BSM = BlockSparseMatrix<Scalar, ColMajor, B, B>;
|
||
BSM bsm;
|
||
SparseMatrix<Scalar> smTmp;
|
||
buildUpperTriPair<Scalar, B>(nB, nnz, 1, bsm, smTmp);
|
||
Matrix<Scalar, Dynamic, 1> x = Matrix<Scalar, Dynamic, 1>::Random(bsm.cols());
|
||
Matrix<Scalar, Dynamic, 1> y(bsm.rows());
|
||
for (auto _ : state) {
|
||
y.noalias() = bsm.template selfadjointView<Upper>() * x;
|
||
benchmark::DoNotOptimize(y.data());
|
||
}
|
||
state.counters["n"] = bsm.rows();
|
||
}
|
||
|
||
// BSM selfadjoint view, Hermitian diagonal blocks, DiagIsSA=false: fills diagonal via selfadjointView<>.
|
||
template <typename Scalar, int B>
|
||
static void BM_BSM_SymmMV_DiagNSA(benchmark::State& state) {
|
||
int nB, nnz;
|
||
parseArgs(state, nB, nnz);
|
||
using BSM = BlockSparseMatrix<Scalar, ColMajor, B, B>;
|
||
BSM bsm;
|
||
SparseMatrix<Scalar> smTmp;
|
||
buildHermDiagUpperTriPair<Scalar, B>(nB, nnz, 1, bsm, smTmp);
|
||
Matrix<Scalar, Dynamic, 1> x = Matrix<Scalar, Dynamic, 1>::Random(bsm.cols());
|
||
Matrix<Scalar, Dynamic, 1> y(bsm.rows());
|
||
for (auto _ : state) {
|
||
y.noalias() = bsm.template selfadjointView<Upper, false>() * x;
|
||
benchmark::DoNotOptimize(y.data());
|
||
}
|
||
state.counters["n"] = bsm.rows();
|
||
}
|
||
|
||
// BSM selfadjoint view, Hermitian diagonal blocks, DiagIsSA=true: full GEMV for diagonal blocks.
|
||
template <typename Scalar, int B>
|
||
static void BM_BSM_SymmMV_DiagSA(benchmark::State& state) {
|
||
int nB, nnz;
|
||
parseArgs(state, nB, nnz);
|
||
using BSM = BlockSparseMatrix<Scalar, ColMajor, B, B>;
|
||
BSM bsm;
|
||
SparseMatrix<Scalar> smTmp;
|
||
buildHermDiagUpperTriPair<Scalar, B>(nB, nnz, 1, bsm, smTmp);
|
||
Matrix<Scalar, Dynamic, 1> x = Matrix<Scalar, Dynamic, 1>::Random(bsm.cols());
|
||
Matrix<Scalar, Dynamic, 1> y(bsm.rows());
|
||
for (auto _ : state) {
|
||
y.noalias() = bsm.template selfadjointView<Upper, true>() * x;
|
||
benchmark::DoNotOptimize(y.data());
|
||
}
|
||
state.counters["n"] = bsm.rows();
|
||
}
|
||
|
||
// Accumulating selfadjoint MV: y += A * x.
|
||
template <typename Scalar, int B>
|
||
static void BM_Sm_SymmMV_Acc(benchmark::State& state) {
|
||
int nB, nnz;
|
||
parseArgs(state, nB, nnz);
|
||
SparseMatrix<Scalar> sm;
|
||
BlockSparseMatrix<Scalar, ColMajor, B, B> tmp;
|
||
buildUpperTriPair<Scalar, B>(nB, nnz, 1, tmp, sm);
|
||
Matrix<Scalar, Dynamic, 1> x = Matrix<Scalar, Dynamic, 1>::Random(sm.cols());
|
||
Matrix<Scalar, Dynamic, 1> y = Matrix<Scalar, Dynamic, 1>::Random(sm.rows());
|
||
for (auto _ : state) {
|
||
y.noalias() += sm.template selfadjointView<Upper>() * x;
|
||
benchmark::DoNotOptimize(y.data());
|
||
}
|
||
state.counters["n"] = sm.rows();
|
||
}
|
||
|
||
template <typename Scalar, int B>
|
||
static void BM_BSM_SymmMV_Acc(benchmark::State& state) {
|
||
int nB, nnz;
|
||
parseArgs(state, nB, nnz);
|
||
using BSM = BlockSparseMatrix<Scalar, ColMajor, B, B>;
|
||
BSM bsm;
|
||
SparseMatrix<Scalar> smTmp;
|
||
buildUpperTriPair<Scalar, B>(nB, nnz, 1, bsm, smTmp);
|
||
Matrix<Scalar, Dynamic, 1> x = Matrix<Scalar, Dynamic, 1>::Random(bsm.cols());
|
||
Matrix<Scalar, Dynamic, 1> y = Matrix<Scalar, Dynamic, 1>::Random(bsm.rows());
|
||
for (auto _ : state) {
|
||
y.noalias() += bsm.template selfadjointView<Upper>() * x;
|
||
benchmark::DoNotOptimize(y.data());
|
||
}
|
||
state.counters["n"] = bsm.rows();
|
||
}
|
||
|
||
template <typename Scalar, int B>
|
||
static void BM_BSM_SymmMV_DiagNSA_Acc(benchmark::State& state) {
|
||
int nB, nnz;
|
||
parseArgs(state, nB, nnz);
|
||
using BSM = BlockSparseMatrix<Scalar, ColMajor, B, B>;
|
||
BSM bsm;
|
||
SparseMatrix<Scalar> smTmp;
|
||
buildHermDiagUpperTriPair<Scalar, B>(nB, nnz, 1, bsm, smTmp);
|
||
Matrix<Scalar, Dynamic, 1> x = Matrix<Scalar, Dynamic, 1>::Random(bsm.cols());
|
||
Matrix<Scalar, Dynamic, 1> y = Matrix<Scalar, Dynamic, 1>::Random(bsm.rows());
|
||
for (auto _ : state) {
|
||
y.noalias() += bsm.template selfadjointView<Upper, false>() * x;
|
||
benchmark::DoNotOptimize(y.data());
|
||
}
|
||
state.counters["n"] = bsm.rows();
|
||
}
|
||
|
||
template <typename Scalar, int B>
|
||
static void BM_BSM_SymmMV_DiagSA_Acc(benchmark::State& state) {
|
||
int nB, nnz;
|
||
parseArgs(state, nB, nnz);
|
||
using BSM = BlockSparseMatrix<Scalar, ColMajor, B, B>;
|
||
BSM bsm;
|
||
SparseMatrix<Scalar> smTmp;
|
||
buildHermDiagUpperTriPair<Scalar, B>(nB, nnz, 1, bsm, smTmp);
|
||
Matrix<Scalar, Dynamic, 1> x = Matrix<Scalar, Dynamic, 1>::Random(bsm.cols());
|
||
Matrix<Scalar, Dynamic, 1> y = Matrix<Scalar, Dynamic, 1>::Random(bsm.rows());
|
||
for (auto _ : state) {
|
||
y.noalias() += bsm.template selfadjointView<Upper, true>() * x;
|
||
benchmark::DoNotOptimize(y.data());
|
||
}
|
||
state.counters["n"] = bsm.rows();
|
||
}
|
||
|
||
// ---------------------------------------------------------------------------
|
||
// Sparse×Dense: Triangular solve
|
||
// ---------------------------------------------------------------------------
|
||
template <typename Scalar, int B>
|
||
static void BM_Sm_TriSolve(benchmark::State& state) {
|
||
int nB, nnz;
|
||
parseArgs(state, nB, nnz);
|
||
SparseMatrix<Scalar> sm;
|
||
BlockSparseMatrix<Scalar, ColMajor, B, B> tmp;
|
||
buildUpperTriPair<Scalar, B>(nB, nnz, 1, tmp, sm, true);
|
||
Matrix<Scalar, Dynamic, 1> rhs = Matrix<Scalar, Dynamic, 1>::Random(sm.cols());
|
||
Matrix<Scalar, Dynamic, 1> x(sm.cols());
|
||
for (auto _ : state) {
|
||
x = rhs;
|
||
sm.template triangularView<Upper>().solveInPlace(x);
|
||
benchmark::DoNotOptimize(x.data());
|
||
}
|
||
state.counters["n"] = sm.rows();
|
||
}
|
||
|
||
template <typename Scalar, int B>
|
||
static void BM_BSM_TriSolve(benchmark::State& state) {
|
||
int nB, nnz;
|
||
parseArgs(state, nB, nnz);
|
||
using BSM = BlockSparseMatrix<Scalar, ColMajor, B, B>;
|
||
BSM bsm;
|
||
SparseMatrix<Scalar> smTmp;
|
||
buildUpperTriPair<Scalar, B>(nB, nnz, 1, bsm, smTmp, true);
|
||
Matrix<Scalar, Dynamic, 1> rhs = Matrix<Scalar, Dynamic, 1>::Random(bsm.cols());
|
||
Matrix<Scalar, Dynamic, 1> x(bsm.cols());
|
||
for (auto _ : state) {
|
||
x = rhs;
|
||
bsm.template triangularView<Upper>().solveInPlace(x);
|
||
benchmark::DoNotOptimize(x.data());
|
||
}
|
||
state.counters["n"] = bsm.rows();
|
||
}
|
||
|
||
// ---------------------------------------------------------------------------
|
||
// Sparse×Sparse: Matrix multiply
|
||
// ---------------------------------------------------------------------------
|
||
template <typename Scalar, int B>
|
||
static void BM_Sm_Sm_Mul(benchmark::State& state) {
|
||
int nB, nnz;
|
||
parseArgs(state, nB, nnz);
|
||
SparseMatrix<Scalar> smA, smB, smC;
|
||
BlockSparseMatrix<Scalar, ColMajor, B, B> tmp;
|
||
buildPair<Scalar, B>(nB, nnz, 1, tmp, smA);
|
||
buildPair<Scalar, B>(nB, nnz, 2, tmp, smB);
|
||
for (auto _ : state) {
|
||
smC = smA * smB;
|
||
benchmark::DoNotOptimize(smC.valuePtr());
|
||
}
|
||
state.counters["n"] = smA.rows();
|
||
}
|
||
|
||
template <typename Scalar, int B>
|
||
static void BM_BSM_BSM_Mul(benchmark::State& state) {
|
||
int nB, nnz;
|
||
parseArgs(state, nB, nnz);
|
||
using BSM = BlockSparseMatrix<Scalar, ColMajor, B, B>;
|
||
BSM bsmA, bsmB, bsmC;
|
||
SparseMatrix<Scalar> smTmp;
|
||
buildPair<Scalar, B>(nB, nnz, 1, bsmA, smTmp);
|
||
buildPair<Scalar, B>(nB, nnz, 2, bsmB, smTmp);
|
||
for (auto _ : state) {
|
||
bsmC = bsmA * bsmB;
|
||
benchmark::DoNotOptimize(bsmC.valuePtr());
|
||
}
|
||
state.counters["n"] = bsmA.rows();
|
||
}
|
||
|
||
// ---------------------------------------------------------------------------
|
||
// Registration
|
||
// ---------------------------------------------------------------------------
|
||
|
||
#define NS benchmark::kNanosecond
|
||
#define US benchmark::kMicrosecond
|
||
// Args: {nB, sparsity%} — nnzPerCol = max(1, pct*nB/100)
|
||
#define REG(fn, S, B) BENCHMARK(fn<S, B>)->Args({200, 1})->Args({200, 5})->Args({200, 10})
|
||
|
||
#define BENCH_TYPE(S, B) \
|
||
REG(BM_Sm_Sm_Add, S, B)->Unit(US); \
|
||
REG(BM_BSM_BSM_Add, S, B)->Unit(US); \
|
||
REG(BM_Sm_SpMV, S, B)->Unit(NS); \
|
||
REG(BM_BSM_SpMV, S, B)->Unit(NS); \
|
||
REG(BM_Sm_SpMV_Acc, S, B)->Unit(NS); \
|
||
REG(BM_BSM_SpMV_Acc, S, B)->Unit(NS); \
|
||
REG(BM_Sm_TriMV, S, B)->Unit(NS); \
|
||
REG(BM_BSM_TriMV, S, B)->Unit(NS); \
|
||
REG(BM_BSM_TriMV_DiagT, S, B)->Unit(NS); \
|
||
REG(BM_Sm_TriMV_Acc, S, B)->Unit(NS); \
|
||
REG(BM_BSM_TriMV_Acc, S, B)->Unit(NS); \
|
||
REG(BM_BSM_TriMV_DiagT_Acc, S, B)->Unit(NS); \
|
||
REG(BM_Sm_SymmMV, S, B)->Unit(NS); \
|
||
REG(BM_BSM_SymmMV, S, B)->Unit(NS); \
|
||
REG(BM_BSM_SymmMV_DiagNSA, S, B)->Unit(NS); \
|
||
REG(BM_BSM_SymmMV_DiagSA, S, B)->Unit(NS); \
|
||
REG(BM_Sm_SymmMV_Acc, S, B)->Unit(NS); \
|
||
REG(BM_BSM_SymmMV_Acc, S, B)->Unit(NS); \
|
||
REG(BM_BSM_SymmMV_DiagNSA_Acc, S, B)->Unit(NS); \
|
||
REG(BM_BSM_SymmMV_DiagSA_Acc, S, B)->Unit(NS); \
|
||
REG(BM_Sm_TriSolve, S, B)->Unit(NS); \
|
||
REG(BM_BSM_TriSolve, S, B)->Unit(NS); \
|
||
REG(BM_Sm_Sm_Mul, S, B)->Unit(US); \
|
||
REG(BM_BSM_BSM_Mul, S, B)->Unit(US);
|
||
|
||
BENCH_TYPE(float, 2)
|
||
BENCH_TYPE(cf, 2)
|
||
BENCH_TYPE(double, 2)
|
||
BENCH_TYPE(cd, 2)
|
||
BENCH_TYPE(float, 3)
|
||
BENCH_TYPE(cf, 3)
|
||
BENCH_TYPE(double, 3)
|
||
BENCH_TYPE(cd, 3)
|
||
BENCH_TYPE(float, 4)
|
||
BENCH_TYPE(cf, 4)
|
||
BENCH_TYPE(double, 4)
|
||
BENCH_TYPE(cd, 4)
|