From cd70e0c4a9efdf9b2f175a3fa5f43319db426ca0 Mon Sep 17 00:00:00 2001 From: Charles Schlosser Date: Thu, 2 Jul 2026 22:20:12 +0000 Subject: [PATCH] BlockSparseMatrix: add block-sparse data structure and operations libeigen/eigen!2644 --- Eigen/SparseCore | 1 + Eigen/src/SparseCore/BlockSparseMatrix.h | 1846 +++++++++++++++++ benchmarks/Sparse/CMakeLists.txt | 1 + benchmarks/Sparse/bench_block_sparse.cpp | 672 ++++++ test/CMakeLists.txt | 1 + test/block_sparse_matrix.cpp | 730 +++++++ .../Eigen/src/SparseExtra/BlockSparseMatrix.h | 932 --------- 7 files changed, 3251 insertions(+), 932 deletions(-) create mode 100644 Eigen/src/SparseCore/BlockSparseMatrix.h create mode 100644 benchmarks/Sparse/bench_block_sparse.cpp create mode 100644 test/block_sparse_matrix.cpp delete mode 100644 unsupported/Eigen/src/SparseExtra/BlockSparseMatrix.h diff --git a/Eigen/SparseCore b/Eigen/SparseCore index 8891a180e..641101348 100644 --- a/Eigen/SparseCore +++ b/Eigen/SparseCore @@ -40,6 +40,7 @@ #include "src/SparseCore/AmbiVector.h" #include "src/SparseCore/SparseCompressedBase.h" #include "src/SparseCore/SparseMatrix.h" +#include "src/SparseCore/BlockSparseMatrix.h" #include "src/SparseCore/SparseMap.h" #include "src/SparseCore/SparseVector.h" #include "src/SparseCore/SparseRef.h" diff --git a/Eigen/src/SparseCore/BlockSparseMatrix.h b/Eigen/src/SparseCore/BlockSparseMatrix.h new file mode 100644 index 000000000..ee5ad8733 --- /dev/null +++ b/Eigen/src/SparseCore/BlockSparseMatrix.h @@ -0,0 +1,1846 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. +// SPDX-FileCopyrightText: The Eigen Authors +// SPDX-License-Identifier: MPL-2.0 + +#ifndef EIGEN_BLOCKSPARSEMATRIX_H +#define EIGEN_BLOCKSPARSEMATRIX_H + +// IWYU pragma: private +#include "./InternalHeaderCheck.h" + +#include +#include +#include + +namespace Eigen { + +// Forward declarations +template +class BlockSparseTriangularView; +template +class BlockSparseSelfAdjointView; +template +class BlockSparseMatrix; + +/** Storage-kind tag for BlockSparseMatrix. */ +struct BlockSparse {}; + +/** Evaluator shape tag for BlockSparseMatrix product dispatch. */ +struct BlockSparseShape { + static std::string debugName() { return "BlockSparseShape"; } +}; + +namespace internal { +// Returns m.adjoint() when Conj==true, m.transpose() otherwise. +// SFINAE overloads keep the return type concrete under C++14 (no if constexpr). +template +std::enable_if_t().adjoint())> adjoint_if(const T& m) { + return m.adjoint(); +} + +template +std::enable_if_t().transpose())> adjoint_if(const T& m) { + return m.transpose(); +} +template <> +struct storage_kind_to_evaluator_kind { + using Kind = IndexBased; +}; + +template <> +struct storage_kind_to_shape { + using Shape = BlockSparseShape; +}; + +template +struct traits> { + using Scalar = Scalar_; + using StorageIndex = StorageIndex_; + using StorageKind = BlockSparse; + using XprKind = MatrixXpr; + + static constexpr Index RowsAtCompileTime = Dynamic; + static constexpr Index ColsAtCompileTime = Dynamic; + static constexpr Index MaxRowsAtCompileTime = Dynamic; + static constexpr Index MaxColsAtCompileTime = Dynamic; + static constexpr int Options = Options_; + static constexpr unsigned int Flags = Options_ | NestByRefBit | LvalueBit; +}; + +} // namespace internal + +/** \class BlockTriplet + * \ingroup SparseCore_Module + * \brief A (blockRow, blockCol, blockValue) triplet for assembling a BlockSparseMatrix. + * + * Coordinates are in \em block space, not element space. + * + * \tparam Scalar_ Numeric scalar type. + * \tparam BlockRows_ Number of rows in each block. + * \tparam BlockCols_ Number of columns in each block. + * \tparam StorageIndex_ Signed integer index type (default: int). + */ +template +class BlockTriplet { + public: + using Scalar = Scalar_; + using StorageIndex = StorageIndex_; + using BlockType = Matrix; + using BlockMapType = Map; + using ConstBlockMapType = Map; + + static constexpr int BlockSize = BlockRows_ * BlockCols_; + + BlockTriplet() = default; + + BlockTriplet(StorageIndex blockRow, StorageIndex blockCol, const BlockType& block) + : m_row(blockRow), m_col(blockCol) { + BlockMapType{m_value} = block; + } + + StorageIndex row() const { return m_row; } + StorageIndex col() const { return m_col; } + // Implicitly usable wherever a MatrixBase expression is expected. + ConstBlockMapType value() const { return ConstBlockMapType(m_value); } + + private: + StorageIndex m_row = 0; + StorageIndex m_col = 0; + // Flat array avoids the alignment padding that a Matrix<> member would incur. + Scalar m_value[BlockSize]; +}; + +/** \class BlockSparseMatrix + * \ingroup SparseCore_Module + * \brief A sparse matrix whose stored nonzeros are fixed-size dense blocks. + * + * Each nonzero entry is a \c BlockRows x \c BlockCols dense matrix. The + * block sparsity pattern is stored in block-level compressed-column + * (ColMajor) or compressed-row (RowMajor) format. + * + * \tparam Scalar_ Numeric scalar type. + * \tparam Options_ ColMajor (0) or RowMajor. Controls both the outer + * iteration direction over blocks and the storage + * layout within each block. Note: vector-shaped blocks + * must use a compatible order — \c BlockCols_=1 requires + * ColMajor, \c BlockRows_=1 requires RowMajor. + * \tparam BlockRows_ Rows per block; must be a fixed positive integer. + * \tparam BlockCols_ Columns per block; must be a fixed positive integer. + * \tparam StorageIndex_ Signed integer type for internal index arrays + * (default: int). + * + * ### Assembly + * Populate the matrix via setFromTriplets(), passing an iterator range of + * BlockTriplet objects in block coordinates. Triplets with the same + * (blockRow, blockCol) are summed. + * + * ### Arithmetic + * Addition, subtraction, and matrix product are supported between compatible + * BlockSparseMatrix instances. Scalar multiplication is also available. + * + * The matrix product \c C = A * B requires + * \c A.BlockCols == B.BlockRows (enforced at compile time by the template + * constraint) and \c A.blockCols() == B.blockRows() (checked at runtime). + * The result has block type \c Matrix. + * + * ### Conversion + * toSparse() converts to a standard SparseMatrix with element-level + * sparsity. fromSparse() reconstructs a BlockSparseMatrix from an + * element-level SparseMatrix whose dimensions are divisible by BlockRows + * and BlockCols. An implicit conversion operator to SparseMatrix is + * also provided. + */ +template +class BlockSparseMatrix + : public EigenBase> { + EIGEN_STATIC_ASSERT(BlockRows_ >= 1, BLOCKROWS_MUST_BE_A_POSITIVE_COMPILE_TIME_SIZE) + EIGEN_STATIC_ASSERT(BlockCols_ >= 1, BLOCKCOLS_MUST_BE_A_POSITIVE_COMPILE_TIME_SIZE) + EIGEN_STATIC_ASSERT(std::is_integral::value&& std::is_signed::value, + STORAGEINDEX_MUST_BE_A_SIGNED_INTEGRAL_TYPE) + // Eigen's Matrix<> requires: a column vector (Cols==1, Rows>1) must be ColMajor; + // a row vector (Rows==1, Cols>1) must be RowMajor. Guard those cases here so + // the error fires at BlockSparseMatrix instantiation rather than inside BlockType. + EIGEN_STATIC_ASSERT(!(BlockCols_ == 1 && BlockRows_ != 1 && bool(Options_ & RowMajorBit)), + INVALID_MATRIX_TEMPLATE_PARAMETERS) + EIGEN_STATIC_ASSERT(!(BlockRows_ == 1 && BlockCols_ != 1 && !bool(Options_ & RowMajorBit)), + INVALID_MATRIX_TEMPLATE_PARAMETERS) + + public: + // ------------------------------------------------------------------------- + // Type aliases & compile-time constants + // ------------------------------------------------------------------------- + using Scalar = Scalar_; + using StorageIndex = StorageIndex_; + using BlockType = Matrix; + using TripletType = BlockTriplet; + + static constexpr int Options = Options_; + static constexpr Index BlockRows = BlockRows_; + static constexpr Index BlockCols = BlockCols_; + static constexpr bool IsRowMajor = Options_ & RowMajorBit; + static constexpr Index BlockSize = BlockRows_ * BlockCols_; + + // If one block occupies a power-of-two number of bytes, and the values array + // is Eigen-allocated (guaranteed aligned to EIGEN_MAX_ALIGN_BYTES), then every + // block pointer is aligned to min(BlockBytes, EIGEN_MAX_ALIGN_BYTES). + static constexpr std::size_t BlockBytes = std::size_t(BlockSize) * sizeof(Scalar); + static constexpr int BlockMapAlignment = ((BlockBytes & (BlockBytes - 1)) == 0 && BlockBytes >= 8) + ? int(numext::mini(BlockBytes, std::size_t(EIGEN_MAX_ALIGN_BYTES))) + : 0; + + using BlockMap = Map; + using ConstBlockMap = Map; + + // ------------------------------------------------------------------------- + // Constructors / copy / move + // ------------------------------------------------------------------------- + + /** Default constructor; creates a 0×0 matrix. */ + BlockSparseMatrix() = default; + + /** Construct a zero matrix with the given number of block-rows and block-columns. */ + BlockSparseMatrix(Index blockRows, Index blockCols) + : m_blockOuterSize(IsRowMajor ? blockRows : blockCols), m_blockInnerSize(IsRowMajor ? blockCols : blockRows) {} + + // ------------------------------------------------------------------------- + // Dimensions + // ------------------------------------------------------------------------- + + /** Total number of element rows. */ + Index rows() const noexcept { return (IsRowMajor ? m_blockOuterSize : m_blockInnerSize) * BlockRows_; } + /** Total number of element columns. */ + Index cols() const noexcept { return (IsRowMajor ? m_blockInnerSize : m_blockOuterSize) * BlockCols_; } + + /** Number of block-rows. */ + Index blockRows() const { return IsRowMajor ? m_blockOuterSize : m_blockInnerSize; } + /** Number of block-columns. */ + Index blockCols() const { return IsRowMajor ? m_blockInnerSize : m_blockOuterSize; } + + /** Outer block dimension (block-cols for ColMajor, block-rows for RowMajor). */ + Index blockOuterSize() const { return m_blockOuterSize; } + /** Inner block dimension (block-rows for ColMajor, block-cols for RowMajor). */ + Index blockInnerSize() const { return m_blockInnerSize; } + + /** Element-level outer size: cols() for ColMajor, rows() for RowMajor. */ + Index outerSize() const { return IsRowMajor ? rows() : cols(); } + /** Element-level inner size: rows() for ColMajor, cols() for RowMajor. */ + Index innerSize() const { return IsRowMajor ? cols() : rows(); } + + /** Number of stored (structurally non-zero) blocks. */ + Index nonZeroBlocks() const { return m_outerIndex(m_blockOuterSize); } + /** Total number of stored scalar coefficients (= nonZeroBlocks() * BlockRows * BlockCols). */ + Index nonZeros() const { return nonZeroBlocks() * BlockSize; } + /** Number of blocks for which storage is currently allocated (capacity). */ + Index allocatedBlocks() const { return Index(m_innerIndex.size()); } + + // ------------------------------------------------------------------------- + // Raw pointer access (for interoperability) + // ------------------------------------------------------------------------- + const StorageIndex* outerIndexPtr() const { return m_outerIndex.data(); } + StorageIndex* outerIndexPtr() { return m_outerIndex.data(); } + const StorageIndex* innerIndexPtr() const { return m_innerIndex.data(); } + StorageIndex* innerIndexPtr() { return m_innerIndex.data(); } + const Scalar* valuePtr() const { return m_values.data(); } + Scalar* valuePtr() { return m_values.data(); } + + // ------------------------------------------------------------------------- + // Block access by sequential nonzero index + // ------------------------------------------------------------------------- + + /** Read-only Map to the \a k-th stored block (block storage follows \c Options_). */ + ConstBlockMap blockRef(Index k) const { return ConstBlockMap(m_values.data() + k * BlockSize); } + /** Mutable Map to the \a k-th stored block. */ + BlockMap blockRef(Index k) { return BlockMap(m_values.data() + k * BlockSize); } + + // ------------------------------------------------------------------------- + // Inner iterator over blocks within one outer vector + // ------------------------------------------------------------------------- + + /** \brief Iterates over stored blocks in outer vector \a outer. + * + * Usage mirrors SparseMatrix::InnerIterator but value() returns a + * Map to a BlockRows×BlockCols matrix, not a scalar. + */ + class InnerIterator { + public: + EIGEN_STRONG_INLINE InnerIterator(const BlockSparseMatrix& mat, Index outer) + : m_mat(mat), m_id(mat.m_outerIndex(outer)), m_end(mat.m_outerIndex(outer + 1)), m_outer(outer) {} + + EIGEN_STRONG_INLINE operator bool() const { return m_id < m_end; } + EIGEN_STRONG_INLINE InnerIterator& operator++() { + ++m_id; + return *this; + } + + /** Current block outer index (block-col for ColMajor, block-row for RowMajor). */ + EIGEN_STRONG_INLINE Index outer() const { return m_outer; } + /** Current block inner index (block-row for ColMajor, block-col for RowMajor). */ + EIGEN_STRONG_INLINE Index index() const { return m_mat.m_innerIndex(m_id); } + /** Block-row of the current block. */ + EIGEN_STRONG_INLINE Index blockRow() const { return IsRowMajor ? m_outer : index(); } + /** Block-column of the current block. */ + EIGEN_STRONG_INLINE Index blockCol() const { return IsRowMajor ? index() : m_outer; } + + /** Read-only Map to the current block value. */ + EIGEN_STRONG_INLINE ConstBlockMap value() const { return m_mat.blockRef(m_id); } + /** Mutable Map to the current block value. */ + EIGEN_STRONG_INLINE BlockMap valueRef() { + return BlockMap(const_cast(m_mat.m_values.data()) + m_id * BlockSize); + } + + private: + const BlockSparseMatrix& m_mat; + Index m_id; + Index m_end; + Index m_outer; + }; + + // ------------------------------------------------------------------------- + // Resize / clear + // ------------------------------------------------------------------------- + + /** Resize to \a blockRows × \a blockCols blocks and set the logical nnz to zero. + * Allocated block storage is retained; call squeeze() to release it. */ + void resize(Index blockRows, Index blockCols) { + m_blockOuterSize = IsRowMajor ? blockRows : blockCols; + m_blockInnerSize = IsRowMajor ? blockCols : blockRows; + m_outerIndex.resize(m_blockOuterSize + 1); + m_outerIndex.setZero(); + } + + /** Clear all stored blocks (logical nnz → 0) while keeping dimensions and allocated storage. */ + void setZero() { + m_outerIndex.resize(m_blockOuterSize + 1); + m_outerIndex.setZero(); + } + + /** Pre-allocate storage for at least \a n blocks without changing the logical sparsity pattern. + * Existing block data is preserved up to min(n, nonZeroBlocks()). */ + void reserve(Index n) { + if (n > Index(m_innerIndex.size())) conservativeResizeBlockStorage_(n); + } + + /** Release any excess allocated block storage so that allocatedBlocks() == nonZeroBlocks(). */ + void squeeze() { + Index nnz = nonZeroBlocks(); + if (nnz < Index(m_innerIndex.size())) conservativeResizeBlockStorage_(nnz); + } + + /** Fill the matrix with the block identity: the min(blockRows,blockCols) diagonal blocks + * are set to the B×B identity; all other blocks are absent. + * + * \pre BlockRows == BlockCols (square blocks). + */ + void setIdentity() { + EIGEN_STATIC_ASSERT(BlockRows_ == BlockCols_, THIS_METHOD_IS_ONLY_FOR_SQUARE_BLOCK_MATRICES) + Index n = (std::min)(m_blockOuterSize, m_blockInnerSize); + m_outerIndex.resize(m_blockOuterSize + 1); + resizeBlockStorage_(n); + for (Index i = 0; i <= m_blockOuterSize; ++i) m_outerIndex(i) = StorageIndex((std::min)(i, n)); + for (StorageIndex i = 0; i < n; ++i) { + m_innerIndex(i) = i; + blockRef(i).setIdentity(); + } + } + + /** Initialize the block structure directly from compressed outer/inner index arrays, + * zero-initializing all block values. + * + * \p outerPtr has size blockCols+1 (ColMajor) or blockRows+1 (RowMajor). + * \p innerPtr has size nnzBlocks. + */ + void setFromOuterInner(Index blockRows, Index blockCols, Index nnzBlocks, const StorageIndex_* outerPtr, + const StorageIndex_* innerPtr) { + m_blockOuterSize = IsRowMajor ? blockRows : blockCols; + m_blockInnerSize = IsRowMajor ? blockCols : blockRows; + m_outerIndex = Map(outerPtr, m_blockOuterSize + 1); + m_innerIndex = Map(innerPtr, nnzBlocks); + resizeBlockStorage_(nnzBlocks); + m_values.setZero(); + } + + // ------------------------------------------------------------------------- + // Assembly + // ------------------------------------------------------------------------- + + /** Fill the matrix from an iterator range of BlockTriplet objects. + * + * Triplet coordinates are in block space. Triplets with the same + * (blockRow, blockCol) pair are summed (their block values are added). + * The input range may be in any order. + * + * \tparam InputIterator Must dereference to a type with \c row(), + * \c col(), and \c value() members, matching + * BlockTriplet's interface. + */ + template + void setFromTriplets(InputIterator begin, InputIterator end); + + // ------------------------------------------------------------------------- + // Conversion to / from SparseMatrix + // ------------------------------------------------------------------------- + + /** Convert to a SparseMatrix with scalar nonzeros. + * + * Each stored block of size BlockRows×BlockCols expands into up to + * BlockRows*BlockCols scalar nonzeros. The resulting SparseMatrix has + * the same storage order as \c *this. + */ + SparseMatrix toSparse() const; + + /** Construct a BlockSparseMatrix from an element-level SparseMatrix. + * + * \pre \c sp.rows() % BlockRows == 0 and \c sp.cols() % BlockCols == 0. + * + * Each scalar entry \c sp(i,j) is placed into position + * \c (i%BlockRows, j%BlockCols) of block \c (i/BlockRows, j/BlockCols). + * Multiple entries mapping to the same element of the same block are + * accumulated with \c +=. + */ + static BlockSparseMatrix fromSparse(const SparseMatrix& sp); + + /** Implicit conversion to SparseMatrix. */ + operator SparseMatrix() const { return toSparse(); } + + // ------------------------------------------------------------------------- + // Element access + // ------------------------------------------------------------------------- + + /** Read element \c (row, col); returns 0 if no block covers that position. */ + Scalar coeff(Index row, Index col) const { + eigen_assert(row >= 0 && row < rows() && col >= 0 && col < cols()); + Index bOuter = IsRowMajor ? (row / BlockRows_) : (col / BlockCols_); + Index bInner = IsRowMajor ? (col / BlockCols_) : (row / BlockRows_); + Index localRow = row % BlockRows_; + Index localCol = col % BlockCols_; + const StorageIndex* beg = m_innerIndex.data() + m_outerIndex(bOuter); + const StorageIndex* fin = m_innerIndex.data() + m_outerIndex(bOuter + 1); + const StorageIndex* it = std::lower_bound(beg, fin, StorageIndex(bInner)); + if (it == fin || *it != bInner) return Scalar(0); + return blockRef(static_cast(it - m_innerIndex.data()))(localRow, localCol); + } + + /** Extract the main scalar diagonal as a dense vector. + * + * Iterates outer slices once and binary-searches for the diagonal block in each + * slice, then copies the relevant entries from that block — one search per outer + * slice (square blocks) or per unique inner-block boundary (non-square blocks), + * vs. one search per scalar element for the coeff-by-coeff approach. + */ + Matrix diagonal() const { + constexpr Index OuterB = IsRowMajor ? BlockRows_ : BlockCols_; + constexpr Index InnerB = IsRowMajor ? BlockCols_ : BlockRows_; + const Index diagSize = numext::mini(rows(), cols()); + Matrix diag = Matrix::Zero(diagSize); + + for (Index out = 0; out < m_blockOuterSize; ++out) { + const Index scalarOuterBegin = out * OuterB; + if (scalarOuterBegin >= diagSize) break; + const Index scalarOuterEnd = numext::mini(scalarOuterBegin + OuterB, diagSize); + + // Group consecutive scalar positions that share the same inner block, then + // binary-search once per group rather than once per scalar element. + // For square blocks this loop runs exactly once per outer slice. + Index i = scalarOuterBegin; + while (i < scalarOuterEnd) { + const Index bInner = i / InnerB; + const Index groupEnd = numext::mini((bInner + 1) * InnerB, scalarOuterEnd); + + const StorageIndex* beg = m_innerIndex.data() + m_outerIndex(out); + const StorageIndex* fin = m_innerIndex.data() + m_outerIndex(out + 1); + const StorageIndex* it = std::lower_bound(beg, fin, StorageIndex(bInner)); + if (it != fin && *it == StorageIndex(bInner)) { + const ConstBlockMap blk = blockRef(static_cast(it - m_innerIndex.data())); + for (Index j = i; j < groupEnd; ++j) { + const Index localRow = IsRowMajor ? (j % OuterB) : (j % InnerB); + const Index localCol = IsRowMajor ? (j % InnerB) : (j % OuterB); + diag(j) = blk(localRow, localCol); + } + } + i = groupEnd; + } + } + return diag; + } + + // ------------------------------------------------------------------------- + // Arithmetic + // ------------------------------------------------------------------------- + + /** Element-wise addition. Both matrices must have the same block dimensions. */ + BlockSparseMatrix operator+(const BlockSparseMatrix& other) const { return disjunctionWith_(other, AddOp_{}); } + + /** Element-wise subtraction. */ + BlockSparseMatrix operator-(const BlockSparseMatrix& other) const { return disjunctionWith_(other, SubOp_{}); } + + /** Element-wise product. Only blocks present in \em both operands contribute to the result. */ + BlockSparseMatrix cwiseProduct(const BlockSparseMatrix& other) const { + return conjunctionWith_(other, CwiseMulOp_{}); + } + + /** Applies a scalar unary functor to every stored nonzero, preserving the sparsity pattern. */ + template + BlockSparseMatrix unaryExpr(ScalarFunc func) const { + return withValues_([&func](const auto& v) { return v.unaryExpr(func); }); + } + + /** Applies a scalar binary functor with union sparsity. + * + * \p func must provide three members: + * \code + * Scalar func(Scalar a, Scalar b) // both present + * Scalar func.lhs(Scalar a) // only lhs present; rhs is implicitly zero + * Scalar func.rhs(Scalar b) // only rhs present; lhs is implicitly zero + * \endcode + */ + template + BlockSparseMatrix disjunctionExpr(const BlockSparseMatrix& other, ScalarFunc func) const { + return disjunctionWith_(other, DisjExprAdapter_{func}); + } + + /** Applies a scalar binary functor with intersection sparsity: only block positions present + * in \em both matrices contribute; \p func is called as \c func(Scalar a, Scalar b). */ + template + BlockSparseMatrix conjunctionExpr(const BlockSparseMatrix& other, ScalarFunc func) const { + return conjunctionWith_(other, [&func](const auto& a, const auto& b) { return a.binaryExpr(b, func); }); + } + + /** Unary negation. */ + BlockSparseMatrix operator-() const { + return withValues_([](const auto& v) { return -v; }); + } + + BlockSparseMatrix& operator+=(const BlockSparseMatrix& other) { return *this = *this + other; } + BlockSparseMatrix& operator-=(const BlockSparseMatrix& other) { return *this = *this - other; } + + /** Scalar multiplication (returns a new matrix). */ + BlockSparseMatrix operator*(const Scalar& s) const { + return withValues_([&s](const auto& v) { return v * s; }); + } + BlockSparseMatrix& operator*=(const Scalar& s) { + m_values.head(nonZeros()) *= s; + return *this; + } + BlockSparseMatrix operator/(const Scalar& s) const { + return withValues_([&s](const auto& v) { return v / s; }); + } + BlockSparseMatrix& operator/=(const Scalar& s) { return *this *= (Scalar(1) / s); } + + /** Scalar-on-left multiplication. */ + friend BlockSparseMatrix operator*(const Scalar& s, const BlockSparseMatrix& m) { return m * s; } + + /** Block-sparse times dense matrix (or vector) product. + * + * Returns a lazy \c Product<> expression evaluated via \c generic_product_impl. + * This enables fused accumulation: + * \code + * b.noalias() = A * x; // no temporary + * b.noalias() += A * x; // fused add + * b.noalias() += alpha * (A * x); // scale then add via evaluator + * \endcode + * + * \pre \c this->cols() == rhs.rows(). + * \pre Scalar types must match. + * \warning The result uses \c AliasFreeProduct, so assignment goes directly + * through \c generic_product_impl::evalTo with no aliasing temporary. + * \c x = A * x silently corrupts; use an explicit temporary if needed. + */ + template + Product operator*(const MatrixBase& rhs) const { + EIGEN_STATIC_ASSERT( + (std::is_same::value), + YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) + return Product(*this, rhs.derived()); + } + + /** Dense matrix (or vector) times block-sparse product (hidden friend). + * + * Returns a lazy \c Product<> expression; evaluated via \c generic_product_impl. + * + * \pre \c lhs.cols() == bsm.rows(). + * \pre Scalar types must match. + * \warning The result uses \c AliasFreeProduct, so assignment goes directly + * through \c generic_product_impl::evalTo with no aliasing temporary. + * \c x = x * A silently corrupts; use an explicit temporary if needed. + */ + template + friend Product operator*(const MatrixBase& lhs, + const BlockSparseMatrix& bsm) { + EIGEN_STATIC_ASSERT( + (std::is_same::value), + YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) + return Product(lhs.derived(), bsm); + } + + /** Block-sparse matrix product. + * + * \tparam RhsBlockCols Block-column count of the right-hand side. + * + * The left-hand side has block type \c Matrix + * and the right-hand side must have block type + * \c Matrix (enforced by the template + * parameter). The result has block type + * \c Matrix. + * + * \pre \c this->blockCols() == rhs.blockRows() (checked at runtime). + */ + template + BlockSparseMatrix operator*( + const BlockSparseMatrix& rhs) const; + + // ------------------------------------------------------------------------- + // Approximate equality (useful for testing) + // ------------------------------------------------------------------------- + bool isApprox(const BlockSparseMatrix& other, + const typename NumTraits::Real& prec = NumTraits::dummy_precision()) const { + using RealScalar = typename NumTraits::Real; + // Frobenius-norm comparison, matching SparseMatrixBase::isApprox semantics + // but computed directly from the block values — no scalar-level SparseMatrix + // materialization of either operand. Explicit zero blocks contribute 0 to + // every norm, so the result is independent of structural differences. + RealScalar n2a = m_values.head(nonZeros()).matrix().squaredNorm(); + RealScalar n2b = other.m_values.head(other.nonZeros()).matrix().squaredNorm(); + BlockSparseMatrix diff = *this - other; + RealScalar d2 = diff.m_values.head(diff.nonZeros()).matrix().squaredNorm(); + return d2 <= prec * prec * numext::mini(n2a, n2b); + } + + // ------------------------------------------------------------------------- + // Transpose / adjoint + // ------------------------------------------------------------------------- + + /** Returns a new BSM whose block dimensions are swapped (BlockRows ↔ BlockCols) + * and each stored block is transposed. Storage order is preserved. */ + BlockSparseMatrix transpose() const; + + /** Returns a new BSM whose block dimensions are swapped and each stored block + * is conjugate-transposed (adjoint). Identical to transpose() for real scalars. */ + BlockSparseMatrix adjoint() const; + + // ------------------------------------------------------------------------- + // View factories + // ------------------------------------------------------------------------- + + /** Returns a block-level triangular view. + * + * \tparam Mode \c Eigen::Upper or \c Eigen::Lower. + * \tparam DiagIsTriangular When \c false (default), diagonal blocks are + * treated as triangular regardless of what is stored in the unused + * triangle — products use \c triangularView on the block and \c eval() + * explicitly zeros the unused triangle. When \c true, the caller + * guarantees that the unused triangle of every diagonal block is already + * zero; products then use a full vectorised GEMV (faster, no zeroing). + */ + template + BlockSparseTriangularView triangularView() const { + return BlockSparseTriangularView(*this); + } + + /** Returns a block-level self-adjoint view. + * + * Only the triangle selected by \p UpLo is read; the opposite triangle is + * reconstructed on-the-fly as the adjoint of each stored off-diagonal block. + * + * \tparam UpLo \c Eigen::Upper or \c Eigen::Lower. + * \tparam DiagIsSelfAdjoint Set to \c true when every diagonal block is + * itself Hermitian \em and both triangles are explicitly stored. + * The dense sub-product then uses a plain product rather than + * \c selfadjointView, which is more efficient for the small + * fixed-size blocks typical here. When \c false (the default), + * only the \p UpLo triangle of each diagonal block is assumed + * valid; \c selfadjointView() reconstructs the full block. + * + * \pre BlockRows == BlockCols (diagonal blocks must be square). + */ + template + BlockSparseSelfAdjointView selfadjointView() const { + EIGEN_STATIC_ASSERT(BlockRows_ == BlockCols_, THIS_METHOD_IS_ONLY_FOR_SQUARE_BLOCK_MATRICES) + return BlockSparseSelfAdjointView(*this); + } + + private: + struct AddOp_ { + template + BlockType operator()(const A& a, const B& b) const { + return a + b; + } + template + BlockType lhs(const A& a) const { + return a; + } + template + BlockType rhs(const B& b) const { + return b; + } + }; + + struct SubOp_ { + template + BlockType operator()(const A& a, const B& b) const { + return a - b; + } + template + BlockType lhs(const A& a) const { + return a; + } + template + BlockType rhs(const B& b) const { + return -b; + } + }; + + // Scalar-to-block adapter for disjunctionExpr: translates scalar functor with lhs/rhs methods + // to the block level. + template + struct DisjExprAdapter_ { + ScalarFunc func_; + template + BlockType operator()(const A& a, const B& b) const { + return a.binaryExpr(b, func_); + } + template + BlockType lhs(const A& a) const { + return a.unaryExpr([this](const Scalar& x) { return func_.lhs(x); }); + } + template + BlockType rhs(const B& b) const { + return b.unaryExpr([this](const Scalar& x) { return func_.rhs(x); }); + } + }; + + struct CwiseMulOp_ { + template + BlockType operator()(const A& a, const B& b) const { + return a.cwiseProduct(b); + } + }; + + // Returns a copy with the same sparsity structure but m_values replaced by f(m_values). + // f receives the flat Eigen Array of the logical (non-zero) coefficients and returns any + // compatible expression. Only the structure is copied — the source values are never + // duplicated, and unused tail capacity is neither copied nor evaluated. + template + BlockSparseMatrix withValues_(F f) const { + Index nnz = nonZeroBlocks(); + BlockSparseMatrix result(blockRows(), blockCols()); + result.m_outerIndex = m_outerIndex; + result.m_innerIndex = m_innerIndex.head(nnz); + result.m_values = f(m_values.head(nnz * BlockSize)); + return result; + } + + // Disjunction (union-pattern): result has a block wherever *this OR other has one. + // lhs-only: block copied from *this unchanged + // rhs-only: op(b) — unary overload of op + // both: op(a, b) — binary overload of op + template + BlockSparseMatrix disjunctionWith_(const BlockSparseMatrix& other, Op op) const { + eigen_assert(blockRows() == other.blockRows() && blockCols() == other.blockCols() && + "BlockSparseMatrix size mismatch"); + BlockSparseMatrix result(blockRows(), blockCols()); + result.resizeBlockStorage_(nonZeroBlocks() + other.nonZeroBlocks()); + Index nnz = 0; + for (Index j = 0; j < m_blockOuterSize; ++j) { + result.m_outerIndex(j) = StorageIndex_(nnz); + Index aId = m_outerIndex(j); + Index aEnd = m_outerIndex(j + 1); + Index bId = other.m_outerIndex(j); + Index bEnd = other.m_outerIndex(j + 1); + while (aId < aEnd || bId < bEnd) { + bool hasA = aId < aEnd; + bool hasB = bId < bEnd; + Index aInner = hasA ? Index(m_innerIndex(aId)) : -1; + Index bInner = hasB ? Index(other.m_innerIndex(bId)) : -1; + // Write the result block in place; the op overloads return a BlockType, + // which is assigned straight into the result's storage Map. + if (hasA && (!hasB || aInner < bInner)) { + result.m_innerIndex(nnz) = StorageIndex_(aInner); + result.blockRef(nnz) = op.lhs(blockRef(aId++)); + } else if (hasB && (!hasA || bInner < aInner)) { + result.m_innerIndex(nnz) = StorageIndex_(bInner); + result.blockRef(nnz) = op.rhs(other.blockRef(bId++)); + } else { + result.m_innerIndex(nnz) = StorageIndex_(aInner); + result.blockRef(nnz) = op(blockRef(aId++), other.blockRef(bId++)); + } + ++nnz; + } + } + result.m_outerIndex(m_blockOuterSize) = StorageIndex_(nnz); + result.conservativeResizeBlockStorage_(nnz); + return result; + } + + // Conjunction (intersection-pattern): result has a block only where *this AND other both have one. + template + BlockSparseMatrix conjunctionWith_(const BlockSparseMatrix& other, BinaryOp func) const { + eigen_assert(blockRows() == other.blockRows() && blockCols() == other.blockCols() && + "BlockSparseMatrix size mismatch"); + BlockSparseMatrix result(blockRows(), blockCols()); + result.resizeBlockStorage_((std::min)(nonZeroBlocks(), other.nonZeroBlocks())); + Index nnz = 0; + for (Index j = 0; j < m_blockOuterSize; ++j) { + result.m_outerIndex(j) = StorageIndex_(nnz); + Index aId = m_outerIndex(j); + Index aEnd = m_outerIndex(j + 1); + Index bId = other.m_outerIndex(j); + Index bEnd = other.m_outerIndex(j + 1); + while (aId < aEnd && bId < bEnd) { + Index aInner = m_innerIndex(aId); + Index bInner = other.m_innerIndex(bId); + if (aInner < bInner) { + ++aId; + } else if (bInner < aInner) { + ++bId; + } else { + result.m_innerIndex(nnz) = StorageIndex_(aInner); + result.blockRef(nnz) = func(blockRef(aId++), other.blockRef(bId++)); + ++nnz; + } + } + } + result.m_outerIndex(m_blockOuterSize) = StorageIndex_(nnz); + result.conservativeResizeBlockStorage_(nnz); + return result; + } + + template + BlockSparseMatrix transposeImpl() const; + + // Resize both block storage arrays non-conservatively and update the capacity counter. + void resizeBlockStorage_(Index n) { + m_innerIndex.resize(n); + m_values.resize(n * BlockSize); + } + + void conservativeResizeBlockStorage_(Index n) { + m_innerIndex.conservativeResize(n); + m_values.conservativeResize(n * BlockSize); + } + + // ------------------------------------------------------------------------- + // Storage + // ------------------------------------------------------------------------- + Index m_blockOuterSize = 0; // block-cols (ColMajor) or block-rows (RowMajor) + Index m_blockInnerSize = 0; // block-rows (ColMajor) or block-cols (RowMajor) + + Array m_outerIndex = // size: m_blockOuterSize + 1 + decltype(m_outerIndex)::Zero(m_blockOuterSize + 1); + Array m_innerIndex; + // Block values stored consecutively; block k occupies + // m_values[k*BlockSize .. (k+1)*BlockSize - 1]. + Array m_values; + + // ------------------------------------------------------------------------- + // MultiInnerIterator + // + // Simultaneously walks BlockOuterSize_ consecutive outer vectors of a + // compressed SparseMatrix, yielding scalar entries one at a time in + // non-decreasing inner-index order. On each step the sub-iterator with + // the smallest current inner index is the "active" one. + // + // BlockOuterSize_ is BlockCols_ for ColMajor (one block-column at a time) + // and BlockRows_ for RowMajor (one block-row at a time). + // ------------------------------------------------------------------------- + template + class MultiInnerIterator { + using StorageIndex = typename SparseMatrixType::StorageIndex; + static constexpr int BlockOuterSize_ = IsRowMajor ? BlockRows_ : BlockCols_; + + public: + MultiInnerIterator(const SparseMatrixType& mat, Index outerBase) + : m_outerPtr(mat.outerIndexPtr()), + m_innerPtr(mat.innerIndexPtr()), + m_valuePtr(mat.valuePtr()), + m_outerBase(outerBase) { + for (int k = 0; k < BlockOuterSize_; ++k) m_pos[k] = m_outerPtr[outerBase + k]; + advance(); + } + + EIGEN_STRONG_INLINE operator bool() const { return m_valid; } + + EIGEN_STRONG_INLINE MultiInnerIterator& operator++() { + ++m_pos[m_active]; + advance(); + return *this; + } + + // Absolute outer index of the current entry. + EIGEN_STRONG_INLINE Index outer() const { return m_outerBase + m_active; } + // Inner index of the current entry. + EIGEN_STRONG_INLINE StorageIndex index() const { return m_innerPtr[m_pos[m_active]]; } + // Scalar value of the current entry. + EIGEN_STRONG_INLINE Scalar value() const { return m_valuePtr[m_pos[m_active]]; } + + private: + // Find the sub-iterator with the smallest inner index. + void advance() { + m_valid = false; + for (int k = 0; k < BlockOuterSize_; ++k) { + if (m_pos[k] < m_outerPtr[m_outerBase + k + 1]) { + if (!m_valid || m_innerPtr[m_pos[k]] < m_innerPtr[m_pos[m_active]]) { + m_active = k; + m_valid = true; + } + } + } + } + + const StorageIndex* m_outerPtr; + const StorageIndex* m_innerPtr; + const Scalar* m_valuePtr; + Index m_outerBase; + StorageIndex m_pos[BlockOuterSize_]; + int m_active = 0; + bool m_valid = false; + }; + + // Allow other instantiations of BlockSparseMatrix to access private members + // (needed by operator*). + template + friend class BlockSparseMatrix; + template + friend class BlockSparseTriangularView; + template + friend class BlockSparseSelfAdjointView; +}; + +// ============================================================================= +// Out-of-line method definitions +// ============================================================================= + +// ----------------------------------------------------------------------------- +// setFromTriplets +// ----------------------------------------------------------------------------- + +template +template +void BlockSparseMatrix::setFromTriplets(InputIterator begin, + InputIterator end) { + Index n = static_cast(std::distance(begin, end)); + + // Copy triplet coordinates and block values into Eigen arrays. + Array tOuter(n), tInner(n); + Array tValues(n * BlockSize); + + Index k = 0; + for (InputIterator it = begin; it != end; ++it, ++k) { + eigen_assert(it->row() >= 0 && it->row() < blockRows() && "setFromTriplets: block row out of range"); + eigen_assert(it->col() >= 0 && it->col() < blockCols() && "setFromTriplets: block col out of range"); + tOuter(k) = IsRowMajor ? StorageIndex(it->row()) : StorageIndex(it->col()); + tInner(k) = IsRowMajor ? StorageIndex(it->col()) : StorageIndex(it->row()); + BlockMap(tValues.data() + k * BlockSize) = it->value(); + } + + // Order the triplet indices by (outer, inner) using a stable LSD radix sort: + // pass 1 buckets by inner, pass 2 by outer. This runs in + // O(n + blockInnerSize + blockOuterSize) with sequential writes, avoiding the + // cache-unfriendly indirect comparison sort. + Array order(n), scratch(n); + { + Array count = Array::Zero(m_blockInnerSize + 1); + for (Index i = 0; i < n; ++i) count(tInner(i) + 1)++; + for (Index i = 0; i < m_blockInnerSize; ++i) count(i + 1) += count(i); + for (Index i = 0; i < n; ++i) order(count(tInner(i))++) = i; + } + { + Array count = Array::Zero(m_blockOuterSize + 1); + for (Index i = 0; i < n; ++i) count(tOuter(i) + 1)++; + for (Index i = 0; i < m_blockOuterSize; ++i) count(i + 1) += count(i); + for (Index i = 0; i < n; ++i) { + Index idx = order(i); + scratch(count(tOuter(idx))++) = idx; + } + order.swap(scratch); + } + + // Reset and pre-allocate (worst case: all n triplets are distinct blocks). + m_outerIndex.resize(m_blockOuterSize + 1); + m_outerIndex.setZero(); + resizeBlockStorage_(n); + + Index nnz = 0; + k = 0; + while (k < n) { + Index pi = order(k); + StorageIndex outer = tOuter(pi); + StorageIndex inner = tInner(pi); + + BlockType block = ConstBlockMap(tValues.data() + pi * BlockSize); + ++k; + + // Accumulate duplicate entries at the same (outer, inner) position. + while (k < n) { + Index pk = order(k); + if (tOuter(pk) != outer || tInner(pk) != inner) break; + block += ConstBlockMap(tValues.data() + pk * BlockSize); + ++k; + } + + m_innerIndex(nnz) = inner; + blockRef(nnz) = block; + m_outerIndex(outer + 1)++; + ++nnz; + } + + // Trim to actual number of unique blocks. + conservativeResizeBlockStorage_(nnz); + + // Convert per-outer block counts to prefix sums. + for (Index j = 0; j < m_blockOuterSize; ++j) { + m_outerIndex(j + 1) += m_outerIndex(j); + } +} + +// ----------------------------------------------------------------------------- +// toSparse +// ----------------------------------------------------------------------------- + +template +SparseMatrix +BlockSparseMatrix::toSparse() const { + SparseMatrix result(rows(), cols()); + result.reserve(nonZeroBlocks() * BlockSize); + + if (!IsRowMajor) { + // ColMajor: outer = block-column j. Emit scalar columns j*BlockCols+c + // in order c = 0..BlockCols-1. Within each scalar column, blocks are + // sorted by bi (block-row), so scalar rows bi*BlockRows+r are increasing. + for (Index j = 0; j < m_blockOuterSize; ++j) { + for (Index c = 0; c < BlockCols_; ++c) { + result.startVec(j * BlockCols_ + c); + for (Index id = m_outerIndex(j); id < m_outerIndex(j + 1); ++id) { + Index bi = m_innerIndex(id); + ConstBlockMap blk = blockRef(id); + for (Index r = 0; r < BlockRows_; ++r) { + result.insertBack(bi * BlockRows_ + r, j * BlockCols_ + c) = blk(r, c); + } + } + } + } + } else { + // RowMajor: outer = block-row bi. Emit scalar rows bi*BlockRows+r + // in order r = 0..BlockRows-1. Within each scalar row, blocks are + // sorted by j (block-col), so scalar cols j*BlockCols+c are increasing. + for (Index bi = 0; bi < m_blockOuterSize; ++bi) { + for (Index r = 0; r < BlockRows_; ++r) { + result.startVec(bi * BlockRows_ + r); + for (Index id = m_outerIndex(bi); id < m_outerIndex(bi + 1); ++id) { + Index j = m_innerIndex(id); + ConstBlockMap blk = blockRef(id); + for (Index c = 0; c < BlockCols_; ++c) { + result.insertBack(bi * BlockRows_ + r, j * BlockCols_ + c) = blk(r, c); + } + } + } + } + } + + result.finalize(); + return result; +} + +// ----------------------------------------------------------------------------- +// fromSparse +// ----------------------------------------------------------------------------- + +template +BlockSparseMatrix +BlockSparseMatrix::fromSparse( + const SparseMatrix& sp) { + eigen_assert(sp.rows() % BlockRows_ == 0 && "matrix rows not divisible by BlockRows"); + eigen_assert(sp.cols() % BlockCols_ == 0 && "matrix cols not divisible by BlockCols"); + eigen_assert(sp.isCompressed() && "fromSparse requires a compressed SparseMatrix"); + + Index bRows = sp.rows() / BlockRows_; + Index bCols = sp.cols() / BlockCols_; + + // BlockOuterSize: how many consecutive outer vectors form one block-outer strip. + // BlockInnerSize: the inner dimension of each block. + constexpr Index BlockOuterSize = IsRowMajor ? BlockRows_ : BlockCols_; + constexpr Index BlockInnerSize = IsRowMajor ? BlockCols_ : BlockRows_; + constexpr StorageIndex_ kEmptyIndex = -1; + + using SpMat = SparseMatrix; + + BlockSparseMatrix result(bRows, bCols); + + // Pass 1: count the number of unique block-inner indices per block-outer, + // by scanning each group of BlockOuterSize consecutive outer vectors together. + for (Index outerBlock = 0; outerBlock < result.m_blockOuterSize; ++outerBlock) { + StorageIndex_ prevInnerBlock = kEmptyIndex; + for (MultiInnerIterator it(sp, outerBlock * BlockOuterSize); it; ++it) { + StorageIndex_ innerBlock = it.index() / StorageIndex_(BlockInnerSize); + if (innerBlock != prevInnerBlock) { + result.m_outerIndex(outerBlock + 1)++; + prevInnerBlock = innerBlock; + } + } + } + + // Prefix sum → result.m_outerIndex becomes the standard CSC/CSR outer pointer. + for (Index j = 0; j < result.m_blockOuterSize; ++j) result.m_outerIndex(j + 1) += result.m_outerIndex(j); + + Index nBlocks = result.m_outerIndex(result.m_blockOuterSize); + result.resizeBlockStorage_(nBlocks); + result.m_values.setZero(); + + // Pass 2: scatter each scalar entry directly into its position within the + // pre-zeroed block value array. + for (Index outerBlock = 0; outerBlock < result.m_blockOuterSize; ++outerBlock) { + Index blockId = result.m_outerIndex(outerBlock) - 1; // incremented on first new block + StorageIndex_ prevInnerBlock = kEmptyIndex; + + for (MultiInnerIterator it(sp, outerBlock * BlockOuterSize); it; ++it) { + Index absOuter = it.outer(); // absolute outer index in sp + StorageIndex_ innerIdx = it.index(); // inner index in sp + StorageIndex_ innerBlock = innerIdx / StorageIndex_(BlockInnerSize); + + if (innerBlock != prevInnerBlock) { + ++blockId; + result.m_innerIndex(blockId) = innerBlock; + prevInnerBlock = innerBlock; + } + + // Scatter into block storage (layout matches the BSM's Options_). + // ColMajor blocks: col * BlockRows_ + row = localOuter * BlockInnerSize + localInner + // RowMajor blocks: row * BlockCols_ + col = localOuter * BlockInnerSize + localInner + Index localOuter = absOuter % BlockOuterSize; + Index localInner = innerIdx % BlockInnerSize; + Index offset = localOuter * BlockInnerSize + localInner; + + result.m_values(blockId * BlockSize + offset) = it.value(); + } + } + + return result; +} + +// ----------------------------------------------------------------------------- +// operator* (block-sparse product) +// ----------------------------------------------------------------------------- + +template +template +BlockSparseMatrix +BlockSparseMatrix::operator*( + const BlockSparseMatrix& rhs) const { + using RhsMatrix = BlockSparseMatrix; + using ResultMatrix = BlockSparseMatrix; + using ResultBlock = Matrix; + constexpr int ResultBlockSize = BlockRows_ * RhsBlockCols; + + eigen_assert(blockCols() == rhs.blockRows() && "BlockSparseMatrix product: lhs.blockCols() != rhs.blockRows()"); + + Index cBlockRows = blockRows(); + Index cBlockCols = rhs.blockCols(); + ResultMatrix result(cBlockRows, cBlockCols); + + // For ColMajor: mask / accum indexed by block-row (size = cBlockRows). + // For RowMajor: mask / accum indexed by block-col (size = cBlockCols). + Index maskSize = IsRowMajor ? cBlockCols : cBlockRows; + Array mask = Array::Zero(maskSize); + Array accumData(maskSize * ResultBlockSize); + Array indices(maskSize); + Index nIndices = 0; + + // Grow result storage geometrically rather than pre-allocating the dense + // worst case (cBlockRows*cBlockCols blocks): the product is typically far + // sparser than that, so the dense bound would blow up peak memory. + // capacity is always kept <= maxResultNnz (the true upper bound), and since + // any single outer emits at most maskSize blocks, the initial estimate of + // maskSize guarantees the first outer fits before the first grow check. + Index cOuterSize = result.m_blockOuterSize; + Index maxResultNnz = cBlockRows * cBlockCols; + Index capacity = numext::mini(maxResultNnz, numext::maxi(maskSize, nonZeroBlocks() + rhs.nonZeroBlocks())); + result.resizeBlockStorage_(capacity); + Index nnz = 0; + + for (Index out = 0; out < cOuterSize; ++out) { + result.m_outerIndex(out) = StorageIndex_(nnz); + + if (!IsRowMajor) { + // ColMajor: out is block-column j of the result. + // For each block B(k,j) and each block A(bi,k): C(bi,j) += A(bi,k)*B(k,j). + Index j = out; + for (Index rhsId = rhs.m_outerIndex(j); rhsId < rhs.m_outerIndex(j + 1); ++rhsId) { + Index k = rhs.m_innerIndex(rhsId); + typename RhsMatrix::ConstBlockMap Bkj = rhs.blockRef(rhsId); + for (Index lhsId = m_outerIndex(k); lhsId < m_outerIndex(k + 1); ++lhsId) { + Index bi = m_innerIndex(lhsId); + if (!mask(bi)) { + mask(bi) = 1; + Map(accumData.data() + bi * ResultBlockSize).noalias() = blockRef(lhsId) * Bkj; + indices(nIndices++) = bi; + } else { + Map(accumData.data() + bi * ResultBlockSize).noalias() += blockRef(lhsId) * Bkj; + } + } + } + } else { + // RowMajor: out is block-row bi of the result. + // For each block A(bi,k) and each block B(k,j): C(bi,j) += A(bi,k)*B(k,j). + Index bi = out; + for (Index lhsId = m_outerIndex(bi); lhsId < m_outerIndex(bi + 1); ++lhsId) { + Index k = m_innerIndex(lhsId); + ConstBlockMap Aik = blockRef(lhsId); + for (Index rhsId = rhs.m_outerIndex(k); rhsId < rhs.m_outerIndex(k + 1); ++rhsId) { + Index j = rhs.m_innerIndex(rhsId); + if (!mask(j)) { + mask(j) = 1; + Map(accumData.data() + j * ResultBlockSize).noalias() = Aik * rhs.blockRef(rhsId); + indices(nIndices++) = j; + } else { + Map(accumData.data() + j * ResultBlockSize).noalias() += Aik * rhs.blockRef(rhsId); + } + } + } + } + + // Sort the accumulated indices so the result's inner index array is sorted. + std::sort(indices.data(), indices.data() + nIndices); + if (nnz + nIndices > capacity) { + capacity = numext::mini(maxResultNnz, numext::maxi(2 * capacity, nnz + nIndices)); + result.conservativeResizeBlockStorage_(capacity); + } + for (Index ki = 0; ki < nIndices; ++ki) { + Index idx = indices(ki); + result.m_innerIndex(nnz) = StorageIndex_(idx); + result.blockRef(nnz) = Map(accumData.data() + idx * ResultBlockSize); + mask(idx) = 0; + ++nnz; + } + nIndices = 0; + } + result.m_outerIndex(cOuterSize) = StorageIndex_(nnz); + + // Trim to actual number of result blocks. + result.conservativeResizeBlockStorage_(nnz); + + return result; +} + +// ============================================================================= +// BlockSparseTriangularView +// ============================================================================= + +/** \class BlockSparseTriangularView + * \ingroup SparseCore_Module + * \brief Lazy block-level triangular view of a BlockSparseMatrix. + * + * Obtained via \c BSM::triangularView() or + * \c BSM::triangularView() (the latter asserts that the unused + * triangle of every diagonal block is already zero, enabling faster vectorised + * products without an explicit triangularView on the block). + * + * By default (\p DiagIsTriangular = \c false) diagonal blocks are treated as + * triangular regardless of what is stored in the unused half: products call + * \c block.triangularView() and \c eval() zeroes the unused + * triangle. Set \p DiagIsTriangular = \c true to skip that overhead when you + * can guarantee the unused triangle is already zero. + */ +template +class BlockSparseTriangularView { + public: + using Scalar = typename BSM::Scalar; + using StorageIndex = typename BSM::StorageIndex; + using BlockType = typename BSM::BlockType; + using BlockMap = typename BSM::BlockMap; + using ConstBlockMap = typename BSM::ConstBlockMap; + static constexpr int BlockRows = BSM::BlockRows; + static constexpr int BlockCols = BSM::BlockCols; + static constexpr int BlockSize = BSM::BlockSize; + static constexpr bool IsRowMajor = BSM::IsRowMajor; + static constexpr bool IsUpper = (Mode & Upper) != 0; + + explicit BlockSparseTriangularView(const BSM& m) : m_matrix(m) {} + + Index rows() const { return m_matrix.rows(); } + Index cols() const { return m_matrix.cols(); } + + // ---- Materialize --------------------------------------------------------- + + /** Copy the triangular blocks into a new BSM; off-triangle blocks are dropped. + * When DiagIsTriangular is false the unused triangle of each diagonal block + * is explicitly zeroed in the output. */ + BSM eval() const { + constexpr int ZeroMode = IsUpper ? StrictlyLower : StrictlyUpper; + const BSM& m = m_matrix; + BSM result(m.blockRows(), m.blockCols()); + result.resizeBlockStorage_(m.nonZeroBlocks()); + Index nnz = 0; + + for (Index out = 0; out < m.m_blockOuterSize; ++out) { + result.m_outerIndex(out) = StorageIndex(nnz); + for (Index id = m.m_outerIndex(out); id < m.m_outerIndex(out + 1); ++id) { + Index inner = m.m_innerIndex(id); + Index bi = IsRowMajor ? out : inner; + Index bj = IsRowMajor ? inner : out; + if (IsUpper ? (bj < bi) : (bj > bi)) continue; + result.m_innerIndex(nnz) = StorageIndex(inner); + result.m_values.template segment(nnz * BlockSize) = + m.m_values.template segment(id * BlockSize); + EIGEN_IF_CONSTEXPR (!DiagIsTriangular) { + if (bi == bj) + BlockMap(result.m_values.data() + nnz * BlockSize).template triangularView().setZero(); + } + ++nnz; + } + } + result.m_outerIndex(m.m_blockOuterSize) = StorageIndex(nnz); + result.conservativeResizeBlockStorage_(nnz); + return result; + } + + /** Convert to a scalar-level SparseMatrix (off-triangle blocks zeroed). */ + SparseMatrix toSparse() const { return eval().toSparse(); } + + // ---- Arithmetic ---------------------------------------------------------- + + BSM operator+(const BlockSparseTriangularView& other) const { return eval() + other.eval(); } + BSM operator-(const BlockSparseTriangularView& other) const { return eval() - other.eval(); } + + /** Tri × BSM product (materialises this view then delegates). */ + template + BlockSparseMatrix operator*( + const BlockSparseMatrix& rhs) const { + return eval() * rhs; + } + + // ---- Dense products (no intermediate materialisation) -------------------- + + template + Matrix operator*(const MatrixBase& rhs) const { + EIGEN_STATIC_ASSERT( + (std::is_same::value), + YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) + eigen_assert(m_matrix.cols() == rhs.rows() && "BlockSparseTriangularView * Dense: dimension mismatch"); + using ResultType = Matrix; + ResultType result = ResultType::Zero(m_matrix.rows(), rhs.cols()); + for (Index out = 0; out < m_matrix.m_blockOuterSize; ++out) { + for (Index id = m_matrix.m_outerIndex(out); id < m_matrix.m_outerIndex(out + 1); ++id) { + Index inner = m_matrix.m_innerIndex(id); + Index bi = IsRowMajor ? out : inner; + Index bj = IsRowMajor ? inner : out; + if (IsUpper ? (bj < bi) : (bj > bi)) continue; + constexpr int DiagMode = IsUpper ? Upper : Lower; + if (!DiagIsTriangular && bi == bj) + result.template middleRows(bi * BlockRows).noalias() += + m_matrix.blockRef(id).template triangularView() * + rhs.template middleRows(bj * BlockCols); + else + result.template middleRows(bi * BlockRows).noalias() += + m_matrix.blockRef(id) * rhs.template middleRows(bj * BlockCols); + } + } + return result; + } + + template + friend Matrix operator*(const MatrixBase& lhs, + const BlockSparseTriangularView& tri) { + EIGEN_STATIC_ASSERT( + (std::is_same::value), + YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) + eigen_assert(lhs.cols() == tri.m_matrix.rows() && "Dense * BlockSparseTriangularView: dimension mismatch"); + constexpr bool isRM = BSM::IsRowMajor; + using ResultType = Matrix; + ResultType result = ResultType::Zero(lhs.rows(), tri.m_matrix.cols()); + for (Index out = 0; out < tri.m_matrix.m_blockOuterSize; ++out) { + for (Index id = tri.m_matrix.m_outerIndex(out); id < tri.m_matrix.m_outerIndex(out + 1); ++id) { + Index inner = tri.m_matrix.m_innerIndex(id); + Index bi = isRM ? out : inner; + Index bj = isRM ? inner : out; + if (IsUpper ? (bj < bi) : (bj > bi)) continue; + constexpr int DiagMode = IsUpper ? Upper : Lower; + if (!DiagIsTriangular && bi == bj) + result.template middleCols(bj * BlockCols).noalias() += + lhs.template middleCols(bi * BlockRows) * + tri.m_matrix.blockRef(id).template triangularView(); + else + result.template middleCols(bj * BlockCols).noalias() += + lhs.template middleCols(bi * BlockRows) * tri.m_matrix.blockRef(id); + } + } + return result; + } + + // ---- Triangular solve ------------------------------------------------------- + + /** Solve T * x = rhs in-place. Requires square blocks. + * ColMajor Lower: forward sub, diagonal first per column. + * ColMajor Upper: backward sub, diagonal last per column. + * RowMajor Lower: forward sub, diagonal last per row. + * RowMajor Upper: backward sub, diagonal first per row. + */ + template + void solveInPlace(MatrixBase& x) const { + doSolveImpl(x.derived()); + } + + /** Proxy returned by transpose(): solveInPlace solves T^T x = b. */ + struct TransposeReturnType { + const BlockSparseTriangularView& m_tri; + template + void solveInPlace(MatrixBase& x) const { + m_tri.template doSolveImpl(x.derived()); + } + }; + + /** Proxy returned by adjoint(): solveInPlace solves T^H x = b. */ + struct AdjointReturnType { + const BlockSparseTriangularView& m_tri; + template + void solveInPlace(MatrixBase& x) const { + m_tri.template doSolveImpl(x.derived()); + } + }; + + TransposeReturnType transpose() const { return {*this}; } + AdjointReturnType adjoint() const { return {*this}; } + + private: + const BSM& m_matrix; + + // Non-transposed solve for both storage orders. + // + // diagFirst = (IsUpper == IsRowMajor): ColMajor Lower→first, ColMajor Upper→last, + // RowMajor Lower→last, RowMajor Upper→first. + // Loop direction: forward for Lower, backward for Upper (same for both storage orders). + // ColMajor: solve diagonal first, then scatter x[inner] -= blk * x[k]. + // RowMajor: gather x[k] -= blk * x[inner] first, then solve diagonal. + template + void doSolveDirect(Derived& x) const { + EIGEN_STATIC_ASSERT(BlockRows == BlockCols, THIS_METHOD_IS_ONLY_FOR_SQUARE_BLOCK_MATRICES) + constexpr int DiagMode = IsUpper ? Upper : Lower; + constexpr bool diagFirst = (IsUpper == BSM::IsRowMajor); + Index nb = m_matrix.blockCols(); // == blockRows() for square matrices + eigen_assert(x.rows() == m_matrix.rows() && "solveInPlace: size mismatch"); + + const StorageIndex* innerPtr = m_matrix.innerIndexPtr(); + const StorageIndex* outerPtr = m_matrix.outerIndexPtr(); + + Index outerStart = IsUpper ? nb - 1 : 0; + Index outerEnd = IsUpper ? -1 : nb; + constexpr Index kStep = IsUpper ? -1 : 1; + + for (Index k = outerStart; k != outerEnd; k += kStep) { + const StorageIndex* beg = innerPtr + outerPtr[k]; + const StorageIndex* end = innerPtr + outerPtr[k + 1]; + if (beg == end) continue; + const StorageIndex* diag_ptr = diagFirst ? beg : end - 1; + const StorageIndex* off_beg = diagFirst ? beg + 1 : beg; + const StorageIndex* off_end = diagFirst ? end : end - 1; + eigen_assert(*diag_ptr == k); + EIGEN_IF_CONSTEXPR (!BSM::IsRowMajor) { + m_matrix.blockRef(diag_ptr - innerPtr) + .template triangularView() + .solveInPlace(x.template middleRows(k * BlockRows)); + for (const StorageIndex* it = off_beg; it != off_end; ++it) + x.template middleRows(*it * BlockRows).noalias() -= + m_matrix.blockRef(it - innerPtr) * x.template middleRows(k * BlockRows); + } else { + for (const StorageIndex* it = off_beg; it != off_end; ++it) + x.template middleRows(k * BlockRows).noalias() -= + m_matrix.blockRef(it - innerPtr) * x.template middleRows(*it * BlockRows); + m_matrix.blockRef(diag_ptr - innerPtr) + .template triangularView() + .solveInPlace(x.template middleRows(k * BlockRows)); + } + } + } + + // Transposed/adjoint solve for both storage orders. + // + // Loop direction: forward for Upper, backward for Lower (same for both storage orders). + // ColMajor: gather x[k] -= adj(blk) * x[inner] first, then solve adj(diagonal). + // RowMajor: solve adj(diagonal) first, then scatter x[inner] -= adj(blk) * x[k]. + template + void doSolveTransposed(Derived& x) const { + EIGEN_STATIC_ASSERT(BlockRows == BlockCols, THIS_METHOD_IS_ONLY_FOR_SQUARE_BLOCK_MATRICES) + constexpr int DiagMode = IsUpper ? Upper : Lower; + constexpr bool diagFirst = (IsUpper == BSM::IsRowMajor); + Index nb = m_matrix.blockCols(); // == blockRows() for square matrices + eigen_assert(x.rows() == m_matrix.rows() && "solveInPlace: size mismatch"); + + const StorageIndex* innerPtr = m_matrix.innerIndexPtr(); + const StorageIndex* outerPtr = m_matrix.outerIndexPtr(); + + Index outerStart = IsUpper ? 0 : nb - 1; + Index outerEnd = IsUpper ? nb : -1; + constexpr Index kStep = IsUpper ? 1 : -1; + + for (Index k = outerStart; k != outerEnd; k += kStep) { + const StorageIndex* beg = innerPtr + outerPtr[k]; + const StorageIndex* end = innerPtr + outerPtr[k + 1]; + if (beg == end) continue; + const StorageIndex* diag_ptr = diagFirst ? beg : end - 1; + const StorageIndex* off_beg = diagFirst ? beg + 1 : beg; + const StorageIndex* off_end = diagFirst ? end : end - 1; + eigen_assert(*diag_ptr == k); + EIGEN_IF_CONSTEXPR (!BSM::IsRowMajor) { + for (const StorageIndex* it = off_beg; it != off_end; ++it) + x.template middleRows(k * BlockRows).noalias() -= + internal::adjoint_if(m_matrix.blockRef(it - innerPtr)) * + x.template middleRows(*it * BlockRows); + internal::adjoint_if(m_matrix.blockRef(diag_ptr - innerPtr).template triangularView()) + .solveInPlace(x.template middleRows(k * BlockRows)); + } else { + internal::adjoint_if(m_matrix.blockRef(diag_ptr - innerPtr).template triangularView()) + .solveInPlace(x.template middleRows(k * BlockRows)); + for (const StorageIndex* it = off_beg; it != off_end; ++it) + x.template middleRows(*it * BlockRows).noalias() -= + internal::adjoint_if(m_matrix.blockRef(it - innerPtr)) * + x.template middleRows(k * BlockRows); + } + } + } + + template + void doSolveImpl(Derived& x) const { + EIGEN_IF_CONSTEXPR (!Transposed) + doSolveDirect(x); + else + doSolveTransposed(x); + } +}; + +// ============================================================================= +// BlockSparseSelfAdjointView +// ============================================================================= + +/** \class BlockSparseSelfAdjointView + * \ingroup SparseCore_Module + * \brief Lazy block-level self-adjoint (Hermitian) view of a BlockSparseMatrix. + * + * Obtained via \c BSM::selfadjointView() or + * \c BSM::selfadjointView() (the latter signals that every + * diagonal block is itself Hermitian, enabling DSYMM/ZHEMM on those blocks). + * + * Only the triangle selected by \p UpLo is read; the other triangle is + * reconstructed on-the-fly as the adjoint of each stored off-diagonal block. + * + * \pre BSM::BlockRows == BSM::BlockCols. + */ +template +class BlockSparseSelfAdjointView { + public: + using Scalar = typename BSM::Scalar; + using StorageIndex = typename BSM::StorageIndex; + using BlockType = typename BSM::BlockType; + using BlockMap = typename BSM::BlockMap; + using ConstBlockMap = typename BSM::ConstBlockMap; + static constexpr int BlockRows = BSM::BlockRows; // == BlockCols + static constexpr int BlockCols = BSM::BlockCols; + static constexpr int BlockSize = BSM::BlockSize; + static constexpr bool IsRowMajor = BSM::IsRowMajor; + static constexpr bool IsUpper = (UpLo & Upper) != 0; + // UpLo passed to Eigen's dense selfadjointView on diagonal blocks: + static constexpr int DiagUpLo = IsUpper ? Upper : Lower; + + explicit BlockSparseSelfAdjointView(const BSM& m) : m_matrix(m) {} + + Index rows() const { return m_matrix.rows(); } + Index cols() const { return m_matrix.cols(); } + + // ---- Materialize --------------------------------------------------------- + + /** Build a full symmetric BSM: stored triangle + adjoint mirror of each + * off-diagonal block. Diagonal blocks: when DiagIsSelfAdjoint is false, + * both triangles are filled from the stored triangle via selfadjointView; + * when true the block is already fully populated and is copied as-is. */ + BSM eval() const { + const BSM& m = m_matrix; + + Index nDiag = 0, nOff = 0; + for (Index out = 0; out < m.m_blockOuterSize; ++out) + for (Index id = m.m_outerIndex(out); id < m.m_outerIndex(out + 1); ++id) { + Index inner = m.m_innerIndex(id); + Index bi = IsRowMajor ? out : inner; + Index bj = IsRowMajor ? inner : out; + if (IsUpper ? (bj < bi) : (bj > bi)) continue; + if (bi == bj) + ++nDiag; + else + ++nOff; + } + + Index nTotal = nDiag + 2 * nOff; + + Array brows(nTotal), bcols(nTotal); + Array bvals(nTotal * BlockSize); + + Index k = 0; + for (Index out = 0; out < m.m_blockOuterSize; ++out) + for (Index id = m.m_outerIndex(out); id < m.m_outerIndex(out + 1); ++id) { + Index inner = m.m_innerIndex(id); + Index bi = IsRowMajor ? out : inner; + Index bj = IsRowMajor ? inner : out; + if (IsUpper ? (bj < bi) : (bj > bi)) continue; + + brows(k) = StorageIndex(bi); + bcols(k) = StorageIndex(bj); + if (!DiagIsSelfAdjoint && bi == bj) + BlockMap(bvals.data() + k * BlockSize) = m.blockRef(id).template selfadjointView(); + else + BlockMap(bvals.data() + k * BlockSize) = m.blockRef(id); + ++k; + + if (bi != bj) { + brows(k) = StorageIndex(bj); + bcols(k) = StorageIndex(bi); + BlockMap(bvals.data() + k * BlockSize) = m.blockRef(id).adjoint(); + ++k; + } + } + + // Sort by (outer, inner) then build BSM directly (no duplicates by construction). + Array perm(nTotal); + std::iota(perm.data(), perm.data() + nTotal, Index(0)); + std::sort(perm.data(), perm.data() + nTotal, [&](Index a, Index b) { + StorageIndex ao = IsRowMajor ? brows(a) : bcols(a); + StorageIndex bo = IsRowMajor ? brows(b) : bcols(b); + if (ao != bo) return ao < bo; + return (IsRowMajor ? bcols(a) : brows(a)) < (IsRowMajor ? bcols(b) : brows(b)); + }); + + BSM result(m.blockRows(), m.blockCols()); + result.resizeBlockStorage_(nTotal); + + for (Index ki = 0; ki < nTotal; ++ki) { + Index pi = perm(ki); + StorageIndex outer = IsRowMajor ? brows(pi) : bcols(pi); + StorageIndex inner = IsRowMajor ? bcols(pi) : brows(pi); + result.m_outerIndex(outer + 1)++; + result.m_innerIndex(ki) = inner; + result.m_values.template segment(ki * BlockSize) = bvals.template segment(pi * BlockSize); + } + for (Index j = 0; j < result.m_blockOuterSize; ++j) result.m_outerIndex(j + 1) += result.m_outerIndex(j); + + return result; + } + + /** Convert to a symmetrised scalar-level SparseMatrix. */ + SparseMatrix toSparse() const { return eval().toSparse(); } + + // ---- Arithmetic ---------------------------------------------------------- + + BSM operator+(const BlockSparseSelfAdjointView& other) const { return eval() + other.eval(); } + BSM operator-(const BlockSparseSelfAdjointView& other) const { return eval() - other.eval(); } + + /** SelfAdj × BSM: materialises the view then uses the general SpGEMM. */ + template + BlockSparseMatrix operator*( + const BlockSparseMatrix& rhs) const { + return eval() * rhs; + } + + // ---- Dense products (no materialisation; exploits both triangles) --------- + + /** SelfAdj × Dense. + * + * Off-diagonal stored block A(bi,bj) contributes: + * result(bi) += A(bi,bj) * rhs(bj) [stored triangle] + * result(bj) += A(bi,bj)^H * rhs(bi) [implicit mirror] + * + * When DiagIsSelfAdjoint is true, both triangles of each diagonal block + * are valid; a plain product is used (faster for small fixed-size blocks). + * Otherwise only the UpLo triangle is assumed valid and selfadjointView + * reconstructs the full diagonal-block product. + */ + template + Matrix operator*(const MatrixBase& rhs) const { + EIGEN_STATIC_ASSERT( + (std::is_same::value), + YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) + eigen_assert(m_matrix.cols() == rhs.rows() && "BlockSparseSelfAdjointView * Dense: dimension mismatch"); + using ResultType = Matrix; + ResultType result = ResultType::Zero(m_matrix.rows(), rhs.cols()); + + for (Index out = 0; out < m_matrix.m_blockOuterSize; ++out) + for (Index id = m_matrix.m_outerIndex(out); id < m_matrix.m_outerIndex(out + 1); ++id) { + Index inner = m_matrix.m_innerIndex(id); + Index bi = IsRowMajor ? out : inner; + Index bj = IsRowMajor ? inner : out; + if (IsUpper ? (bj < bi) : (bj > bi)) continue; + + if (bi == bj) { + EIGEN_IF_CONSTEXPR (DiagIsSelfAdjoint) { + result.template middleRows(bi * BlockRows).noalias() += + m_matrix.blockRef(id) * rhs.template middleRows(bj * BlockCols); + } else { + // Materialize the tiny diagonal block as a fixed-size Hermitian matrix, then use the + // ordinary (coeff-based for a vector rhs) product. This avoids the runtime-sized, + // EIGEN_DONT_INLINE selfadjoint_matrix_vector_product kernel, which is tuned for large + // matrices and is pure overhead for a 2-4 row block. + BlockType diag = m_matrix.blockRef(id).template selfadjointView(); + result.template middleRows(bi * BlockRows).noalias() += + diag * rhs.template middleRows(bj * BlockCols); + } + } else { + result.template middleRows(bi * BlockRows).noalias() += + m_matrix.blockRef(id) * rhs.template middleRows(bj * BlockCols); + result.template middleRows(bj * BlockRows).noalias() += + m_matrix.blockRef(id).adjoint() * rhs.template middleRows(bi * BlockRows); + } + } + return result; + } + + /** Dense × SelfAdj: lhs * A == (A^H * lhs^H)^H == (A * lhs^H)^H for Hermitian A. */ + template + friend Matrix operator*(const MatrixBase& lhs, + const BlockSparseSelfAdjointView& view) { + EIGEN_STATIC_ASSERT( + (std::is_same::value), + YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) + return (view * lhs.adjoint()).adjoint(); + } + + private: + const BSM& m_matrix; +}; + +// ============================================================================= +// BlockSparseMatrix::transposeImpl / transpose / adjoint (out-of-line) +// ============================================================================= + +template +template +BlockSparseMatrix +BlockSparseMatrix::transposeImpl() const { + using ResultType = BlockSparseMatrix; + ResultType result(blockCols(), blockRows()); + + // Count entries per new outer (= old inner). + for (Index id = 0; id < nonZeroBlocks(); ++id) result.m_outerIndex(m_innerIndex(id) + 1)++; + + // Prefix sum. + for (Index j = 0; j < result.m_blockOuterSize; ++j) result.m_outerIndex(j + 1) += result.m_outerIndex(j); + + Index nnz = nonZeroBlocks(); + result.resizeBlockStorage_(nnz); + + // One insertion cursor per new outer; start at the prefix-sum boundary. + // Because we iterate oldOuter in increasing order, for each newOuter = oldInner + // the emitted newInner = oldOuter values are automatically sorted. + Array pos = result.m_outerIndex.head(result.m_blockOuterSize); + + for (Index oldOuter = 0; oldOuter < m_blockOuterSize; ++oldOuter) { + for (Index id = m_outerIndex(oldOuter); id < m_outerIndex(oldOuter + 1); ++id) { + Index newOuter = m_innerIndex(id); + Index insertAt = pos(newOuter)++; + result.m_innerIndex(insertAt) = StorageIndex_(oldOuter); + result.blockRef(insertAt) = internal::adjoint_if(blockRef(id)); + } + } + return result; +} + +template +BlockSparseMatrix +BlockSparseMatrix::transpose() const { + return transposeImpl(); +} + +template +BlockSparseMatrix +BlockSparseMatrix::adjoint() const { + return transposeImpl(); +} + +namespace internal { + +// --------------------------------------------------------------------------- +// generic_product_impl: BlockSparse × Dense → Dense +// Provides evalTo / addTo / subTo / scaleAndAddTo via generic_product_impl_base. +// --------------------------------------------------------------------------- +template +struct generic_product_impl + : generic_product_impl_base> { + using Scalar = typename Product::Scalar; + + template + static void scaleAndAddTo(Dst& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha) { + constexpr bool IsRM = (Lhs::Options & RowMajorBit) != 0; + constexpr int BR = Lhs::BlockRows; + constexpr int BC = Lhs::BlockCols; + const typename Lhs::StorageIndex* outerPtr = lhs.outerIndexPtr(); + const typename Lhs::StorageIndex* innerPtr = lhs.innerIndexPtr(); + // Branch on alpha before the loop: alpha==1 and alpha==-1 avoid creating a + // CwiseUnaryOp, which defeats SIMD for complex scalars. + bool a1 = (alpha == Scalar(1)); + bool am1 = (alpha == Scalar(-1)); + for (Eigen::Index out = 0; out < lhs.blockOuterSize(); ++out) { + for (Eigen::Index id = outerPtr[out]; id < outerPtr[out + 1]; ++id) { + Eigen::Index inner = innerPtr[id]; + Eigen::Index bi = IsRM ? out : inner; + Eigen::Index bj = IsRM ? inner : out; + auto dst_seg = dst.template middleRows
(bi * BR); + auto rhs_seg = rhs.template middleRows(bj * BC); + if (EIGEN_PREDICT_TRUE(a1)) + dst_seg.noalias() += lhs.blockRef(id) * rhs_seg; + else if (am1) + dst_seg.noalias() -= lhs.blockRef(id) * rhs_seg; + else { + // Materialize block×rhs_seg into a small fixed-size stack buffer, then + // scale by alpha. Keeps the B×B block as a plain Map for vectorization. + typedef Matrix TmpType; + TmpType tmp(BR, rhs.cols()); + tmp.noalias() = lhs.blockRef(id) * rhs_seg; + dst_seg.noalias() += alpha * tmp; + } + } + } + } +}; + +// --------------------------------------------------------------------------- +// generic_product_impl: Dense × BlockSparse → Dense +// --------------------------------------------------------------------------- +template +struct generic_product_impl + : generic_product_impl_base> { + using Scalar = typename Product::Scalar; + + template + static void scaleAndAddTo(Dst& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha) { + constexpr bool IsRM = (Rhs::Options & RowMajorBit) != 0; + constexpr int BR = Rhs::BlockRows; + constexpr int BC = Rhs::BlockCols; + const typename Rhs::StorageIndex* outerPtr = rhs.outerIndexPtr(); + const typename Rhs::StorageIndex* innerPtr = rhs.innerIndexPtr(); + bool a1 = (alpha == Scalar(1)); + bool am1 = (alpha == Scalar(-1)); + for (Eigen::Index out = 0; out < rhs.blockOuterSize(); ++out) { + for (Eigen::Index id = outerPtr[out]; id < outerPtr[out + 1]; ++id) { + Eigen::Index inner = innerPtr[id]; + Eigen::Index bi = IsRM ? out : inner; + Eigen::Index bj = IsRM ? inner : out; + auto dst_seg = dst.template middleCols(bj * BC); + auto lhs_seg = lhs.template middleCols
(bi * BR); + if (EIGEN_PREDICT_TRUE(a1)) + dst_seg.noalias() += lhs_seg * rhs.blockRef(id); + else if (am1) + dst_seg.noalias() -= lhs_seg * rhs.blockRef(id); + else { + typedef Matrix TmpType; + TmpType tmp(lhs.rows(), BC); + tmp.noalias() = lhs_seg * rhs.blockRef(id); + dst_seg.noalias() += alpha * tmp; + } + } + } + } +}; + +} // namespace internal + +} // end namespace Eigen + +#endif // EIGEN_BLOCKSPARSEMATRIX_H diff --git a/benchmarks/Sparse/CMakeLists.txt b/benchmarks/Sparse/CMakeLists.txt index eba60bfc3..19b25b95c 100644 --- a/benchmarks/Sparse/CMakeLists.txt +++ b/benchmarks/Sparse/CMakeLists.txt @@ -8,3 +8,4 @@ eigen_add_benchmark(bench_sparseview_assign bench_sparseview_assign.cpp) eigen_add_benchmark(bench_sparse_solvers bench_sparse_solvers.cpp) eigen_add_benchmark(bench_sparseqr_lookahead bench_sparseqr_lookahead.cpp) eigen_add_benchmark(bench_threaded_spmv bench_threaded_spmv.cpp LIBRARIES Threads::Threads) +eigen_add_benchmark(bench_block_sparse bench_block_sparse.cpp) diff --git a/benchmarks/Sparse/bench_block_sparse.cpp b/benchmarks/Sparse/bench_block_sparse.cpp new file mode 100644 index 000000000..e089e70d7 --- /dev/null +++ b/benchmarks/Sparse/bench_block_sparse.cpp @@ -0,0 +1,672 @@ +// 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 +#include + +#include +#include +#include + +using namespace Eigen; +using cd = std::complex; +using cf = std::complex; + +// --------------------------------------------------------------------------- +// 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, complex. +// --------------------------------------------------------------------------- +template +static typename std::enable_if::value, Scalar>::type randVal( + std::mt19937& rng, std::normal_distribution& d) { + return Scalar(d(rng)); +} + +template +static typename std::enable_if::value, Scalar>::type randVal( + std::mt19937& rng, std::normal_distribution& d) { + using R = typename Scalar::value_type; + return Scalar{R(d(rng)), R(d(rng))}; +} + +// --------------------------------------------------------------------------- +// Build a general (full) block-sparse pair. +// --------------------------------------------------------------------------- +template +static void buildPair(int nB, int nnzPerCol, unsigned seed, BlockSparseMatrix& bsm, + SparseMatrix& sm) { + using BSM = BlockSparseMatrix; + using BT = typename BSM::BlockType; + using Triplet = typename BSM::TripletType; + + std::mt19937 rng(seed); + std::uniform_int_distribution rowDist(0, nB - 1); + std::normal_distribution vd; + + std::vector triplets; + triplets.reserve(nB * nnzPerCol); + for (int j = 0; j < nB; ++j) { + std::set 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(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 +static void buildUpperTriPair(int nB, int nnzPerCol, unsigned seed, BlockSparseMatrix& bsm, + SparseMatrix& sm, bool forSolve = false) { + using BSM = BlockSparseMatrix; + using BT = typename BSM::BlockType; + using Triplet = typename BSM::TripletType; + + std::mt19937 rng(seed); + std::uniform_int_distribution rowDist(0, nB - 1); + std::normal_distribution vd; + + std::vector 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(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 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(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 +static void buildActuallyTriPair(int nB, int nnzPerCol, unsigned seed, BlockSparseMatrix& bsm, + SparseMatrix& sm) { + using BSM = BlockSparseMatrix; + using BT = typename BSM::BlockType; + using Triplet = typename BSM::TripletType; + + std::mt19937 rng(seed); + std::uniform_int_distribution rowDist(0, nB - 1); + std::normal_distribution vd; + + std::vector 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(rng, vd); + triplets.emplace_back(j, j, diag); + + std::set 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(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 +static void buildHermDiagUpperTriPair(int nB, int nnzPerCol, unsigned seed, + BlockSparseMatrix& bsm, SparseMatrix& sm) { + using BSM = BlockSparseMatrix; + using BT = typename BSM::BlockType; + using Triplet = typename BSM::TripletType; + + std::mt19937 rng(seed); + std::uniform_int_distribution rowDist(0, nB - 1); + std::normal_distribution vd; + + std::vector 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(rng, vd); + BT diag = (raw + raw.adjoint()) / Scalar(2); + triplets.emplace_back(j, j, diag); + + std::set 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(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 +static void BM_Sm_Sm_Add(benchmark::State& state) { + int nB, nnz; + parseArgs(state, nB, nnz); + SparseMatrix smA, smB, smC; + BlockSparseMatrix tmp; + buildPair(nB, nnz, 1, tmp, smA); + buildPair(nB, nnz, 2, tmp, smB); + for (auto _ : state) { + smC = smA + smB; + benchmark::DoNotOptimize(smC.valuePtr()); + } + state.counters["n"] = smA.rows(); +} + +template +static void BM_BSM_BSM_Add(benchmark::State& state) { + int nB, nnz; + parseArgs(state, nB, nnz); + using BSM = BlockSparseMatrix; + BSM bsmA, bsmB, bsmC; + SparseMatrix smTmp; + buildPair(nB, nnz, 1, bsmA, smTmp); + buildPair(nB, nnz, 2, bsmB, smTmp); + for (auto _ : state) { + bsmC = bsmA + bsmB; + benchmark::DoNotOptimize(bsmC.valuePtr()); + } + state.counters["n"] = bsmA.rows(); +} + +// --------------------------------------------------------------------------- +// Sparse×Dense: GEMV +// --------------------------------------------------------------------------- +template +static void BM_Sm_SpMV(benchmark::State& state) { + int nB, nnz; + parseArgs(state, nB, nnz); + SparseMatrix sm; + BlockSparseMatrix tmp; + buildPair(nB, nnz, 1, tmp, sm); + Matrix x = Matrix::Random(sm.cols()); + Matrix y(sm.rows()); + for (auto _ : state) { + y.noalias() = sm * x; + benchmark::DoNotOptimize(y.data()); + } + state.counters["n"] = sm.rows(); +} + +template +static void BM_BSM_SpMV(benchmark::State& state) { + int nB, nnz; + parseArgs(state, nB, nnz); + using BSM = BlockSparseMatrix; + BSM bsm; + SparseMatrix smTmp; + buildPair(nB, nnz, 1, bsm, smTmp); + Matrix x = Matrix::Random(bsm.cols()); + Matrix 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 +static void BM_Sm_SpMV_Acc(benchmark::State& state) { + int nB, nnz; + parseArgs(state, nB, nnz); + SparseMatrix sm; + BlockSparseMatrix tmp; + buildPair(nB, nnz, 1, tmp, sm); + Matrix x = Matrix::Random(sm.cols()); + Matrix y = Matrix::Random(sm.rows()); + for (auto _ : state) { + y.noalias() += sm * x; + benchmark::DoNotOptimize(y.data()); + } + state.counters["n"] = sm.rows(); +} + +template +static void BM_BSM_SpMV_Acc(benchmark::State& state) { + int nB, nnz; + parseArgs(state, nB, nnz); + using BSM = BlockSparseMatrix; + BSM bsm; + SparseMatrix smTmp; + buildPair(nB, nnz, 1, bsm, smTmp); + Matrix x = Matrix::Random(bsm.cols()); + Matrix y = Matrix::Random(bsm.rows()); + for (auto _ : state) { + y.noalias() += bsm * x; + benchmark::DoNotOptimize(y.data()); + } + state.counters["n"] = bsm.rows(); +} + +// --------------------------------------------------------------------------- +// Sparse×Dense: Triangular MV +// --------------------------------------------------------------------------- +template +static void BM_Sm_TriMV(benchmark::State& state) { + int nB, nnz; + parseArgs(state, nB, nnz); + SparseMatrix sm; + BlockSparseMatrix tmp; + buildUpperTriPair(nB, nnz, 1, tmp, sm); + Matrix x = Matrix::Random(sm.cols()); + Matrix y(sm.rows()); + for (auto _ : state) { + y.noalias() = sm.template triangularView() * x; + benchmark::DoNotOptimize(y.data()); + } + state.counters["n"] = sm.rows(); +} + +// BSM triangular view, DiagIsTriangular=false: diagonal blocks treated as triangular via triangularView<>. +template +static void BM_BSM_TriMV(benchmark::State& state) { + int nB, nnz; + parseArgs(state, nB, nnz); + using BSM = BlockSparseMatrix; + BSM bsm; + SparseMatrix smTmp; + buildUpperTriPair(nB, nnz, 1, bsm, smTmp); + Matrix x = Matrix::Random(bsm.cols()); + Matrix y(bsm.rows()); + for (auto _ : state) { + y.noalias() = bsm.template triangularView() * x; + benchmark::DoNotOptimize(y.data()); + } + state.counters["n"] = bsm.rows(); +} + +// BSM triangular view, DiagIsTriangular=true: diagonal blocks are actually triangular, uses full GEMV. +template +static void BM_BSM_TriMV_DiagT(benchmark::State& state) { + int nB, nnz; + parseArgs(state, nB, nnz); + using BSM = BlockSparseMatrix; + BSM bsm; + SparseMatrix smTmp; + buildActuallyTriPair(nB, nnz, 1, bsm, smTmp); + Matrix x = Matrix::Random(bsm.cols()); + Matrix y(bsm.rows()); + for (auto _ : state) { + y.noalias() = bsm.template triangularView() * x; + benchmark::DoNotOptimize(y.data()); + } + state.counters["n"] = bsm.rows(); +} + +// Accumulating triangular MV: y += A * x. +template +static void BM_Sm_TriMV_Acc(benchmark::State& state) { + int nB, nnz; + parseArgs(state, nB, nnz); + SparseMatrix sm; + BlockSparseMatrix tmp; + buildUpperTriPair(nB, nnz, 1, tmp, sm); + Matrix x = Matrix::Random(sm.cols()); + Matrix y = Matrix::Random(sm.rows()); + for (auto _ : state) { + y.noalias() += sm.template triangularView() * x; + benchmark::DoNotOptimize(y.data()); + } + state.counters["n"] = sm.rows(); +} + +template +static void BM_BSM_TriMV_Acc(benchmark::State& state) { + int nB, nnz; + parseArgs(state, nB, nnz); + using BSM = BlockSparseMatrix; + BSM bsm; + SparseMatrix smTmp; + buildUpperTriPair(nB, nnz, 1, bsm, smTmp); + Matrix x = Matrix::Random(bsm.cols()); + Matrix y = Matrix::Random(bsm.rows()); + for (auto _ : state) { + y.noalias() += bsm.template triangularView() * x; + benchmark::DoNotOptimize(y.data()); + } + state.counters["n"] = bsm.rows(); +} + +template +static void BM_BSM_TriMV_DiagT_Acc(benchmark::State& state) { + int nB, nnz; + parseArgs(state, nB, nnz); + using BSM = BlockSparseMatrix; + BSM bsm; + SparseMatrix smTmp; + buildActuallyTriPair(nB, nnz, 1, bsm, smTmp); + Matrix x = Matrix::Random(bsm.cols()); + Matrix y = Matrix::Random(bsm.rows()); + for (auto _ : state) { + y.noalias() += bsm.template triangularView() * x; + benchmark::DoNotOptimize(y.data()); + } + state.counters["n"] = bsm.rows(); +} + +// --------------------------------------------------------------------------- +// Sparse×Dense: Selfadjoint MV +// --------------------------------------------------------------------------- +template +static void BM_Sm_SymmMV(benchmark::State& state) { + int nB, nnz; + parseArgs(state, nB, nnz); + SparseMatrix sm; + BlockSparseMatrix tmp; + buildUpperTriPair(nB, nnz, 1, tmp, sm); + Matrix x = Matrix::Random(sm.cols()); + Matrix y(sm.rows()); + for (auto _ : state) { + y.noalias() = sm.template selfadjointView() * x; + benchmark::DoNotOptimize(y.data()); + } + state.counters["n"] = sm.rows(); +} + +// BSM selfadjoint view, general (non-Hermitian) diagonal blocks, DiagIsSA=false. +template +static void BM_BSM_SymmMV(benchmark::State& state) { + int nB, nnz; + parseArgs(state, nB, nnz); + using BSM = BlockSparseMatrix; + BSM bsm; + SparseMatrix smTmp; + buildUpperTriPair(nB, nnz, 1, bsm, smTmp); + Matrix x = Matrix::Random(bsm.cols()); + Matrix y(bsm.rows()); + for (auto _ : state) { + y.noalias() = bsm.template selfadjointView() * x; + benchmark::DoNotOptimize(y.data()); + } + state.counters["n"] = bsm.rows(); +} + +// BSM selfadjoint view, Hermitian diagonal blocks, DiagIsSA=false: fills diagonal via selfadjointView<>. +template +static void BM_BSM_SymmMV_DiagNSA(benchmark::State& state) { + int nB, nnz; + parseArgs(state, nB, nnz); + using BSM = BlockSparseMatrix; + BSM bsm; + SparseMatrix smTmp; + buildHermDiagUpperTriPair(nB, nnz, 1, bsm, smTmp); + Matrix x = Matrix::Random(bsm.cols()); + Matrix y(bsm.rows()); + for (auto _ : state) { + y.noalias() = bsm.template selfadjointView() * x; + benchmark::DoNotOptimize(y.data()); + } + state.counters["n"] = bsm.rows(); +} + +// BSM selfadjoint view, Hermitian diagonal blocks, DiagIsSA=true: full GEMV for diagonal blocks. +template +static void BM_BSM_SymmMV_DiagSA(benchmark::State& state) { + int nB, nnz; + parseArgs(state, nB, nnz); + using BSM = BlockSparseMatrix; + BSM bsm; + SparseMatrix smTmp; + buildHermDiagUpperTriPair(nB, nnz, 1, bsm, smTmp); + Matrix x = Matrix::Random(bsm.cols()); + Matrix y(bsm.rows()); + for (auto _ : state) { + y.noalias() = bsm.template selfadjointView() * x; + benchmark::DoNotOptimize(y.data()); + } + state.counters["n"] = bsm.rows(); +} + +// Accumulating selfadjoint MV: y += A * x. +template +static void BM_Sm_SymmMV_Acc(benchmark::State& state) { + int nB, nnz; + parseArgs(state, nB, nnz); + SparseMatrix sm; + BlockSparseMatrix tmp; + buildUpperTriPair(nB, nnz, 1, tmp, sm); + Matrix x = Matrix::Random(sm.cols()); + Matrix y = Matrix::Random(sm.rows()); + for (auto _ : state) { + y.noalias() += sm.template selfadjointView() * x; + benchmark::DoNotOptimize(y.data()); + } + state.counters["n"] = sm.rows(); +} + +template +static void BM_BSM_SymmMV_Acc(benchmark::State& state) { + int nB, nnz; + parseArgs(state, nB, nnz); + using BSM = BlockSparseMatrix; + BSM bsm; + SparseMatrix smTmp; + buildUpperTriPair(nB, nnz, 1, bsm, smTmp); + Matrix x = Matrix::Random(bsm.cols()); + Matrix y = Matrix::Random(bsm.rows()); + for (auto _ : state) { + y.noalias() += bsm.template selfadjointView() * x; + benchmark::DoNotOptimize(y.data()); + } + state.counters["n"] = bsm.rows(); +} + +template +static void BM_BSM_SymmMV_DiagNSA_Acc(benchmark::State& state) { + int nB, nnz; + parseArgs(state, nB, nnz); + using BSM = BlockSparseMatrix; + BSM bsm; + SparseMatrix smTmp; + buildHermDiagUpperTriPair(nB, nnz, 1, bsm, smTmp); + Matrix x = Matrix::Random(bsm.cols()); + Matrix y = Matrix::Random(bsm.rows()); + for (auto _ : state) { + y.noalias() += bsm.template selfadjointView() * x; + benchmark::DoNotOptimize(y.data()); + } + state.counters["n"] = bsm.rows(); +} + +template +static void BM_BSM_SymmMV_DiagSA_Acc(benchmark::State& state) { + int nB, nnz; + parseArgs(state, nB, nnz); + using BSM = BlockSparseMatrix; + BSM bsm; + SparseMatrix smTmp; + buildHermDiagUpperTriPair(nB, nnz, 1, bsm, smTmp); + Matrix x = Matrix::Random(bsm.cols()); + Matrix y = Matrix::Random(bsm.rows()); + for (auto _ : state) { + y.noalias() += bsm.template selfadjointView() * x; + benchmark::DoNotOptimize(y.data()); + } + state.counters["n"] = bsm.rows(); +} + +// --------------------------------------------------------------------------- +// Sparse×Dense: Triangular solve +// --------------------------------------------------------------------------- +template +static void BM_Sm_TriSolve(benchmark::State& state) { + int nB, nnz; + parseArgs(state, nB, nnz); + SparseMatrix sm; + BlockSparseMatrix tmp; + buildUpperTriPair(nB, nnz, 1, tmp, sm, true); + Matrix rhs = Matrix::Random(sm.cols()); + Matrix x(sm.cols()); + for (auto _ : state) { + x = rhs; + sm.template triangularView().solveInPlace(x); + benchmark::DoNotOptimize(x.data()); + } + state.counters["n"] = sm.rows(); +} + +template +static void BM_BSM_TriSolve(benchmark::State& state) { + int nB, nnz; + parseArgs(state, nB, nnz); + using BSM = BlockSparseMatrix; + BSM bsm; + SparseMatrix smTmp; + buildUpperTriPair(nB, nnz, 1, bsm, smTmp, true); + Matrix rhs = Matrix::Random(bsm.cols()); + Matrix x(bsm.cols()); + for (auto _ : state) { + x = rhs; + bsm.template triangularView().solveInPlace(x); + benchmark::DoNotOptimize(x.data()); + } + state.counters["n"] = bsm.rows(); +} + +// --------------------------------------------------------------------------- +// Sparse×Sparse: Matrix multiply +// --------------------------------------------------------------------------- +template +static void BM_Sm_Sm_Mul(benchmark::State& state) { + int nB, nnz; + parseArgs(state, nB, nnz); + SparseMatrix smA, smB, smC; + BlockSparseMatrix tmp; + buildPair(nB, nnz, 1, tmp, smA); + buildPair(nB, nnz, 2, tmp, smB); + for (auto _ : state) { + smC = smA * smB; + benchmark::DoNotOptimize(smC.valuePtr()); + } + state.counters["n"] = smA.rows(); +} + +template +static void BM_BSM_BSM_Mul(benchmark::State& state) { + int nB, nnz; + parseArgs(state, nB, nnz); + using BSM = BlockSparseMatrix; + BSM bsmA, bsmB, bsmC; + SparseMatrix smTmp; + buildPair(nB, nnz, 1, bsmA, smTmp); + buildPair(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)->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) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 80ed9eb35..081c915e2 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -310,6 +310,7 @@ ei_add_test(stdlist) ei_add_test(stdlist_overload) ei_add_test(stddeque) ei_add_test(stddeque_overload) +ei_add_test(block_sparse_matrix) ei_add_test(sparse_basic) ei_add_test(sparse_block) ei_add_test(sparse_vector) diff --git a/test/block_sparse_matrix.cpp b/test/block_sparse_matrix.cpp new file mode 100644 index 000000000..f20cae826 --- /dev/null +++ b/test/block_sparse_matrix.cpp @@ -0,0 +1,730 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. +// SPDX-FileCopyrightText: The Eigen Authors +// SPDX-License-Identifier: MPL-2.0 + +#include "sparse.h" + +using namespace Eigen; + +// --------------------------------------------------------------------------- +// Helper: build a BlockSparseMatrix from a dense matrix +// --------------------------------------------------------------------------- + +// Treats each BlockRows x BlockCols tile of `dense` as one block whenever the +// tile is non-zero. +template +BlockSparseMatrix denseToBlock( + const Matrix& dense) { + using BSM = BlockSparseMatrix; + using Triplet = typename BSM::TripletType; + + Index bRows = dense.rows() / BlockRows; + Index bCols = dense.cols() / BlockCols; + + std::vector triplets; + for (Index bi = 0; bi < bRows; ++bi) { + for (Index bj = 0; bj < bCols; ++bj) { + Matrix tile = dense.block(bi * BlockRows, bj * BlockCols, BlockRows, BlockCols); + if (tile.squaredNorm() > 0) triplets.emplace_back(StorageIndex(bi), StorageIndex(bj), tile); + } + } + + BSM bsm(bRows, bCols); + bsm.setFromTriplets(triplets.begin(), triplets.end()); + return bsm; +} + +// --------------------------------------------------------------------------- +// Core test driver templated on block size, storage order, and scalar type +// --------------------------------------------------------------------------- + +template +void test_block_sparse(int bRows, int bCols) { + using StorageIndex = int; + using BSM = BlockSparseMatrix; + using SpMat = SparseMatrix; + using DenseMat = Matrix; + + int rows = bRows * BlockRows; + int cols = bCols * BlockCols; + + // Build two random dense matrices with a coarse block structure. + DenseMat dA = DenseMat::Zero(rows, cols); + DenseMat dB = DenseMat::Zero(rows, cols); + for (int bi = 0; bi < bRows; ++bi) { + for (int bj = 0; bj < bCols; ++bj) { + if (internal::random(0.0, 1.0) < 0.4) { + dA.block(bi * BlockRows, bj * BlockCols, BlockRows, BlockCols) = DenseMat::Random(BlockRows, BlockCols); + } + if (internal::random(0.0, 1.0) < 0.4) { + dB.block(bi * BlockRows, bj * BlockCols, BlockRows, BlockCols) = DenseMat::Random(BlockRows, BlockCols); + } + } + } + + BSM A = denseToBlock(dA); + BSM B = denseToBlock(dB); + + // ---- toSparse / fromSparse round-trip ----------------------------------- + { + SpMat spA = A.toSparse(); + DenseMat dense_spA(spA); + VERIFY_IS_APPROX(dense_spA, dA); + + BSM A2 = BSM::fromSparse(spA); + VERIFY_IS_APPROX(A2.toSparse(), spA); + + // Implicit conversion operator + SpMat spA3 = A; + VERIFY_IS_APPROX(DenseMat(spA3), dA); + } + + // ---- Addition ----------------------------------------------------------- + { + BSM C = A + B; + DenseMat dC(C.toSparse()); + VERIFY_IS_APPROX(dC, dA + dB); + + BSM D = A; + D += B; + VERIFY_IS_APPROX(DenseMat(D.toSparse()), dA + dB); + } + + // ---- Subtraction -------------------------------------------------------- + { + BSM C = A - B; + VERIFY_IS_APPROX(DenseMat(C.toSparse()), dA - dB); + + BSM D = A; + D -= B; + VERIFY_IS_APPROX(DenseMat(D.toSparse()), dA - dB); + } + + // ---- Unary minus -------------------------------------------------------- + { + BSM C = -A; + VERIFY_IS_APPROX(DenseMat(C.toSparse()), -dA); + } + + // ---- cwiseProduct (conjunction / intersection) -------------------------- + { + BSM C = A.cwiseProduct(B); + // Result must only have blocks where both A and B had blocks. + DenseMat dC = dA.cwiseProduct(dB); + VERIFY_IS_APPROX(DenseMat(C.toSparse()), dC); + } + + // ---- unaryExpr ---------------------------------------------------------- + { + BSM C = A.unaryExpr([](const Scalar& x) { return x * x; }); + VERIFY_IS_APPROX(DenseMat(C.toSparse()), dA.array().square().matrix()); + } + + // ---- disjunctionExpr (union sparsity) ----------------------------------- + { + // sum via disjunctionExpr should match operator+ + struct AddExpr { + Scalar operator()(const Scalar& a, const Scalar& b) const { return a + b; } + Scalar lhs(const Scalar& a) const { return a; } + Scalar rhs(const Scalar& b) const { return b; } + }; + BSM C = A.disjunctionExpr(B, AddExpr{}); + VERIFY_IS_APPROX(DenseMat(C.toSparse()), dA + dB); + } + + // ---- conjunctionExpr (intersection sparsity) ---------------------------- + { + BSM C = A.conjunctionExpr(B, [](const Scalar& a, const Scalar& b) { return a * b; }); + VERIFY_IS_APPROX(DenseMat(C.toSparse()), dA.cwiseProduct(dB)); + } + + // ---- Scalar multiplication ---------------------------------------------- + { + Scalar s = Scalar(3.14); + BSM C = A * s; + VERIFY_IS_APPROX(DenseMat(C.toSparse()), dA * s); + + BSM D = s * A; + VERIFY_IS_APPROX(DenseMat(D.toSparse()), s * dA); + + BSM E = A; + E *= s; + VERIFY_IS_APPROX(DenseMat(E.toSparse()), dA * s); + } + + // ---- Element access (coeff) --------------------------------------------- + { + for (int i = 0; i < rows; ++i) { + for (int j = 0; j < cols; ++j) { + VERIFY_IS_APPROX(A.coeff(i, j), dA(i, j)); + } + } + } + + // ---- setIdentity ---------------------------------------------------------- + EIGEN_IF_CONSTEXPR (BlockRows == BlockCols) { + BSM Id(bRows, bCols); + Id.setIdentity(); + VERIFY_IS_APPROX(DenseMat(Id.toSparse()), DenseMat::Identity(rows, cols)); + if (bRows == bCols) { + VERIFY_IS_APPROX(DenseMat((Id * A).toSparse()), dA); + VERIFY_IS_APPROX(DenseMat((A * Id).toSparse()), dA); + } + } + + // ---- setFromTriplets with duplicate blocks (accumulation) --------------- + { + using Trip = typename BSM::TripletType; + using BlockMat = Matrix; + BlockMat half = BlockMat::Ones() * Scalar(0.5); + + std::vector trips; + trips.emplace_back(StorageIndex(0), StorageIndex(0), half); + trips.emplace_back(StorageIndex(0), StorageIndex(0), half); // duplicate -> sum + + BSM M(bRows, bCols); + M.setFromTriplets(trips.begin(), trips.end()); + + VERIFY(M.nonZeroBlocks() == 1); + VERIFY_IS_APPROX(DenseMat(M.blockRef(0)), DenseMat(BlockMat::Ones())); + } +} + +// --------------------------------------------------------------------------- +// Square-block product test +// --------------------------------------------------------------------------- + +template +void test_block_sparse_product(int bM, int bK, int bN) { + using StorageIndex = int; + using BSMA = BlockSparseMatrix; + using DenseMat = Matrix; + + int rowsA = bM * B, colsA = bK * B; + int colsB = bN * B; + + DenseMat dA = DenseMat::Zero(rowsA, colsA); + DenseMat dB = DenseMat::Zero(colsA, colsB); + + for (int bi = 0; bi < bM; ++bi) + for (int bk = 0; bk < bK; ++bk) + if (internal::random(0.0, 1.0) < 0.4) dA.block(bi * B, bk * B, B, B) = DenseMat::Random(B, B); + + for (int bk = 0; bk < bK; ++bk) + for (int bj = 0; bj < bN; ++bj) + if (internal::random(0.0, 1.0) < 0.4) dB.block(bk * B, bj * B, B, B) = DenseMat::Random(B, B); + + BSMA A = denseToBlock(dA); + BSMA Bmat = denseToBlock(dB); + + BSMA C = A * Bmat; + DenseMat dC(C.toSparse()); + VERIFY_IS_APPROX(dC, dA * dB); +} + +// --------------------------------------------------------------------------- +// Block-sparse * dense and dense * block-sparse products +// --------------------------------------------------------------------------- + +template +void test_block_sparse_dense_product(int bRows, int bCols) { + using StorageIndex = int; + using BSM = BlockSparseMatrix; + using DenseMat = Matrix; + using DenseVec = Matrix; + using RowVec = Matrix; + + int rows = bRows * BlockRows; + int cols = bCols * BlockCols; + + DenseMat dA = DenseMat::Zero(rows, cols); + for (int bi = 0; bi < bRows; ++bi) + for (int bj = 0; bj < bCols; ++bj) + if (internal::random(0.0, 1.0) < 0.5) + dA.block(bi * BlockRows, bj * BlockCols, BlockRows, BlockCols) = DenseMat::Random(BlockRows, BlockCols); + + BSM A = denseToBlock(dA); + + // BSM * dense matrix + { + DenseMat rhs = DenseMat::Random(cols, 5); + VERIFY_IS_APPROX(A * rhs, dA * rhs); + } + + // BSM * column vector + { + DenseVec v = DenseVec::Random(cols); + VERIFY_IS_APPROX(A * v, dA * v); + } + + // dense matrix * BSM + { + DenseMat lhs = DenseMat::Random(5, rows); + VERIFY_IS_APPROX(lhs * A, lhs * dA); + } + + // row vector * BSM + { + RowVec v = RowVec::Random(rows); + VERIFY_IS_APPROX(v * A, v * dA); + } +} + +// --------------------------------------------------------------------------- +// Non-square block product: A(BR x BC) * B(BC x BC2) -> C(BR x BC2) +// --------------------------------------------------------------------------- + +void test_nonsquare_block_product() { + using Scalar = double; + using StorageIndex = int; + constexpr int BR = 2, BC = 3, BC2 = 4; + using BSMA = BlockSparseMatrix; + using BSMB = BlockSparseMatrix; + using DenseMat = Matrix; + + int bM = 4, bK = 3, bN = 5; + DenseMat dA = DenseMat::Zero(bM * BR, bK * BC); + DenseMat dB = DenseMat::Zero(bK * BC, bN * BC2); + + for (int bi = 0; bi < bM; ++bi) + for (int bk = 0; bk < bK; ++bk) + if (internal::random() > 0.0) dA.block(bi * BR, bk * BC, BR, BC) = DenseMat::Random(BR, BC); + + for (int bk = 0; bk < bK; ++bk) + for (int bj = 0; bj < bN; ++bj) + if (internal::random() > 0.0) dB.block(bk * BC, bj * BC2, BC, BC2) = DenseMat::Random(BC, BC2); + + BSMA A = denseToBlock(dA); + BSMB Bmat = denseToBlock(dB); + + BlockSparseMatrix C = A * Bmat; + DenseMat dC(C.toSparse()); + VERIFY_IS_APPROX(dC, dA * dB); +} + +// --------------------------------------------------------------------------- +// Transpose and adjoint +// --------------------------------------------------------------------------- + +template +void test_block_sparse_transpose(int bRows, int bCols) { + using StorageIndex = int; + using BSM = BlockSparseMatrix; + using BSMT = BlockSparseMatrix; + using DenseMat = Matrix; + + int rows = bRows * BlockRows, cols = bCols * BlockCols; + DenseMat dA = DenseMat::Zero(rows, cols); + for (int bi = 0; bi < bRows; ++bi) + for (int bj = 0; bj < bCols; ++bj) + if (internal::random(0.0, 1.0) < 0.5) + dA.block(bi * BlockRows, bj * BlockCols, BlockRows, BlockCols) = DenseMat::Random(BlockRows, BlockCols); + + BSM A = denseToBlock(dA); + + // transpose + BSMT At = A.transpose(); + VERIFY_IS_APPROX(DenseMat(At.toSparse()), dA.transpose()); + + // adjoint (conjugate transpose for complex, same as transpose for real) + BSMT Ah = A.adjoint(); + VERIFY_IS_APPROX(DenseMat(Ah.toSparse()), dA.adjoint()); + + // (A^T)^T == A + BSM AtT = At.transpose(); + VERIFY_IS_APPROX(DenseMat(AtT.toSparse()), dA); +} + +// --------------------------------------------------------------------------- +// Triangular view: eval, +/-, dense products, DiagIsTriangular path +// --------------------------------------------------------------------------- + +template +void test_block_sparse_triangular(int bN) { + using StorageIndex = int; + using BSM = BlockSparseMatrix; + using DenseMat = Matrix; + + int N = bN * B; + DenseMat dA = DenseMat::Zero(N, N); + for (int bi = 0; bi < bN; ++bi) + for (int bj = 0; bj < bN; ++bj) + if (internal::random(0.0, 1.0) < 0.5) dA.block(bi * B, bj * B, B, B) = DenseMat::Random(B, B); + + BSM A = denseToBlock(dA); + + // Upper eval + { + BSM Au = A.template triangularView().eval(); + DenseMat dAu = dA; + for (int bi = 0; bi < bN; ++bi) { + for (int bj = 0; bj < bi; ++bj) dAu.block(bi * B, bj * B, B, B).setZero(); + dAu.block(bi * B, bi * B, B, B).template triangularView().setZero(); + } + VERIFY_IS_APPROX(DenseMat(Au.toSparse()), dAu); + } + + // Lower eval + { + BSM Al = A.template triangularView().eval(); + DenseMat dAl = dA; + for (int bi = 0; bi < bN; ++bi) { + for (int bj = bi + 1; bj < bN; ++bj) dAl.block(bi * B, bj * B, B, B).setZero(); + dAl.block(bi * B, bi * B, B, B).template triangularView().setZero(); + } + VERIFY_IS_APPROX(DenseMat(Al.toSparse()), dAl); + } + + // Tri * dense and dense * Tri + { + DenseMat dAu = dA; + for (int bi = 0; bi < bN; ++bi) { + for (int bj = 0; bj < bi; ++bj) dAu.block(bi * B, bj * B, B, B).setZero(); + dAu.block(bi * B, bi * B, B, B).template triangularView().setZero(); + } + + DenseMat rhs = DenseMat::Random(N, 4); + DenseMat lhs = DenseMat::Random(3, N); + + VERIFY_IS_APPROX(A.template triangularView() * rhs, dAu * rhs); + VERIFY_IS_APPROX(lhs * A.template triangularView(), lhs * dAu); + } + + // Tri + Tri + { + DenseMat dB = DenseMat::Zero(N, N); + for (int bi = 0; bi < bN; ++bi) + for (int bj = 0; bj < bN; ++bj) + if (internal::random(0.0, 1.0) < 0.5) dB.block(bi * B, bj * B, B, B) = DenseMat::Random(B, B); + BSM Bmat = denseToBlock(dB); + + DenseMat dAu = dA, dBu = dB; + for (int bi = 0; bi < bN; ++bi) { + for (int bj = 0; bj < bi; ++bj) { + dAu.block(bi * B, bj * B, B, B).setZero(); + dBu.block(bi * B, bj * B, B, B).setZero(); + } + dAu.block(bi * B, bi * B, B, B).template triangularView().setZero(); + dBu.block(bi * B, bi * B, B, B).template triangularView().setZero(); + } + BSM C = A.template triangularView() + Bmat.template triangularView(); + VERIFY_IS_APPROX(DenseMat(C.toSparse()), dAu + dBu); + } + + // DiagIsTriangular=true product: build a BSM whose diagonal blocks are already + // upper-triangular in storage, and verify the product matches DiagIsTriangular=false. + { + DenseMat dAu = DenseMat::Zero(N, N); + for (int bi = 0; bi < bN; ++bi) { + DenseMat blk = DenseMat::Random(B, B); + blk.template triangularView().setZero(); + dAu.block(bi * B, bi * B, B, B) = blk; + for (int bj = bi + 1; bj < bN; ++bj) + if (internal::random(0.0, 1.0) < 0.5) dAu.block(bi * B, bj * B, B, B) = DenseMat::Random(B, B); + } + BSM Au = denseToBlock(dAu); + DenseMat rhs = DenseMat::Random(N, 3); + DenseMat r1 = Au.template triangularView() * rhs; + DenseMat r2 = Au.template triangularView() * rhs; + VERIFY_IS_APPROX(r1, r2); + } +} + +// --------------------------------------------------------------------------- +// Triangular solve: forward/backward x direct/transposed/adjoint x both layouts +// --------------------------------------------------------------------------- + +template +void test_block_sparse_triangular_solve(int bN) { + using StorageIndex = int; + using BSM = BlockSparseMatrix; + using DenseMat = Matrix; + using RealScalar = typename NumTraits::Real; + + int N = bN * B; + + // Build a dense lower-block-triangular matrix. + // Diagonal blocks are lower triangular with non-zero diagonal. + // Off-diagonal (lower) blocks are random dense. + auto makeDenseLower = [&]() { + DenseMat dL = DenseMat::Zero(N, N); + for (int bi = 0; bi < bN; ++bi) { + DenseMat blk = DenseMat::Random(B, B); + blk.template triangularView().setZero(); + for (int k = 0; k < B; ++k) blk(k, k) = Scalar(RealScalar(B + k + 1)); + dL.block(bi * B, bi * B, B, B) = blk; + for (int bj = 0; bj < bi; ++bj) + if (internal::random(0.0, 1.0) < 0.6) dL.block(bi * B, bj * B, B, B) = DenseMat::Random(B, B); + } + return dL; + }; + + auto makeDenseUpper = [&]() { + DenseMat dU = DenseMat::Zero(N, N); + for (int bi = 0; bi < bN; ++bi) { + DenseMat blk = DenseMat::Random(B, B); + blk.template triangularView().setZero(); + for (int k = 0; k < B; ++k) blk(k, k) = Scalar(RealScalar(B + k + 1)); + dU.block(bi * B, bi * B, B, B) = blk; + for (int bj = bi + 1; bj < bN; ++bj) + if (internal::random(0.0, 1.0) < 0.6) dU.block(bi * B, bj * B, B, B) = DenseMat::Random(B, B); + } + return dU; + }; + + // Lower triangular: direct, transposed, adjoint + { + DenseMat dL = makeDenseLower(); + BSM L = denseToBlock(dL); + + // L x = b + { + DenseMat b = DenseMat::Random(N, 3); + DenseMat x = b; + L.template triangularView().solveInPlace(x); + VERIFY_IS_APPROX(dL * x, b); + } + // L^T x = b + { + DenseMat b = DenseMat::Random(N, 3); + DenseMat x = b; + L.template triangularView().transpose().solveInPlace(x); + VERIFY_IS_APPROX(dL.transpose() * x, b); + } + // L^H x = b + { + DenseMat b = DenseMat::Random(N, 3); + DenseMat x = b; + L.template triangularView().adjoint().solveInPlace(x); + VERIFY_IS_APPROX(dL.adjoint() * x, b); + } + } + + // Upper triangular: direct, transposed, adjoint + { + DenseMat dU = makeDenseUpper(); + BSM U = denseToBlock(dU); + + // U x = b + { + DenseMat b = DenseMat::Random(N, 3); + DenseMat x = b; + U.template triangularView().solveInPlace(x); + VERIFY_IS_APPROX(dU * x, b); + } + // U^T x = b + { + DenseMat b = DenseMat::Random(N, 3); + DenseMat x = b; + U.template triangularView().transpose().solveInPlace(x); + VERIFY_IS_APPROX(dU.transpose() * x, b); + } + // U^H x = b + { + DenseMat b = DenseMat::Random(N, 3); + DenseMat x = b; + U.template triangularView().adjoint().solveInPlace(x); + VERIFY_IS_APPROX(dU.adjoint() * x, b); + } + } + + // DiagIsTriangular=true: diagonal blocks are already properly triangular in storage; + // result must match DiagIsTriangular=false (which zeroes the unused triangle first). + { + DenseMat dL = makeDenseLower(); + BSM L = denseToBlock(dL); + DenseMat b = DenseMat::Random(N, 2); + DenseMat x1 = b, x2 = b; + { + auto tri_false = L.template triangularView(); + tri_false.solveInPlace(x1); + } + { + auto tri_true = L.template triangularView(); + tri_true.solveInPlace(x2); + } + VERIFY_IS_APPROX(x1, x2); + } +} + +// --------------------------------------------------------------------------- +// Self-adjoint view: eval, +/-, dense products +// --------------------------------------------------------------------------- + +template +void test_block_sparse_selfadjoint(int bN) { + using StorageIndex = int; + using BSM = BlockSparseMatrix; + using DenseMat = Matrix; + + int N = bN * B; + + // Build a Hermitian dense matrix (stored upper triangle only). + DenseMat dFull = DenseMat::Zero(N, N); + for (int bi = 0; bi < bN; ++bi) { + // Diagonal block: Hermitian (symmetric for real, conjugate-symmetric for complex). + DenseMat blk = DenseMat::Random(B, B); + blk = (blk + blk.adjoint()).eval(); + dFull.block(bi * B, bi * B, B, B) = blk; + for (int bj = bi + 1; bj < bN; ++bj) + if (internal::random(0.0, 1.0) < 0.5) { + DenseMat offblk = DenseMat::Random(B, B); + dFull.block(bi * B, bj * B, B, B) = offblk; + dFull.block(bj * B, bi * B, B, B) = offblk.adjoint(); // conjugate transpose mirror + } + } + + // Build BSM from the upper triangle only. + DenseMat dUpper = DenseMat::Zero(N, N); + for (int bi = 0; bi < bN; ++bi) + for (int bj = bi; bj < bN; ++bj) dUpper.block(bi * B, bj * B, B, B) = dFull.block(bi * B, bj * B, B, B); + + BSM A = denseToBlock(dUpper); + + // eval() must reproduce the full Hermitian matrix + { + BSM Asym = A.template selfadjointView().eval(); + VERIFY_IS_APPROX(DenseMat(Asym.toSparse()), dFull); + } + + // selfadjointView * dense + { + DenseMat rhs = DenseMat::Random(N, 5); + VERIFY_IS_APPROX(A.template selfadjointView() * rhs, dFull * rhs); + } + + // dense * selfadjointView + { + DenseMat lhs = DenseMat::Random(4, N); + VERIFY_IS_APPROX(lhs * A.template selfadjointView(), lhs * dFull); + } + + // selfadjointView + selfadjointView + { + BSM Bmat = denseToBlock(dUpper * Scalar(2)); + BSM C = A.template selfadjointView() + Bmat.template selfadjointView(); + VERIFY_IS_APPROX(DenseMat(C.toSparse()), dFull * Scalar(3)); + } + + // DiagIsSelfAdjoint path: diagonal blocks ARE Hermitian, product should match + { + DenseMat rhs = DenseMat::Random(N, 3); + DenseMat result = A.template selfadjointView() * rhs; + VERIFY_IS_APPROX(result, dFull * rhs); + } +} + +// --------------------------------------------------------------------------- +// BlockTriplet type-trait checks +// --------------------------------------------------------------------------- + +void test_block_triplet_traits() { + // Flat scalar array means BlockTriplet should be trivially copyable and + // standard-layout for any trivially-copyable Scalar and StorageIndex. + EIGEN_STATIC_ASSERT((std::is_trivially_copyable>::value), + BLOCKTRIPLET_MUST_BE_TRIVIALLY_COPYABLE) + EIGEN_STATIC_ASSERT((std::is_trivially_copyable>::value), + BLOCKTRIPLET_MUST_BE_TRIVIALLY_COPYABLE) + EIGEN_STATIC_ASSERT((std::is_trivially_copyable>::value), + BLOCKTRIPLET_MUST_BE_TRIVIALLY_COPYABLE) + EIGEN_STATIC_ASSERT((std::is_standard_layout>::value), BLOCKTRIPLET_MUST_BE_STANDARD_LAYOUT) + EIGEN_STATIC_ASSERT((std::is_standard_layout>::value), + BLOCKTRIPLET_MUST_BE_STANDARD_LAYOUT) + + // No alignment padding: size must equal 2*sizeof(StorageIndex) + BlockSize*sizeof(Scalar). + EIGEN_STATIC_ASSERT((sizeof(BlockTriplet) == 2 * sizeof(int) + 4 * sizeof(float)), + BLOCKTRIPLET_MUST_HAVE_NO_ALIGNMENT_PADDING) + EIGEN_STATIC_ASSERT((sizeof(BlockTriplet) == 2 * sizeof(int) + 4 * sizeof(double)), + BLOCKTRIPLET_MUST_HAVE_NO_ALIGNMENT_PADDING) + EIGEN_STATIC_ASSERT((sizeof(BlockTriplet) == 2 * sizeof(int) + 16 * sizeof(float)), + BLOCKTRIPLET_MUST_HAVE_NO_ALIGNMENT_PADDING) + EIGEN_STATIC_ASSERT((sizeof(BlockTriplet) == 2 * sizeof(int) + 16 * sizeof(double)), + BLOCKTRIPLET_MUST_HAVE_NO_ALIGNMENT_PADDING) +} + +// --------------------------------------------------------------------------- +// Main entry point +// --------------------------------------------------------------------------- + +EIGEN_DECLARE_TEST(block_sparse_matrix) { + // ColMajor, real double, various block sizes and matrix sizes + CALL_SUBTEST_1((test_block_sparse<1, 1, ColMajor>(6, 8))); + CALL_SUBTEST_2((test_block_sparse<2, 2, ColMajor>(4, 6))); + CALL_SUBTEST_3((test_block_sparse<3, 3, ColMajor>(4, 5))); + CALL_SUBTEST_4((test_block_sparse<4, 4, ColMajor>(3, 3))); + CALL_SUBTEST_5((test_block_sparse<2, 3, ColMajor>(5, 4))); + + // RowMajor, real double + CALL_SUBTEST_6((test_block_sparse<2, 2, RowMajor>(4, 6))); + CALL_SUBTEST_7((test_block_sparse<3, 3, RowMajor>(4, 5))); + CALL_SUBTEST_8((test_block_sparse<2, 3, RowMajor>(5, 4))); + + // Products (ColMajor) + CALL_SUBTEST_9((test_block_sparse_product<2, ColMajor>(4, 5, 3))); + CALL_SUBTEST_9((test_block_sparse_product<3, ColMajor>(3, 4, 5))); + + // Products (RowMajor) + CALL_SUBTEST_10((test_block_sparse_product<2, RowMajor>(4, 5, 3))); + CALL_SUBTEST_10((test_block_sparse_product<3, RowMajor>(3, 4, 5))); + + // Non-square block product + CALL_SUBTEST_11(test_nonsquare_block_product()); + + // Block-sparse * dense and dense * block-sparse (ColMajor) + CALL_SUBTEST_12((test_block_sparse_dense_product<2, 2, ColMajor>(4, 5))); + CALL_SUBTEST_12((test_block_sparse_dense_product<3, 3, ColMajor>(3, 4))); + CALL_SUBTEST_12((test_block_sparse_dense_product<2, 3, ColMajor>(4, 3))); + + // Block-sparse * dense and dense * block-sparse (RowMajor) + CALL_SUBTEST_13((test_block_sparse_dense_product<2, 2, RowMajor>(4, 5))); + CALL_SUBTEST_13((test_block_sparse_dense_product<3, 3, RowMajor>(3, 4))); + CALL_SUBTEST_13((test_block_sparse_dense_product<2, 3, RowMajor>(4, 3))); + + // Transpose / adjoint (ColMajor and RowMajor, square and non-square blocks) + CALL_SUBTEST_14((test_block_sparse_transpose<2, 2, ColMajor>(4, 5))); + CALL_SUBTEST_14((test_block_sparse_transpose<2, 3, ColMajor>(5, 4))); + CALL_SUBTEST_14((test_block_sparse_transpose<2, 2, RowMajor>(4, 5))); + CALL_SUBTEST_14((test_block_sparse_transpose<3, 2, RowMajor>(4, 5))); + + // Triangular view (ColMajor and RowMajor) + CALL_SUBTEST_15((test_block_sparse_triangular<2, ColMajor>(5))); + CALL_SUBTEST_15((test_block_sparse_triangular<3, ColMajor>(4))); + CALL_SUBTEST_15((test_block_sparse_triangular<2, RowMajor>(5))); + + // Self-adjoint view (ColMajor and RowMajor) + CALL_SUBTEST_16((test_block_sparse_selfadjoint<2, ColMajor>(5))); + CALL_SUBTEST_16((test_block_sparse_selfadjoint<3, ColMajor>(4))); + CALL_SUBTEST_16((test_block_sparse_selfadjoint<2, RowMajor>(5))); + + // BlockTriplet type traits + CALL_SUBTEST_17(test_block_triplet_traits()); + + // Complex scalar coverage: conjugation paths in adjoint, selfadjoint, triangular products + CALL_SUBTEST_18((test_block_sparse<2, 2, ColMajor, std::complex>(4, 6))); + CALL_SUBTEST_18((test_block_sparse<3, 3, ColMajor, std::complex>(4, 5))); + CALL_SUBTEST_18((test_block_sparse<2, 2, RowMajor, std::complex>(4, 6))); + CALL_SUBTEST_18((test_block_sparse_dense_product<2, 2, ColMajor, std::complex>(4, 5))); + CALL_SUBTEST_18((test_block_sparse_dense_product<2, 2, RowMajor, std::complex>(4, 5))); + CALL_SUBTEST_18((test_block_sparse_transpose<2, 2, ColMajor, std::complex>(4, 5))); + CALL_SUBTEST_18((test_block_sparse_transpose<2, 3, ColMajor, std::complex>(5, 4))); + CALL_SUBTEST_18((test_block_sparse_selfadjoint<2, ColMajor, std::complex>(5))); + CALL_SUBTEST_18((test_block_sparse_selfadjoint<3, ColMajor, std::complex>(4))); + CALL_SUBTEST_18((test_block_sparse_selfadjoint<2, RowMajor, std::complex>(5))); + CALL_SUBTEST_18((test_block_sparse_triangular<2, ColMajor, std::complex>(5))); + CALL_SUBTEST_18((test_block_sparse_triangular<2, RowMajor, std::complex>(5))); + + // Triangular solve: forward/backward x direct/transposed/adjoint x ColMajor+RowMajor + CALL_SUBTEST_19((test_block_sparse_triangular_solve<2, ColMajor, double>(5))); + CALL_SUBTEST_19((test_block_sparse_triangular_solve<3, ColMajor, double>(4))); + CALL_SUBTEST_19((test_block_sparse_triangular_solve<2, RowMajor, double>(5))); + CALL_SUBTEST_19((test_block_sparse_triangular_solve<3, RowMajor, double>(4))); + CALL_SUBTEST_19((test_block_sparse_triangular_solve<2, ColMajor, std::complex>(5))); + CALL_SUBTEST_19((test_block_sparse_triangular_solve<3, ColMajor, std::complex>(4))); + CALL_SUBTEST_19((test_block_sparse_triangular_solve<2, RowMajor, std::complex>(5))); + CALL_SUBTEST_19((test_block_sparse_triangular_solve<3, RowMajor, std::complex>(4))); +} diff --git a/unsupported/Eigen/src/SparseExtra/BlockSparseMatrix.h b/unsupported/Eigen/src/SparseExtra/BlockSparseMatrix.h deleted file mode 100644 index d1ef21673..000000000 --- a/unsupported/Eigen/src/SparseExtra/BlockSparseMatrix.h +++ /dev/null @@ -1,932 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2013 Desire Nuentsa -// Copyright (C) 2013 Gael Guennebaud -// -// This Source Code Form is subject to the terms of the Mozilla -// Public License v. 2.0. If a copy of the MPL was not distributed -// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. -// SPDX-License-Identifier: MPL-2.0 - -#ifndef EIGEN_SPARSEBLOCKMATRIX_H -#define EIGEN_SPARSEBLOCKMATRIX_H - -// IWYU pragma: private -#include "./InternalHeaderCheck.h" - -namespace Eigen { -/** \ingroup SparseCore_Module - * - * \class BlockSparseMatrix - * - * \brief A versatile sparse matrix representation where each element is a block - * - * This class provides routines to manipulate block sparse matrices stored in a - * BSR-like representation. There are two main types : - * - * 1. All blocks have the same number of rows and columns, called block size - * in the following. In this case, if this block size is known at compile time, - * it can be given as a template parameter like - * \code - * BlockSparseMatrix bmat(b_rows, b_cols); - * \endcode - * Here, bmat is a b_rows x b_cols block sparse matrix - * where each coefficient is a 3x3 dense matrix. - * If the block size is fixed but will be given at runtime, - * \code - * BlockSparseMatrix bmat(b_rows, b_cols); - * bmat.setBlockSize(block_size); - * \endcode - * - * 2. The second case is for variable-block sparse matrices. - * Here each block has its own dimensions. The only restriction is that all the blocks - * in a row (resp. a column) should have the same number of rows (resp. of columns). - * It is thus required in this case to describe the layout of the matrix by calling - * setBlockLayout(rowBlocks, colBlocks). - * - * In any of the previous cases, the matrix can be filled by calling setFromTriplets(). - * A regular sparse matrix can be converted to a block sparse matrix and vice versa. - * It is obviously required to describe the block layout beforehand by calling either - * setBlockSize() for fixed-size blocks or setBlockLayout for variable-size blocks. - * - * \tparam Scalar_ The Scalar type - * \tparam _BlockAtCompileTime The block layout option. It takes the following values - * Dynamic : block size known at runtime - * a numeric number : fixed-size block known at compile time - */ -template -class BlockSparseMatrix; - -template -class BlockSparseMatrixView; - -namespace internal { -template -struct traits > { - typedef Scalar_ Scalar; - typedef Index_ Index; - typedef Sparse StorageKind; // FIXME: determine where StorageKind is used. - typedef MatrixXpr XprKind; - enum { - RowsAtCompileTime = Dynamic, - ColsAtCompileTime = Dynamic, - MaxRowsAtCompileTime = Dynamic, - MaxColsAtCompileTime = Dynamic, - BlockSize = _BlockAtCompileTime, - Flags = Options_ | NestByRefBit | LvalueBit, - CoeffReadCost = NumTraits::ReadCost, - SupportedAccessPatterns = InnerRandomAccessPattern - }; -}; -template -struct traits > { - typedef Ref< - Matrix > - Scalar; - typedef Ref< - Matrix > - RealScalar; -}; - -// Function object to sort a triplet list -template -struct TripletComp { - typedef typename Iterator::value_type Triplet; - bool operator()(const Triplet& a, const Triplet& b) { - EIGEN_IF_CONSTEXPR (IsColMajor) { - return ((a.col() == b.col() && a.row() < b.row()) || (a.col() < b.col())); - } else { - return ((a.row() == b.row() && a.col() < b.col()) || (a.row() < b.row())); - } - } -}; -} // end namespace internal - -/* Proxy to view the block sparse matrix as a regular sparse matrix */ -template -class BlockSparseMatrixView : public SparseMatrixBase { - public: - typedef Ref Scalar; - typedef Ref RealScalar; - typedef typename BlockSparseMatrixT::Index Index; - typedef BlockSparseMatrixT Nested; - enum { - Flags = BlockSparseMatrixT::Options, - Options = BlockSparseMatrixT::Options, - RowsAtCompileTime = BlockSparseMatrixT::RowsAtCompileTime, - ColsAtCompileTime = BlockSparseMatrixT::ColsAtCompileTime, - MaxColsAtCompileTime = BlockSparseMatrixT::MaxColsAtCompileTime, - MaxRowsAtCompileTime = BlockSparseMatrixT::MaxRowsAtCompileTime - }; - - public: - BlockSparseMatrixView(const BlockSparseMatrixT& spblockmat) : m_spblockmat(spblockmat) {} - - Index outerSize() const { return (Flags & RowMajorBit) == 1 ? this->rows() : this->cols(); } - Index cols() const { return m_spblockmat.blockCols(); } - Index rows() const { return m_spblockmat.blockRows(); } - Scalar coeff(Index row, Index col) { return m_spblockmat.coeff(row, col); } - Scalar coeffRef(Index row, Index col) { return m_spblockmat.coeffRef(row, col); } - // Wrapper to iterate over all blocks - class InnerIterator : public BlockSparseMatrixT::BlockInnerIterator { - public: - InnerIterator(const BlockSparseMatrixView& mat, Index outer) - : BlockSparseMatrixT::BlockInnerIterator(mat.m_spblockmat, outer) {} - }; - - protected: - const BlockSparseMatrixT& m_spblockmat; -}; - -// Proxy to view a regular vector as a block vector -template -class BlockVectorView { - public: - enum { - BlockSize = BlockSparseMatrixT::BlockSize, - ColsAtCompileTime = VectorType::ColsAtCompileTime, - RowsAtCompileTime = VectorType::RowsAtCompileTime, - Flags = VectorType::Flags - }; - typedef Ref > - Scalar; - typedef typename BlockSparseMatrixT::Index Index; - - public: - BlockVectorView(const BlockSparseMatrixT& spblockmat, const VectorType& vec) : m_spblockmat(spblockmat), m_vec(vec) {} - inline Index cols() const { return m_vec.cols(); } - inline Index size() const { return m_spblockmat.blockRows(); } - inline Scalar coeff(Index bi) const { - Index startRow = m_spblockmat.blockRowsIndex(bi); - Index rowSize = m_spblockmat.blockRowsIndex(bi + 1) - startRow; - return m_vec.middleRows(startRow, rowSize); - } - inline Scalar coeff(Index bi, Index j) const { - Index startRow = m_spblockmat.blockRowsIndex(bi); - Index rowSize = m_spblockmat.blockRowsIndex(bi + 1) - startRow; - return m_vec.block(startRow, j, rowSize, 1); - } - - protected: - const BlockSparseMatrixT& m_spblockmat; - const VectorType& m_vec; -}; - -template -class BlockVectorReturn; - -// Proxy to view a regular vector as a block vector -template -class BlockVectorReturn { - public: - enum { - ColsAtCompileTime = VectorType::ColsAtCompileTime, - RowsAtCompileTime = VectorType::RowsAtCompileTime, - Flags = VectorType::Flags - }; - typedef Ref > Scalar; - typedef typename BlockSparseMatrixT::Index Index; - - public: - BlockVectorReturn(const BlockSparseMatrixT& spblockmat, VectorType& vec) : m_spblockmat(spblockmat), m_vec(vec) {} - inline Index size() const { return m_spblockmat.blockRows(); } - inline Scalar coeffRef(Index bi) { - Index startRow = m_spblockmat.blockRowsIndex(bi); - Index rowSize = m_spblockmat.blockRowsIndex(bi + 1) - startRow; - return m_vec.middleRows(startRow, rowSize); - } - inline Scalar coeffRef(Index bi, Index j) { - Index startRow = m_spblockmat.blockRowsIndex(bi); - Index rowSize = m_spblockmat.blockRowsIndex(bi + 1) - startRow; - return m_vec.block(startRow, j, rowSize, 1); - } - - protected: - const BlockSparseMatrixT& m_spblockmat; - VectorType& m_vec; -}; - -// Block version of the sparse dense product -template -class BlockSparseTimeDenseProduct; - -namespace internal { - -template -struct traits > { - typedef Dense StorageKind; - typedef MatrixXpr XprKind; - typedef typename BlockSparseMatrixT::Scalar Scalar; - typedef typename BlockSparseMatrixT::Index Index; - enum { - RowsAtCompileTime = Dynamic, - ColsAtCompileTime = Dynamic, - MaxRowsAtCompileTime = Dynamic, - MaxColsAtCompileTime = Dynamic, - Flags = 0, - CoeffReadCost = internal::traits::CoeffReadCost - }; -}; -} // end namespace internal - -template -class BlockSparseTimeDenseProduct : public ProductBase, Lhs, Rhs> { - public: - EIGEN_PRODUCT_PUBLIC_INTERFACE(BlockSparseTimeDenseProduct) - - BlockSparseTimeDenseProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs, rhs) {} - - template - void scaleAndAddTo(Dest& dest, const typename Rhs::Scalar& alpha) const { - BlockVectorReturn tmpDest(m_lhs, dest); - internal::sparse_time_dense_product(BlockSparseMatrixView(m_lhs), BlockVectorView(m_lhs, m_rhs), - tmpDest, alpha); - } - - private: - BlockSparseTimeDenseProduct& operator=(const BlockSparseTimeDenseProduct&); -}; - -template -class BlockSparseMatrix - : public SparseMatrixBase > { - public: - typedef Scalar_ Scalar; - typedef typename NumTraits::Real RealScalar; - typedef StorageIndex_ StorageIndex; - typedef - typename internal::ref_selector >::type - Nested; - - enum { - Options = Options_, - Flags = Options, - BlockSize = _BlockAtCompileTime, - RowsAtCompileTime = Dynamic, - ColsAtCompileTime = Dynamic, - MaxRowsAtCompileTime = Dynamic, - MaxColsAtCompileTime = Dynamic, - IsVectorAtCompileTime = 0, - IsColMajor = Flags & RowMajorBit ? 0 : 1 - }; - typedef Matrix BlockScalar; - typedef Matrix - BlockRealScalar; - typedef std::conditional_t<_BlockAtCompileTime == Dynamic, Scalar, BlockScalar> BlockScalarReturnType; - typedef BlockSparseMatrix PlainObject; - - public: - // Default constructor - BlockSparseMatrix() - : m_innerBSize(0), - m_outerBSize(0), - m_innerOffset(0), - m_outerOffset(0), - m_nonzerosblocks(0), - m_values(0), - m_blockPtr(0), - m_indices(0), - m_outerIndex(0), - m_blockSize(BlockSize) {} - - /** - * \brief Construct and resize - * - */ - BlockSparseMatrix(Index brow, Index bcol) - : m_innerBSize(IsColMajor ? brow : bcol), - m_outerBSize(IsColMajor ? bcol : brow), - m_innerOffset(0), - m_outerOffset(0), - m_nonzerosblocks(0), - m_values(0), - m_blockPtr(0), - m_indices(0), - m_outerIndex(0), - m_blockSize(BlockSize) {} - - /** - * \brief Copy-constructor - */ - BlockSparseMatrix(const BlockSparseMatrix& other) - : m_innerBSize(other.m_innerBSize), - m_outerBSize(other.m_outerBSize), - m_nonzerosblocks(other.m_nonzerosblocks), - m_nonzeros(other.m_nonzeros), - m_blockPtr(0), - m_blockSize(other.m_blockSize) { - // TODO: decide whether to allow copying between variable-size and fixed-size blocks. - eigen_assert(m_blockSize == BlockSize && " CAN NOT COPY BETWEEN FIXED-SIZE AND VARIABLE-SIZE BLOCKS"); - - std::copy(other.m_innerOffset, other.m_innerOffset + m_innerBSize + 1, m_innerOffset); - std::copy(other.m_outerOffset, other.m_outerOffset + m_outerBSize + 1, m_outerOffset); - std::copy(other.m_values, other.m_values + m_nonzeros, m_values); - - if (m_blockSize != Dynamic) std::copy(other.m_blockPtr, other.m_blockPtr + m_nonzerosblocks, m_blockPtr); - - std::copy(other.m_indices, other.m_indices + m_nonzerosblocks, m_indices); - std::copy(other.m_outerIndex, other.m_outerIndex + m_outerBSize, m_outerIndex); - } - - friend void swap(BlockSparseMatrix& first, BlockSparseMatrix& second) { - std::swap(first.m_innerBSize, second.m_innerBSize); - std::swap(first.m_outerBSize, second.m_outerBSize); - std::swap(first.m_innerOffset, second.m_innerOffset); - std::swap(first.m_outerOffset, second.m_outerOffset); - std::swap(first.m_nonzerosblocks, second.m_nonzerosblocks); - std::swap(first.m_nonzeros, second.m_nonzeros); - std::swap(first.m_values, second.m_values); - std::swap(first.m_blockPtr, second.m_blockPtr); - std::swap(first.m_indices, second.m_indices); - std::swap(first.m_outerIndex, second.m_outerIndex); - std::swap(first.m_BlockSize, second.m_blockSize); - } - - BlockSparseMatrix& operator=(BlockSparseMatrix other) { - // Copy-and-swap paradigm ... avoid leaked data if thrown - swap(*this, other); - return *this; - } - - // Destructor - ~BlockSparseMatrix() { - delete[] m_outerIndex; - delete[] m_innerOffset; - delete[] m_outerOffset; - delete[] m_indices; - delete[] m_blockPtr; - delete[] m_values; - } - - /** - * \brief Constructor from a sparse matrix - * - */ - template - inline BlockSparseMatrix(const MatrixType& spmat) : m_blockSize(BlockSize) { - EIGEN_STATIC_ASSERT((m_blockSize != Dynamic), THIS_METHOD_IS_ONLY_FOR_FIXED_SIZE); - - *this = spmat; - } - - /** - * \brief Assignment from a sparse matrix with the same storage order - * - * Convert from a sparse matrix to block sparse matrix. - * \warning Before calling this function, it is necessary to call - * either setBlockLayout() (matrices with variable-size blocks) - * or setBlockSize() (for fixed-size blocks). - */ - template - inline BlockSparseMatrix& operator=(const MatrixType& spmat) { - eigen_assert((m_innerBSize != 0 && m_outerBSize != 0) && - "Trying to assign to a zero-size matrix, call resize() first"); - eigen_assert(((MatrixType::Options & RowMajorBit) != IsColMajor) && "Wrong storage order"); - typedef SparseMatrix MatrixPatternType; - MatrixPatternType blockPattern(blockRows(), blockCols()); - m_nonzeros = 0; - - // First, compute the number of nonzero blocks and their locations - for (StorageIndex bj = 0; bj < m_outerBSize; ++bj) { - // Browse each outer block and compute the structure - std::vector nzblocksFlag(m_innerBSize, false); // Record the existing blocks - blockPattern.startVec(bj); - for (StorageIndex j = blockOuterIndex(bj); j < blockOuterIndex(bj + 1); ++j) { - typename MatrixType::InnerIterator it_spmat(spmat, j); - for (; it_spmat; ++it_spmat) { - StorageIndex bi = innerToBlock(it_spmat.index()); // Index of the current nonzero block - if (!nzblocksFlag[bi]) { - // Save the index of this nonzero block - nzblocksFlag[bi] = true; - blockPattern.insertBackByOuterInnerUnordered(bj, bi) = true; - // Compute the total number of nonzeros (including explicit zeros in blocks) - m_nonzeros += blockOuterSize(bj) * blockInnerSize(bi); - } - } - } // end current outer block - } - blockPattern.finalize(); - - // Allocate the internal arrays - setBlockStructure(blockPattern); - - for (StorageIndex nz = 0; nz < m_nonzeros; ++nz) m_values[nz] = Scalar(0); - for (StorageIndex bj = 0; bj < m_outerBSize; ++bj) { - // Now copy the values - for (StorageIndex j = blockOuterIndex(bj); j < blockOuterIndex(bj + 1); ++j) { - // Browse the outer block column by column (for column-major matrices) - typename MatrixType::InnerIterator it_spmat(spmat, j); - for (; it_spmat; ++it_spmat) { - StorageIndex idx = 0; // Position of this block in the column block - StorageIndex bi = innerToBlock(it_spmat.index()); // Index of the current nonzero block - // Go to the inner block where this element belongs to - while (bi > m_indices[m_outerIndex[bj] + idx]) ++idx; // Not expensive for ordered blocks - StorageIndex idxVal; // Get the right position in the array of values for this element - if (m_blockSize == Dynamic) { - // Offset from all blocks before ... - idxVal = m_blockPtr[m_outerIndex[bj] + idx]; - // ... and offset inside the block - idxVal += (j - blockOuterIndex(bj)) * blockOuterSize(bj) + it_spmat.index() - m_innerOffset[bi]; - } else { - // All blocks before - idxVal = (m_outerIndex[bj] + idx) * m_blockSize * m_blockSize; - // inside the block - idxVal += (j - blockOuterIndex(bj)) * m_blockSize + (it_spmat.index() % m_blockSize); - } - // Insert the value - m_values[idxVal] = it_spmat.value(); - } // end of this column - } // end of this block - } // end of this outer block - - return *this; - } - - /** - * \brief Set the nonzero block pattern of the matrix - * - * Given a sparse matrix describing the nonzero block pattern, - * this function prepares the internal pointers for values. - * After calling this function, any *nonzero* block (bi, bj) can be set - * with a simple call to coeffRef(bi,bj). - * - * - * \warning Before calling this function, it is necessary to call - * either setBlockLayout() (matrices with variable-size blocks) - * or setBlockSize() (for fixed-size blocks). - * - * \param blockPattern Sparse matrix of boolean elements describing the block structure - * - * \sa setBlockLayout() \sa setBlockSize() - */ - template - void setBlockStructure(const MatrixType& blockPattern) { - resize(blockPattern.rows(), blockPattern.cols()); - reserve(blockPattern.nonZeros()); - - // Browse the block pattern and set up the various pointers - m_outerIndex[0] = 0; - if (m_blockSize == Dynamic) m_blockPtr[0] = 0; - for (StorageIndex nz = 0; nz < m_nonzeros; ++nz) m_values[nz] = Scalar(0); - for (StorageIndex bj = 0; bj < m_outerBSize; ++bj) { - // Browse each outer block - - // First, copy and save the indices of nonzero blocks - // FIXME : find a way to avoid this ... - std::vector nzBlockIdx; - typename MatrixType::InnerIterator it(blockPattern, bj); - for (; it; ++it) { - nzBlockIdx.push_back(it.index()); - } - std::sort(nzBlockIdx.begin(), nzBlockIdx.end()); - - // Now, fill block indices and (eventually) pointers to blocks - for (StorageIndex idx = 0; idx < nzBlockIdx.size(); ++idx) { - StorageIndex offset = m_outerIndex[bj] + idx; // offset in m_indices - m_indices[offset] = nzBlockIdx[idx]; - if (m_blockSize == Dynamic) - m_blockPtr[offset] = m_blockPtr[offset - 1] + blockInnerSize(nzBlockIdx[idx]) * blockOuterSize(bj); - // There is no blockPtr for fixed-size blocks; not needed. - } - // Save the pointer to the next outer block - m_outerIndex[bj + 1] = m_outerIndex[bj] + nzBlockIdx.size(); - } - } - - /** - * \brief Set the number of rows and columns blocks - */ - inline void resize(Index brow, Index bcol) { - m_innerBSize = IsColMajor ? brow : bcol; - m_outerBSize = IsColMajor ? bcol : brow; - } - - /** - * \brief set the block size at runtime for fixed-size block layout - * - * Call this only for fixed-size blocks - */ - inline void setBlockSize(Index blockSize) { m_blockSize = blockSize; } - - /** - * \brief Set the row and column block layouts, - * - * This function set the size of each row and column block. - * So this function should be used only for blocks with variable size. - * \param rowBlocks : Number of rows per row block - * \param colBlocks : Number of columns per column block - * \sa resize(), setBlockSize() - */ - inline void setBlockLayout(const VectorXi& rowBlocks, const VectorXi& colBlocks) { - const VectorXi& innerBlocks = IsColMajor ? rowBlocks : colBlocks; - const VectorXi& outerBlocks = IsColMajor ? colBlocks : rowBlocks; - eigen_assert(m_innerBSize == innerBlocks.size() && "CHECK THE NUMBER OF ROW OR COLUMN BLOCKS"); - eigen_assert(m_outerBSize == outerBlocks.size() && "CHECK THE NUMBER OF ROW OR COLUMN BLOCKS"); - m_outerBSize = outerBlocks.size(); - // starting index of blocks... cumulative sums - m_innerOffset = new StorageIndex[m_innerBSize + 1]; - m_outerOffset = new StorageIndex[m_outerBSize + 1]; - m_innerOffset[0] = 0; - m_outerOffset[0] = 0; - std::partial_sum(&innerBlocks[0], &innerBlocks[m_innerBSize - 1] + 1, &m_innerOffset[1]); - std::partial_sum(&outerBlocks[0], &outerBlocks[m_outerBSize - 1] + 1, &m_outerOffset[1]); - - // Compute the total number of nonzeros - m_nonzeros = 0; - for (StorageIndex bj = 0; bj < m_outerBSize; ++bj) - for (StorageIndex bi = 0; bi < m_innerBSize; ++bi) m_nonzeros += outerBlocks[bj] * innerBlocks[bi]; - } - - /** - * \brief Allocate the internal array of pointers to blocks and their inner indices - * - * \note For fixed-size blocks, call setBlockSize() to set the block. - * And for variable-size blocks, call setBlockLayout() before using this function - * - * \param nonzerosblocks Number of nonzero blocks. The total number of nonzeros - * is computed in setBlockLayout() for variable-size blocks - * \sa setBlockSize() - */ - inline void reserve(const Index nonzerosblocks) { - eigen_assert((m_innerBSize != 0 && m_outerBSize != 0) && - "TRYING TO RESERVE ZERO-SIZE MATRICES, CALL resize() first"); - - // FIXME: Should free if already allocated. - m_outerIndex = new StorageIndex[m_outerBSize + 1]; - - m_nonzerosblocks = nonzerosblocks; - if (m_blockSize != Dynamic) { - m_nonzeros = nonzerosblocks * (m_blockSize * m_blockSize); - m_blockPtr = 0; - } else { - // m_nonzeros is already computed in setBlockLayout() - m_blockPtr = new StorageIndex[m_nonzerosblocks + 1]; - } - m_indices = new StorageIndex[m_nonzerosblocks + 1]; - m_values = new Scalar[m_nonzeros]; - } - - /** - * \brief Fill values in a matrix from a triplet list. - * - * Each triplet item has a block stored in an Eigen dense matrix. - * The InputIterator class should provide the functions row(), col() and value() - * - * \note For fixed-size blocks, call setBlockSize() before this function. - * - * FIXME: Do not accept duplicates. - */ - template - void setFromTriplets(const InputIterator& begin, const InputIterator& end) { - eigen_assert((m_innerBSize != 0 && m_outerBSize != 0) && "ZERO BLOCKS, PLEASE CALL resize() before"); - - /* First, sort the triplet list - * FIXME: This can be unnecessarily expensive since only the inner indices have to be sorted. - * The best approach is like in SparseMatrix::setFromTriplets() - */ - internal::TripletComp tripletcomp; - std::sort(begin, end, tripletcomp); - - /* Count the number of rows and column blocks, - * and the number of nonzero blocks per outer dimension - */ - VectorXi rowBlocks(m_innerBSize); // Size of each block row - VectorXi colBlocks(m_outerBSize); // Size of each block column - rowBlocks.setZero(); - colBlocks.setZero(); - VectorXi nzblock_outer(m_outerBSize); // Number of nz blocks per outer vector - VectorXi nz_outer(m_outerBSize); // Number of nz per outer vector...for variable-size blocks - nzblock_outer.setZero(); - nz_outer.setZero(); - for (InputIterator it(begin); it != end; ++it) { - eigen_assert(it->row() >= 0 && it->row() < this->blockRows() && it->col() >= 0 && it->col() < this->blockCols()); - eigen_assert((it->value().rows() == it->value().cols() && (it->value().rows() == m_blockSize)) || - (m_blockSize == Dynamic)); - - if (m_blockSize == Dynamic) { - eigen_assert((rowBlocks[it->row()] == 0 || rowBlocks[it->row()] == it->value().rows()) && - "NON CORRESPONDING SIZES FOR ROW BLOCKS"); - eigen_assert((colBlocks[it->col()] == 0 || colBlocks[it->col()] == it->value().cols()) && - "NON CORRESPONDING SIZES FOR COLUMN BLOCKS"); - rowBlocks[it->row()] = it->value().rows(); - colBlocks[it->col()] = it->value().cols(); - } - nz_outer(IsColMajor ? it->col() : it->row()) += it->value().rows() * it->value().cols(); - nzblock_outer(IsColMajor ? it->col() : it->row())++; - } - // Allocate member arrays - if (m_blockSize == Dynamic) setBlockLayout(rowBlocks, colBlocks); - StorageIndex nzblocks = nzblock_outer.sum(); - reserve(nzblocks); - - // Temporary markers - VectorXi block_id(m_outerBSize); // To be used as a block marker during insertion - - // Setup outer index pointers and markers - m_outerIndex[0] = 0; - if (m_blockSize == Dynamic) m_blockPtr[0] = 0; - for (StorageIndex bj = 0; bj < m_outerBSize; ++bj) { - m_outerIndex[bj + 1] = m_outerIndex[bj] + nzblock_outer(bj); - block_id(bj) = m_outerIndex[bj]; - if (m_blockSize == Dynamic) { - m_blockPtr[m_outerIndex[bj + 1]] = m_blockPtr[m_outerIndex[bj]] + nz_outer(bj); - } - } - - // Fill the matrix - for (InputIterator it(begin); it != end; ++it) { - StorageIndex outer = IsColMajor ? it->col() : it->row(); - StorageIndex inner = IsColMajor ? it->row() : it->col(); - m_indices[block_id(outer)] = inner; - StorageIndex block_size = it->value().rows() * it->value().cols(); - StorageIndex nz_marker = blockPtr(block_id[outer]); - memcpy(&(m_values[nz_marker]), it->value().data(), block_size * sizeof(Scalar)); - if (m_blockSize == Dynamic) { - m_blockPtr[block_id(outer) + 1] = m_blockPtr[block_id(outer)] + block_size; - } - block_id(outer)++; - } - } - - /** - * \returns the number of rows - */ - inline Index rows() const { return (IsColMajor ? innerSize() : outerSize()); } - - /** - * \returns the number of cols - */ - inline Index cols() const { return (IsColMajor ? outerSize() : innerSize()); } - - inline Index innerSize() const { - if (m_blockSize == Dynamic) - return m_innerOffset[m_innerBSize]; - else - return (m_innerBSize * m_blockSize); - } - - inline Index outerSize() const { - if (m_blockSize == Dynamic) - return m_outerOffset[m_outerBSize]; - else - return (m_outerBSize * m_blockSize); - } - /** \returns the number of rows grouped by blocks */ - inline Index blockRows() const { return (IsColMajor ? m_innerBSize : m_outerBSize); } - /** \returns the number of columns grouped by blocks */ - inline Index blockCols() const { return (IsColMajor ? m_outerBSize : m_innerBSize); } - - inline Index outerBlocks() const { return m_outerBSize; } - inline Index innerBlocks() const { return m_innerBSize; } - - /** \returns the block index where outer belongs to */ - inline Index outerToBlock(Index outer) const { - eigen_assert(outer < outerSize() && "OUTER INDEX OUT OF BOUNDS"); - - if (m_blockSize != Dynamic) return (outer / m_blockSize); // Integer division - - StorageIndex b_outer = 0; - while (m_outerOffset[b_outer] <= outer) ++b_outer; - return b_outer - 1; - } - /** \returns the block index where inner belongs to */ - inline Index innerToBlock(Index inner) const { - eigen_assert(inner < innerSize() && "OUTER INDEX OUT OF BOUNDS"); - - if (m_blockSize != Dynamic) return (inner / m_blockSize); // Integer division - - StorageIndex b_inner = 0; - while (m_innerOffset[b_inner] <= inner) ++b_inner; - return b_inner - 1; - } - - /** - *\returns a reference to the (i,j) block as an Eigen Dense Matrix - */ - Ref coeffRef(Index brow, Index bcol) { - eigen_assert(brow < blockRows() && "BLOCK ROW INDEX OUT OF BOUNDS"); - eigen_assert(bcol < blockCols() && "BLOCK nzblocksFlagCOLUMN OUT OF BOUNDS"); - - StorageIndex rsize = IsColMajor ? blockInnerSize(brow) : blockOuterSize(bcol); - StorageIndex csize = IsColMajor ? blockOuterSize(bcol) : blockInnerSize(brow); - StorageIndex inner = IsColMajor ? brow : bcol; - StorageIndex outer = IsColMajor ? bcol : brow; - StorageIndex offset = m_outerIndex[outer]; - while (offset < m_outerIndex[outer + 1] && m_indices[offset] != inner) offset++; - if (m_indices[offset] == inner) { - return Map(&(m_values[blockPtr(offset)]), rsize, csize); - } else { - // FIXME: The block does not exist; insert it. - eigen_assert("DYNAMIC INSERTION IS NOT YET SUPPORTED"); - } - } - - /** - * \returns the value of the (i,j) block as an Eigen Dense Matrix - */ - Map coeff(Index brow, Index bcol) const { - eigen_assert(brow < blockRows() && "BLOCK ROW INDEX OUT OF BOUNDS"); - eigen_assert(bcol < blockCols() && "BLOCK COLUMN OUT OF BOUNDS"); - - StorageIndex rsize = IsColMajor ? blockInnerSize(brow) : blockOuterSize(bcol); - StorageIndex csize = IsColMajor ? blockOuterSize(bcol) : blockInnerSize(brow); - StorageIndex inner = IsColMajor ? brow : bcol; - StorageIndex outer = IsColMajor ? bcol : brow; - StorageIndex offset = m_outerIndex[outer]; - while (offset < m_outerIndex[outer + 1] && m_indices[offset] != inner) offset++; - if (m_indices[offset] == inner) { - return Map(&(m_values[blockPtr(offset)]), rsize, csize); - } else - eigen_assert("NOT YET SUPPORTED"); - } - - // Block Matrix times vector product - template - BlockSparseTimeDenseProduct operator*(const VecType& lhs) const { - return BlockSparseTimeDenseProduct(*this, lhs); - } - - /** \returns the number of nonzero blocks */ - inline Index nonZerosBlocks() const { return m_nonzerosblocks; } - /** \returns the total number of nonzero elements, including eventual explicit zeros in blocks */ - inline Index nonZeros() const { return m_nonzeros; } - - inline BlockScalarReturnType* valuePtr() { return static_cast(m_values); } - inline StorageIndex* innerIndexPtr() { return m_indices; } - inline const StorageIndex* innerIndexPtr() const { return m_indices; } - inline StorageIndex* outerIndexPtr() { return m_outerIndex; } - inline const StorageIndex* outerIndexPtr() const { return m_outerIndex; } - - /** \brief for compatibility purposes with the SparseMatrix class */ - inline bool isCompressed() const { return true; } - /** - * \returns the starting index of the bi row block - */ - inline Index blockRowsIndex(Index bi) const { return IsColMajor ? blockInnerIndex(bi) : blockOuterIndex(bi); } - - /** - * \returns the starting index of the bj col block - */ - inline Index blockColsIndex(Index bj) const { return IsColMajor ? blockOuterIndex(bj) : blockInnerIndex(bj); } - - inline Index blockOuterIndex(Index bj) const { - return (m_blockSize == Dynamic) ? m_outerOffset[bj] : (bj * m_blockSize); - } - inline Index blockInnerIndex(Index bi) const { - return (m_blockSize == Dynamic) ? m_innerOffset[bi] : (bi * m_blockSize); - } - - // Not needed. - inline Index blockInnerSize(Index bi) const { - return (m_blockSize == Dynamic) ? (m_innerOffset[bi + 1] - m_innerOffset[bi]) : m_blockSize; - } - inline Index blockOuterSize(Index bj) const { - return (m_blockSize == Dynamic) ? (m_outerOffset[bj + 1] - m_outerOffset[bj]) : m_blockSize; - } - - /** - * \brief Browse the matrix by outer index - */ - class InnerIterator; // Browse column by column - - /** - * \brief Browse the matrix by block outer index - */ - class BlockInnerIterator; // Browse block by block - - friend std::ostream& operator<<(std::ostream& s, const BlockSparseMatrix& m) { - for (StorageIndex j = 0; j < m.outerBlocks(); ++j) { - BlockInnerIterator itb(m, j); - for (; itb; ++itb) { - s << "(" << itb.row() << ", " << itb.col() << ")\n"; - s << itb.value() << "\n"; - } - } - s << std::endl; - return s; - } - - /** - * \returns the starting position of the block \p id in the array of values - */ - Index blockPtr(Index id) const { - if (m_blockSize == Dynamic) - return m_blockPtr[id]; - else - return id * m_blockSize * m_blockSize; - } - - protected: - // To be implemented - // Insert a block at a particular location... need to make a room for that - Map insert(Index brow, Index bcol); - - Index m_innerBSize; // Number of block rows - Index m_outerBSize; // Number of block columns - StorageIndex* m_innerOffset; // Starting index of each inner block (size m_innerBSize+1) - StorageIndex* m_outerOffset; // Starting index of each outer block (size m_outerBSize+1) - Index m_nonzerosblocks; // Total nonzeros blocks (lower than m_innerBSize x m_outerBSize) - Index m_nonzeros; // Total nonzeros elements - Scalar* m_values; // Values stored block column after block column (size m_nonzeros) - StorageIndex* m_blockPtr; // Pointer to the beginning of each block in m_values, size m_nonzerosblocks ... null for - // fixed-size blocks - StorageIndex* m_indices; // Inner block indices, size m_nonzerosblocks ... OK - StorageIndex* m_outerIndex; // Starting pointer of each block column in m_indices (size m_outerBSize)... OK - Index m_blockSize; // Size of a block for fixed-size blocks, otherwise -1 -}; - -template -class BlockSparseMatrix::BlockInnerIterator { - public: - enum { Flags = Options_ }; - - BlockInnerIterator(const BlockSparseMatrix& mat, const Index outer) - : m_mat(mat), m_outer(outer), m_id(mat.m_outerIndex[outer]), m_end(mat.m_outerIndex[outer + 1]) {} - - inline BlockInnerIterator& operator++() { - m_id++; - return *this; - } - - inline const Map value() const { - return Map(&(m_mat.m_values[m_mat.blockPtr(m_id)]), rows(), cols()); - } - inline Map valueRef() { - return Map(&(m_mat.m_values[m_mat.blockPtr(m_id)]), rows(), cols()); - } - // Block inner index - inline Index index() const { return m_mat.m_indices[m_id]; } - inline Index outer() const { return m_outer; } - // block row index - inline Index row() const { return index(); } - // block column index - inline Index col() const { return outer(); } - // FIXME: Number of rows in the current block. - inline Index rows() const { - return (m_mat.m_blockSize == Dynamic) ? (m_mat.m_innerOffset[index() + 1] - m_mat.m_innerOffset[index()]) - : m_mat.m_blockSize; - } - // Number of columns in the current block ... - inline Index cols() const { - return (m_mat.m_blockSize == Dynamic) ? (m_mat.m_outerOffset[m_outer + 1] - m_mat.m_outerOffset[m_outer]) - : m_mat.m_blockSize; - } - inline operator bool() const { return (m_id < m_end); } - - protected: - const BlockSparseMatrix& m_mat; - const Index m_outer; - Index m_id; - Index m_end; -}; - -template -class BlockSparseMatrix::InnerIterator { - public: - InnerIterator(const BlockSparseMatrix& mat, Index outer) - : m_mat(mat), - m_outerB(mat.outerToBlock(outer)), - m_outer(outer), - itb(mat, mat.outerToBlock(outer)), - m_offset(outer - mat.blockOuterIndex(m_outerB)) { - if (itb) { - m_id = m_mat.blockInnerIndex(itb.index()); - m_start = m_id; - m_end = m_mat.blockInnerIndex(itb.index() + 1); - } - } - inline InnerIterator& operator++() { - m_id++; - if (m_id >= m_end) { - ++itb; - if (itb) { - m_id = m_mat.blockInnerIndex(itb.index()); - m_start = m_id; - m_end = m_mat.blockInnerIndex(itb.index() + 1); - } - } - return *this; - } - inline const Scalar& value() const { return itb.value().coeff(m_id - m_start, m_offset); } - inline Scalar& valueRef() { return itb.valueRef().coeff(m_id - m_start, m_offset); } - inline Index index() const { return m_id; } - inline Index outer() const { return m_outer; } - inline Index col() const { return outer(); } - inline Index row() const { return index(); } - inline operator bool() const { return itb; } - - protected: - const BlockSparseMatrix& m_mat; - const Index m_outer; - const Index m_outerB; - BlockInnerIterator itb; // Iterator through the blocks - const Index m_offset; // Position of this column in the block - Index m_start; // starting inner index of this block - Index m_id; // current inner index in the block - Index m_end; // starting inner index of the next block -}; -} // end namespace Eigen - -#endif // EIGEN_SPARSEBLOCKMATRIX_H