diff --git a/Eigen/src/SparseCore/BlockSparseMatrix.h b/Eigen/src/SparseCore/BlockSparseMatrix.h index ee5ad8733..44b1dc287 100644 --- a/Eigen/src/SparseCore/BlockSparseMatrix.h +++ b/Eigen/src/SparseCore/BlockSparseMatrix.h @@ -1350,21 +1350,26 @@ class BlockSparseTriangularView { 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); + // BSM befriends BlockSparseTriangularView, but GCC <= 11 and MSVC do not + // extend that friendship to the body of a friend function defined inside + // the view (CWG 1699), so only BSM's public interface may be used here. + const BSM& m = tri.m_matrix; + const StorageIndex* outerPtr = m.outerIndexPtr(); + const StorageIndex* innerPtr = m.innerIndexPtr(); + ResultType result = ResultType::Zero(lhs.rows(), m.cols()); + for (Index out = 0; out < m.blockOuterSize(); ++out) { + for (Index id = outerPtr[out]; id < outerPtr[out + 1]; ++id) { + Index inner = innerPtr[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(); + lhs.template middleCols(bi * BlockRows) * m.blockRef(id).template triangularView(); else result.template middleCols(bj * BlockCols).noalias() += - lhs.template middleCols(bi * BlockRows) * tri.m_matrix.blockRef(id); + lhs.template middleCols(bi * BlockRows) * m.blockRef(id); } } return result; diff --git a/test/block_sparse_matrix.cpp b/test/block_sparse_matrix.cpp index f20cae826..f4150694e 100644 --- a/test/block_sparse_matrix.cpp +++ b/test/block_sparse_matrix.cpp @@ -39,6 +39,33 @@ BlockSparseMatrix denseToBl return bsm; } +template +struct BlockSparseIdentityTester { + template + static void run(int, int, const BlockSparseMatrix&, + const Matrix&) {} +}; + +template <> +struct BlockSparseIdentityTester { + template + static void run(int bRows, int bCols, const BlockSparseMatrix& A, + const Matrix& dA) { + using BSM = BlockSparseMatrix; + using DenseMat = Matrix; + + const int rows = bRows * BlockRows; + const int cols = bCols * 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); + } + } +}; + // --------------------------------------------------------------------------- // Core test driver templated on block size, storage order, and scalar type // --------------------------------------------------------------------------- @@ -167,15 +194,8 @@ void test_block_sparse(int bRows, int bCols) { } // ---- 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); - } - } + BlockSparseIdentityTester::template run( + bRows, bCols, A, dA); // ---- setFromTriplets with duplicate blocks (accumulation) --------------- {