SparseCore: SparsityPatternRef view; pattern-only AMD ordering & ILUT prep libeigen/eigen!2489 Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>
diff --git a/Eigen/SparseCore b/Eigen/SparseCore index 6020e42..5d941e6 100644 --- a/Eigen/SparseCore +++ b/Eigen/SparseCore
@@ -59,6 +59,7 @@ #include "src/SparseCore/SparsePermutation.h" #include "src/SparseCore/SparseFuzzy.h" #include "src/SparseCore/SparseSolverBase.h" +#include "src/SparseCore/SparsityPatternRef.h" // IWYU pragma: end_exports #include "src/Core/util/ReenableStupidWarnings.h"
diff --git a/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h b/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h index 61e8d3c..5fbb5f2 100644 --- a/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h +++ b/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h
@@ -195,20 +195,6 @@ template <typename MatrixType> Index computeRowMatching(const MatrixType& amat); - // Produce a column-major CSC sparsity pattern for `amat` (integers only — - // the scalar values are never read or copied). When `amat` is already a - // compressed column-major SparseMatrix, `outer`/`inner` point directly into - // its index storage; otherwise the pattern is materialized into the local - // `outer_buf`/`inner_buf` arrays and pointers are set into them. - static void patternColMajor(const SparseMatrix<Scalar, ColMajor, StorageIndex>& amat, - Matrix<StorageIndex, Dynamic, 1>& outer_buf, Matrix<StorageIndex, Dynamic, 1>& inner_buf, - const StorageIndex*& outer, const StorageIndex*& inner); - - template <typename MatrixType> - static void patternColMajor(const MatrixType& amat, Matrix<StorageIndex, Dynamic, 1>& outer_buf, - Matrix<StorageIndex, Dynamic, 1>& inner_buf, const StorageIndex*& outer, - const StorageIndex*& inner); - protected: FactorType m_lu; RealScalar m_droptol; @@ -262,58 +248,6 @@ return m_lu.template triangularView<Upper>(); } -// Specialization: amat is already a column-major SparseMatrix. -// Share its index storage directly when compressed; otherwise materialize the -// indices (without copying any scalar values) into outer_buf/inner_buf. -template <typename Scalar, typename StorageIndex> -void IncompleteLUT<Scalar, StorageIndex>::patternColMajor(const SparseMatrix<Scalar, ColMajor, StorageIndex>& amat, - Matrix<StorageIndex, Dynamic, 1>& outer_buf, - Matrix<StorageIndex, Dynamic, 1>& inner_buf, - const StorageIndex*& outer, const StorageIndex*& inner) { - if (amat.isCompressed()) { - outer = amat.outerIndexPtr(); - inner = amat.innerIndexPtr(); - return; - } - const Index n = amat.cols(); - const StorageIndex* a_outer = amat.outerIndexPtr(); - const StorageIndex* a_inner_nz = amat.innerNonZeroPtr(); - const StorageIndex* a_inner = amat.innerIndexPtr(); - outer_buf.resize(n + 1); - outer_buf(0) = 0; - for (Index j = 0; j < n; ++j) outer_buf(j + 1) = outer_buf(j) + a_inner_nz[j]; - inner_buf.resize(outer_buf(n)); - for (Index j = 0; j < n; ++j) { - const StorageIndex* src = a_inner + a_outer[j]; - std::copy(src, src + a_inner_nz[j], inner_buf.data() + outer_buf(j)); - } - outer = outer_buf.data(); - inner = inner_buf.data(); -} - -// Generic fallback: any other sparse input (row-major, expressions). Build a -// column-major pattern via inner iterators — no scalar values are read. -template <typename Scalar, typename StorageIndex> -template <typename MatrixType_> -void IncompleteLUT<Scalar, StorageIndex>::patternColMajor(const MatrixType_& amat, - Matrix<StorageIndex, Dynamic, 1>& outer_buf, - Matrix<StorageIndex, Dynamic, 1>& inner_buf, - const StorageIndex*& outer, const StorageIndex*& inner) { - using internal::convert_index; - const Index n = amat.cols(); - outer_buf.setZero(n + 1); - for (Index i = 0; i < amat.outerSize(); ++i) - for (typename MatrixType_::InnerIterator it(amat, i); it; ++it) ++outer_buf(it.col() + 1); - for (Index j = 0; j < n; ++j) outer_buf(j + 1) += outer_buf(j); - inner_buf.resize(outer_buf(n)); - Matrix<StorageIndex, Dynamic, 1> head = outer_buf.head(n); - for (Index i = 0; i < amat.outerSize(); ++i) - for (typename MatrixType_::InnerIterator it(amat, i); it; ++it) - inner_buf(head(it.col())++) = convert_index<StorageIndex>(it.row()); - outer = outer_buf.data(); - inner = inner_buf.data(); -} - // Compute a row permutation m_Pr such that (m_Pr * amat) has a structurally // nonzero diagonal wherever one exists. Returns the number of matched columns. // Uses a maximum bipartite cardinality matching on the sparsity pattern, with @@ -324,14 +258,15 @@ Index IncompleteLUT<Scalar, StorageIndex>::computeRowMatching(const MatrixType_& amat) { using internal::convert_index; const Index n = amat.rows(); - // We only need the column-major sparsity pattern; never read scalar values. - // Share amat's index storage when it is already a compressed column-major - // SparseMatrix, otherwise build a value-free pattern into local arrays. + // We only need amat's column-major sparsity pattern; never read scalar + // values. The pattern view aliases amat's index storage when amat is + // already a column-major SparseMatrix, and otherwise materializes a CSC + // pattern into the scratch buffers. Matrix<StorageIndex, Dynamic, 1> outer_buf; Matrix<StorageIndex, Dynamic, 1> inner_buf; - const StorageIndex* outer = nullptr; - const StorageIndex* inner = nullptr; - patternColMajor(amat, outer_buf, inner_buf, outer, inner); + internal::SparsityPatternRef<StorageIndex> pat = internal::make_col_major_pattern_ref(amat, outer_buf, inner_buf); + const StorageIndex* outer = pat.outer; + const StorageIndex* inner = pat.inner; const StorageIndex kUnmatched = StorageIndex(-1); // match_row[j] = original row matched to column j; match_col[i] = column matched to row i. @@ -343,7 +278,8 @@ // the same analysis is reusable for any matrix sharing this stored pattern. // Phase 1: greedy diagonal preference. for (Index j = 0; j < n; ++j) { - for (Index k = outer[j]; k < outer[j + 1]; ++k) { + const Index col_end = outer[j] + pat.nonZeros(j); + for (Index k = outer[j]; k < col_end; ++k) { if (Index(inner[k]) == j) { match_row[j] = convert_index<StorageIndex>(j); match_col[j] = convert_index<StorageIndex>(j); @@ -354,7 +290,8 @@ // Phase 2: greedy off-diagonal pickup of any free row. for (Index j = 0; j < n; ++j) { if (match_row[j] != kUnmatched) continue; - for (Index k = outer[j]; k < outer[j + 1]; ++k) { + const Index col_end = outer[j] + pat.nonZeros(j); + for (Index k = outer[j]; k < col_end; ++k) { Index i = inner[k]; if (match_col[i] == kUnmatched) { match_row[j] = convert_index<StorageIndex>(i); @@ -387,7 +324,7 @@ while (!stack_col.empty()) { Index j = stack_col.back(); Index pos = stack_pos.back(); - Index col_end = outer[j + 1]; + Index col_end = outer[j] + pat.nonZeros(j); bool advanced = false; while (pos < col_end) { @@ -463,16 +400,21 @@ computeRowMatching(amat); // 2. Compute the Fill-reducing permutation on the row-permuted matrix. - // Since ILUT does not perform any numerical pivoting, - // it is highly preferable to keep the diagonal through symmetric permutations. - // To this end, let's symmetrize the pattern and perform AMD on it. - SparseMatrix<Scalar, ColMajor, StorageIndex> mat1 = m_Pr * amat; - SparseMatrix<Scalar, ColMajor, StorageIndex> mat2 = mat1.transpose(); - // FIXME: for a nearly symmetric pattern, mat2+mat1 is appropriate; - // for a highly non-symmetric pattern, mat2*mat1 should be preferred. - SparseMatrix<Scalar, ColMajor, StorageIndex> AtA = mat2 + mat1; + // Since ILUT does not perform any numerical pivoting, it is highly + // preferable to keep the diagonal through symmetric permutations. AMD + // computes a fill-reducing ordering for a symmetric matrix and only reads + // the sparsity pattern; build a value-free, row-permuted representation + // (1-byte placeholder Scalar, indices remapped through m_Pr) and feed that + // to AMDOrdering, avoiding the previous mat1/mat2/AtA value copies. + SparseMatrix<signed char, ColMajor, StorageIndex> permuted_pattern; + { + Matrix<StorageIndex, Dynamic, 1> outer_buf; + Matrix<StorageIndex, Dynamic, 1> inner_buf; + internal::SparsityPatternRef<StorageIndex> pat = internal::make_col_major_pattern_ref(amat, outer_buf, inner_buf); + internal::materialize_col_major_pattern(pat, m_Pr.indices().data(), permuted_pattern); + } AMDOrdering<StorageIndex> ordering; - ordering(AtA, m_P); + ordering(permuted_pattern, m_P); m_Pinv = m_P.inverse(); // cache the inverse permutation // Cache the composition m_Pinv * m_Pr so _solve_impl applies a single // permutation to the RHS instead of two.
diff --git a/Eigen/src/OrderingMethods/Ordering.h b/Eigen/src/OrderingMethods/Ordering.h index 8ede4de..8dcc0aa 100644 --- a/Eigen/src/OrderingMethods/Ordering.h +++ b/Eigen/src/OrderingMethods/Ordering.h
@@ -20,14 +20,14 @@ /** \internal * \ingroup OrderingMethods_Module - * \param[in] A the input non-symmetric matrix - * \param[out] symmat the symmetric pattern A^T+A from the input matrix \a A. - * FIXME: only the sparsity pattern should be used here; values should be ignored. + * Build the symmetric sparsity pattern \c symmat = pattern(A^T + A) from a + * column-major input \a A. Only the sparsity pattern is read or written — + * scalar values are placeholders and are not meaningful. */ template <typename MatrixType> void ordering_helper_at_plus_a(const MatrixType& A, MatrixType& symmat) { MatrixType C; - C = A.transpose(); // NOTE: Could be costly + C = A.transpose(); for (int i = 0; i < C.rows(); i++) { for (typename MatrixType::InnerIterator it(C, i); it; ++it) it.valueRef() = typename MatrixType::Scalar(0); } @@ -40,7 +40,8 @@ * \class AMDOrdering * * Functor computing the \em approximate \em minimum \em degree ordering - * If the matrix is not structurally symmetric, an ordering of A^T+A is computed + * If the matrix is not structurally symmetric, an ordering of A^T+A is computed. + * Only the sparsity pattern of the input is read — scalar values are not. * \tparam StorageIndex The type of indices of the matrix * \sa COLAMDOrdering */ @@ -49,28 +50,47 @@ public: typedef PermutationMatrix<Dynamic, Dynamic, StorageIndex> PermutationType; - /** Compute the permutation vector from a sparse matrix - * This routine is much faster if the input matrix is column-major + /** Compute the permutation vector from a sparse matrix. + * Only the sparsity pattern of \a mat is read; scalar values are not. + * This routine is much faster if the input matrix is column-major. */ template <typename MatrixType> void operator()(const MatrixType& mat, PermutationType& perm) const { - // Compute the symmetric pattern - SparseMatrix<typename MatrixType::Scalar, ColMajor, StorageIndex> symm; - internal::ordering_helper_at_plus_a(mat, symm); - - // Call the AMD routine - // m_mat.prune(keep_diag()); + // AMD only reads the sparsity pattern. Project mat to a column-major + // SparseMatrix<signed char> (1-byte placeholder values; never read), then + // symmetrize and run minimum-degree ordering — independent of the source's + // Scalar type or storage order, this avoids the O(nnz) Scalar-value copy + // the previous implementation paid to build the symmetric pattern. + SparseMatrix<signed char, ColMajor, StorageIndex> A; + { + Matrix<StorageIndex, Dynamic, 1> outer_buf; + Matrix<StorageIndex, Dynamic, 1> inner_buf; + internal::SparsityPatternRef<StorageIndex> pat = internal::make_col_major_pattern_ref(mat, outer_buf, inner_buf); + internal::materialize_col_major_pattern(pat, A); + } + SparseMatrix<signed char, ColMajor, StorageIndex> symm; + internal::ordering_helper_at_plus_a(A, symm); internal::minimum_degree_ordering(symm, perm); } - /** Compute the permutation with a selfadjoint matrix */ + /** Compute the permutation with a selfadjoint matrix. + * Only the sparsity pattern is used; scalar values are not. + */ template <typename SrcType, unsigned int SrcUpLo> void operator()(const SparseSelfAdjointView<SrcType, SrcUpLo>& mat, PermutationType& perm) const { - SparseMatrix<typename SrcType::Scalar, ColMajor, StorageIndex> C; - C = mat; - - // Call the AMD routine - // m_mat.prune(keep_diag()); //Remove the diagonal elements + // Materialize the underlying triangle's pattern as a SparseMatrix<signed char> + // (works for any source Scalar including complex — no value read), then + // expand to a full symmetric SparseMatrix<signed char>. + SparseMatrix<signed char, ColMajor, StorageIndex> sc_src; + { + Matrix<StorageIndex, Dynamic, 1> outer_buf; + Matrix<StorageIndex, Dynamic, 1> inner_buf; + internal::SparsityPatternRef<StorageIndex> pat = + internal::make_col_major_pattern_ref(mat.matrix(), outer_buf, inner_buf); + internal::materialize_col_major_pattern(pat, sc_src); + } + SparseMatrix<signed char, ColMajor, StorageIndex> C; + C = sc_src.template selfadjointView<SrcUpLo>(); internal::minimum_degree_ordering(C, perm); } };
diff --git a/Eigen/src/SparseCore/SparsityPatternRef.h b/Eigen/src/SparseCore/SparsityPatternRef.h new file mode 100644 index 0000000..73e03b4 --- /dev/null +++ b/Eigen/src/SparseCore/SparsityPatternRef.h
@@ -0,0 +1,180 @@ +// 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/. + +#ifndef EIGEN_SPARSITY_PATTERN_REF_H +#define EIGEN_SPARSITY_PATTERN_REF_H + +// IWYU pragma: private +#include "./InternalHeaderCheck.h" + +namespace Eigen { +namespace internal { + +/** \internal + * \ingroup SparseCore_Module + * + * Non-owning view of a column-major sparsity pattern. The pattern is encoded + * as CSC index arrays: + * - \c outer (size \c outerSize+1 when compressed) gives the start of each + * column's row-index range in \c inner. + * - \c inner holds the row indices. + * - \c innerNonZero is non-null only for an uncompressed source: in that + * case column \c j's row indices live in + * [outer[j], outer[j] + innerNonZero[j]) instead of [outer[j], outer[j+1]). + * + * Use the helper \ref nonZeros to iterate uniformly over both layouts: + * \code + * const Index nz = pattern.nonZeros(j); + * for (Index k = pattern.outer[j], end = k + nz; k < end; ++k) { + * Index row = pattern.inner[k]; + * } + * \endcode + * + * Construct via \ref make_col_major_pattern_ref, which shares the source + * storage when possible (already column-major SparseMatrix) and otherwise + * materializes the pattern into caller-supplied scratch buffers without ever + * reading the source's scalar values. + */ +template <typename StorageIndex> +struct SparsityPatternRef { + const StorageIndex* outer = nullptr; + const StorageIndex* inner = nullptr; + // null ⇒ source is compressed: column j is [outer[j], outer[j+1]). + // otherwise: column j is [outer[j], outer[j] + innerNonZero[j]). + const StorageIndex* innerNonZero = nullptr; + Index outerSize = 0; // number of columns + Index innerSize = 0; // number of rows + + EIGEN_DEVICE_FUNC inline bool isCompressed() const { return innerNonZero == nullptr; } + + EIGEN_DEVICE_FUNC inline Index nonZeros(Index j) const { + return isCompressed() ? Index(outer[j + 1] - outer[j]) : Index(innerNonZero[j]); + } +}; + +/** \internal + * Build a column-major sparsity-pattern view of \c amat without copying any + * scalar values. The returned view aliases \c amat's index storage when + * \c amat is already a column-major \c SparseMatrix; otherwise the pattern is + * materialized into the caller-supplied \c outer_buf / \c inner_buf and the + * returned view points into them. + * + * The buffers must outlive every use of the returned view. + * + * Specialization for column-major \c SparseMatrix (the no-copy fast path). + */ +template <typename Scalar, typename StorageIndex> +SparsityPatternRef<StorageIndex> make_col_major_pattern_ref(const SparseMatrix<Scalar, ColMajor, StorageIndex>& amat, + Matrix<StorageIndex, Dynamic, 1>& /*outer_buf*/, + Matrix<StorageIndex, Dynamic, 1>& /*inner_buf*/) { + SparsityPatternRef<StorageIndex> p; + p.outer = amat.outerIndexPtr(); + p.inner = amat.innerIndexPtr(); + p.innerNonZero = amat.isCompressed() ? nullptr : amat.innerNonZeroPtr(); + p.outerSize = amat.cols(); + p.innerSize = amat.rows(); + return p; +} + +/** \internal + * Generic fallback for any other sparse expression (row-major, products, + * permutations, etc.). The pattern is materialized into \c outer_buf and + * \c inner_buf via a value-free counting transpose driven by the expression + * evaluator. Some evaluators may materialize internally. + */ +template <typename Derived> +SparsityPatternRef<typename Derived::StorageIndex> make_col_major_pattern_ref( + const SparseMatrixBase<Derived>& amat_base, Matrix<typename Derived::StorageIndex, Dynamic, 1>& outer_buf, + Matrix<typename Derived::StorageIndex, Dynamic, 1>& inner_buf) { + typedef typename Derived::StorageIndex StorageIndex; + const Derived& amat = amat_base.derived(); + internal::evaluator<Derived> amat_eval(amat); + const Index n_cols = amat.cols(); + const Index n_outer = amat.outerSize(); + + outer_buf.setZero(n_cols + 1); + for (Index i = 0; i < n_outer; ++i) + for (typename internal::evaluator<Derived>::InnerIterator it(amat_eval, i); it; ++it) ++outer_buf(it.col() + 1); + for (Index j = 0; j < n_cols; ++j) outer_buf(j + 1) += outer_buf(j); + inner_buf.resize(outer_buf(n_cols)); + Matrix<StorageIndex, Dynamic, 1> head = outer_buf.head(n_cols); + for (Index i = 0; i < n_outer; ++i) + for (typename internal::evaluator<Derived>::InnerIterator it(amat_eval, i); it; ++it) + inner_buf(head(it.col())++) = convert_index<StorageIndex>(it.row()); + + SparsityPatternRef<StorageIndex> p; + p.outer = outer_buf.data(); + p.inner = inner_buf.data(); + p.innerNonZero = nullptr; + p.outerSize = n_cols; + p.innerSize = amat.rows(); + return p; +} + +/** \internal + * Materialize a column-major sparsity pattern view as a compressed + * \c SparseMatrix with a 1-byte placeholder \c Scalar. Values are filled with + * a fixed nonzero sentinel — downstream pattern consumers (\c transpose(), + * sparse \c operator+, \c selfadjointView expansion) read coefficients even + * when the algorithm only cares about the pattern, so leaving them + * uninitialized would be UB. + * + * Intended for fill-reducing ordering algorithms (AMD, COLAMD) that read only + * the pattern. \b Precondition: \c out's storage must not alias \c pat — + * \c resize / \c resizeNonZeros may reallocate and invalidate the view. + * + * If \c row_perm is non-null, each inner row index is remapped through it + * (entry \c (i, j) of \c pat becomes \c (row_perm[i], j) in the result), + * implementing a left multiplication by a row permutation without copying any + * scalar values from the source. + * + * The result preserves SparseMatrix's invariant that inner indices are sorted + * within each column. + */ +template <typename StorageIndex> +void materialize_col_major_pattern(const SparsityPatternRef<StorageIndex>& pat, const StorageIndex* row_perm, + SparseMatrix<signed char, ColMajor, StorageIndex>& out) { + const Index n_outer = pat.outerSize; + const Index n_inner = pat.innerSize; + out.resize(n_inner, n_outer); + + Index total = 0; + for (Index j = 0; j < n_outer; ++j) total += pat.nonZeros(j); + out.resizeNonZeros(total); + + StorageIndex* o_outer = out.outerIndexPtr(); + StorageIndex* o_inner = out.innerIndexPtr(); + o_outer[0] = 0; + for (Index j = 0; j < n_outer; ++j) { + const Index nz = pat.nonZeros(j); + o_outer[j + 1] = convert_index<StorageIndex>(o_outer[j] + nz); + const StorageIndex* src = pat.inner + pat.outer[j]; + StorageIndex* dst = o_inner + o_outer[j]; + if (row_perm == nullptr) { + std::copy(src, src + nz, dst); + } else { + for (Index k = 0; k < nz; ++k) dst[k] = row_perm[src[k]]; + std::sort(dst, dst + nz); + } + } + // Fill values with a fixed nonzero sentinel so downstream consumers that + // read coefficients (transpose, sparse +, selfadjointView expansion) do not + // observe uninitialized memory. + std::fill_n(out.valuePtr(), total, static_cast<signed char>(1)); +} + +/** \internal Convenience overload — identity row permutation. */ +template <typename StorageIndex> +void materialize_col_major_pattern(const SparsityPatternRef<StorageIndex>& pat, + SparseMatrix<signed char, ColMajor, StorageIndex>& out) { + materialize_col_major_pattern(pat, static_cast<const StorageIndex*>(nullptr), out); +} + +} // namespace internal +} // namespace Eigen + +#endif // EIGEN_SPARSITY_PATTERN_REF_H
diff --git a/test/sparse_basic.cpp b/test/sparse_basic.cpp index d583d9d..5ea6843 100644 --- a/test/sparse_basic.cpp +++ b/test/sparse_basic.cpp
@@ -1109,6 +1109,56 @@ VERIFY_IS_EQUAL(vec.coeff(0), 0.0); } +template <typename Pattern, typename SparseMatrixType> +void verify_sparsity_pattern_ref_matches(const Pattern& pattern, const SparseMatrixType& expected) { + VERIFY_IS_EQUAL(pattern.outerSize, expected.cols()); + VERIFY_IS_EQUAL(pattern.innerSize, expected.rows()); + + Index pattern_nonzeros = 0; + for (Index j = 0; j < pattern.outerSize; ++j) { + const Index nz = pattern.nonZeros(j); + pattern_nonzeros += nz; + for (Index k = pattern.outer[j], end = k + nz; k < end; ++k) { + VERIFY(expected.coeff(pattern.inner[k], j) != typename SparseMatrixType::Scalar(0)); + } + } + VERIFY_IS_EQUAL(pattern_nonzeros, expected.nonZeros()); +} + +template <int> +void sparsity_pattern_ref_sparse_expressions() { + typedef std::complex<double> Scalar; + typedef SparseMatrix<Scalar, ColMajor, int> SparseMatrixType; + typedef SparseMatrix<double, ColMajor, int> RealSparseMatrixType; + typedef SparseMatrix<double, RowMajor, int> RowMajorSparseMatrixType; + typedef Matrix<int, Dynamic, 1> VectorI; + + SparseMatrixType a(3, 3); + std::vector<Triplet<Scalar, int>> triplets; + triplets.emplace_back(0, 0, Scalar(1.0, 1.0)); + triplets.emplace_back(1, 0, Scalar(2.0, -1.0)); + triplets.emplace_back(2, 1, Scalar(3.0, 2.0)); + triplets.emplace_back(2, 2, Scalar(4.0, -3.0)); + a.setFromTriplets(triplets.begin(), triplets.end()); + + PermutationMatrix<Dynamic, Dynamic, int> p(3); + p.indices() << 1, 2, 0; + + VectorI outer, inner; + SparseMatrixType permuted = p * a; + const internal::SparsityPatternRef<int> permuted_pattern = internal::make_col_major_pattern_ref(p * a, outer, inner); + verify_sparsity_pattern_ref_matches(permuted_pattern, permuted); + + RealSparseMatrixType real_part = a.real(); + const internal::SparsityPatternRef<int> real_pattern = internal::make_col_major_pattern_ref(a.real(), outer, inner); + verify_sparsity_pattern_ref_matches(real_pattern, real_part); + + RowMajorSparseMatrixType row_major = real_part; + const internal::SparsityPatternRef<int> row_major_pattern = + internal::make_col_major_pattern_ref(row_major, outer, inner); + verify_sparsity_pattern_ref_matches(row_major_pattern, real_part); +} + template <int> void bug1105() { // Regression test for bug 1105 @@ -1160,5 +1210,6 @@ CALL_SUBTEST_5(bug1105<0>()); CALL_SUBTEST_1(sparse_sub_assign_eigenbase<0>()); CALL_SUBTEST_1(ambivector_coeff<0>()); + CALL_SUBTEST_1(sparsity_pattern_ref_sparse_expressions<0>()); } #endif