Core: Add in-place self-adjoint permutations and share Cholesky pivots libeigen/eigen!3395 Closes #3174 Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>
diff --git a/CHANGELOG.md b/CHANGELOG.md index 2a6ef15..972e692 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md
@@ -9,6 +9,7 @@ - Support building Eigen's in-tree BLAS library with 64-bit (`ILP64`) integers (`eigen_blas_ilp64` and `eigen_blas_ilp64_static`, controlled by `EIGEN_BUILD_BLAS_ILP64`) - `PreconditionSquareMatrix` option, which applies JacobiSVD's QR preconditioners to square matrices. BDCSVD forwards this option to JacobiSVD for matrices smaller than the switch size [!3212] - `DenseBase::replicate` and `VectorwiseOp::replicate` accept `Eigen::fix<N>` factors [!3262, !3282] +- `SelfAdjointView::applyPermutation()` and `SelfAdjointView::applyTransposition()`, which replace the viewed matrix `S` by `P S P^T` in place, reading and writing only the stored triangle [#3174, !3395] Bug fixes: - Row-major triangular matrix-vector products could return incorrect results or double-destroy temporaries for scalars requiring initialization, and could leak some large right-hand-side buffers [!2818]
diff --git a/Eigen/src/Cholesky/BunchKaufman.h b/Eigen/src/Cholesky/BunchKaufman.h index d35554a..63462cb 100644 --- a/Eigen/src/Cholesky/BunchKaufman.h +++ b/Eigen/src/Cholesky/BunchKaufman.h
@@ -414,31 +414,6 @@ template <> struct bunch_kaufman<Lower> { - // Interchange rows and columns kk and kp (kp >= kk) in the lower triangle of the Hermitian matrix - // `mat`, including the already-computed factor columns to the left of column `kfirst` (so that the - // stored unit triangular factor stays consistent with a single up-front permutation P). `kfirst` is - // the first column of the pivot block (== kk for a 1x1 pivot, == kk-1 for a 2x2 pivot). - template <typename MatrixType> - static void apply_symmetric_pivot(MatrixType& mat, Index kfirst, Index kk, Index kp, Index kstep) { - using Scalar = typename MatrixType::Scalar; - const Index n = mat.rows(); - const Index s = n - kp - 1; - if (s > 0) mat.col(kk).tail(s).swap(mat.col(kp).tail(s)); - for (Index i = kk + 1; i < kp; ++i) { - Scalar tmp = mat.coeff(i, kk); - mat.coeffRef(i, kk) = numext::conj(mat.coeff(kp, i)); - mat.coeffRef(kp, i) = numext::conj(tmp); - } - numext::swap(mat.coeffRef(kk, kk), mat.coeffRef(kp, kp)); - EIGEN_IF_CONSTEXPR (NumTraits<Scalar>::IsComplex) { - mat.coeffRef(kp, kk) = numext::conj(mat.coeff(kp, kk)); - } - if (kfirst > 0) mat.row(kk).head(kfirst).swap(mat.row(kp).head(kfirst)); - if (kstep == 2) { - numext::swap(mat.coeffRef(kfirst + 1, kfirst), mat.coeffRef(kp, kfirst)); - } - } - // Unblocked (level-2 BLAS) Bunch-Kaufman factorization of the lower triangle of `mat`, in place. // Columns [k0, n) are factorized. On output: // - the strictly-lower triangle holds the unit lower factor L (with explicit zeros at the @@ -501,7 +476,7 @@ } const Index kk = k + kstep - 1; // column of the pivot block to interchange with kp - if (kp != kk) apply_symmetric_pivot(mat, k, kk, kp, kstep); + if (kp != kk) selfadjoint_swap_lower(mat, kk, kp); if (kstep == 1) { transpositions.coeffRef(k) = StorageIndex(kp); @@ -674,7 +649,7 @@ const Index kk = jc + kstep - 1; if (kp != kk) { - apply_symmetric_pivot(mat, jc, kk, kp, kstep); + selfadjoint_swap_lower(mat, kk, kp); const Index nwc = j + kstep; // number of populated W columns W.row(kk).head(nwc).swap(W.row(kp).head(nwc)); } @@ -1086,15 +1061,13 @@ template <typename MatrixType, int UpLo_> MatrixType BunchKaufman<MatrixType, UpLo_>::reconstructedMatrix() const { eigen_assert(m_isInitialized && "BunchKaufman is not initialized."); - const Index size = m_matrix.rows(); - MatrixType res(size, size); - - res.setIdentity(); - res = transpositionsP() * res; // P - res = matrixU() * res; // U P = L^* P - applyD(res); // D L^* P - res = matrixL() * res; // L D L^* P - res = transpositionsP().transpose() * res; // P^T L D L^* P + MatrixType res = matrixU(); + applyD(res); + res = matrixL() * res; + // P^T (L D L^*) P; the stored pivots cannot alias the new result. + auto lower_triangle = res.template selfadjointView<Lower>(); + for (Index k = m_matrix.rows() - 1; k >= 0; --k) lower_triangle.applyTransposition(k, m_transpositions.coeff(k)); + res.template triangularView<StrictlyUpper>() = res.adjoint(); return res; }
diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index deca15e..5594e61 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h
@@ -383,19 +383,7 @@ transpositions.coeffRef(k) = IndexType(index_of_biggest_in_corner); if (k != index_of_biggest_in_corner) { - // apply the transposition while taking care to consider only - // the lower triangular part - Index s = size - index_of_biggest_in_corner - 1; // trailing size after the biggest element - mat.row(k).head(k).swap(mat.row(index_of_biggest_in_corner).head(k)); - mat.col(k).tail(s).swap(mat.col(index_of_biggest_in_corner).tail(s)); - std::swap(mat.coeffRef(k, k), mat.coeffRef(index_of_biggest_in_corner, index_of_biggest_in_corner)); - for (Index i = k + 1; i < index_of_biggest_in_corner; ++i) { - Scalar tmp = mat.coeffRef(i, k); - mat.coeffRef(i, k) = numext::conj(mat.coeffRef(index_of_biggest_in_corner, i)); - mat.coeffRef(index_of_biggest_in_corner, i) = numext::conj(tmp); - } - EIGEN_IF_CONSTEXPR (NumTraits<Scalar>::IsComplex) - mat.coeffRef(index_of_biggest_in_corner, k) = numext::conj(mat.coeff(index_of_biggest_in_corner, k)); + selfadjoint_swap_lower(mat, k, index_of_biggest_in_corner); } // partition the matrix: @@ -725,20 +713,13 @@ template <typename MatrixType, int UpLo_> MatrixType LDLT<MatrixType, UpLo_>::reconstructedMatrix() const { eigen_assert(m_isInitialized && "LDLT is not initialized."); - const Index size = m_matrix.rows(); - MatrixType res(size, size); - - // P - res.setIdentity(); - res = transpositionsP() * res; - // L^* P - res = matrixU() * res; - // D(L^*P) + MatrixType res = matrixU(); res = vectorD().real().asDiagonal() * res; - // L(DL^*P) res = matrixL() * res; - // P^T (LDL^*P) - res = transpositionsP().transpose() * res; + // P^T (L D L^*) P; the stored pivots cannot alias the new result. + auto lower_triangle = res.template selfadjointView<Lower>(); + for (Index k = m_matrix.rows() - 1; k >= 0; --k) lower_triangle.applyTransposition(k, m_transpositions.coeff(k)); + res.template triangularView<StrictlyUpper>() = res.adjoint(); return res; }
diff --git a/Eigen/src/Core/PermutationMatrix.h b/Eigen/src/Core/PermutationMatrix.h index 212f8f6..9e95178 100644 --- a/Eigen/src/Core/PermutationMatrix.h +++ b/Eigen/src/Core/PermutationMatrix.h
@@ -226,25 +226,25 @@ * * \note \blank \note_try_to_help_rvo */ - inline InverseReturnType inverse() const { return InverseReturnType(derived()); } + EIGEN_DEVICE_FUNC inline InverseReturnType inverse() const { return InverseReturnType(derived()); } /** \returns the transpose permutation matrix. * * \note \blank \note_try_to_help_rvo */ - inline InverseReturnType transpose() const { return InverseReturnType(derived()); } + EIGEN_DEVICE_FUNC inline InverseReturnType transpose() const { return InverseReturnType(derived()); } /** \returns the adjoint of the permutation matrix. Its entries are real and it is orthogonal, so this equals * transpose() and inverse(). * * \note \blank \note_try_to_help_rvo */ - InverseReturnType adjoint() const { return InverseReturnType(derived()); } + EIGEN_DEVICE_FUNC InverseReturnType adjoint() const { return InverseReturnType(derived()); } /**** multiplication helpers to hopefully get RVO ****/ #ifndef EIGEN_PARSED_BY_DOXYGEN protected: template <typename OtherDerived> - void assignTranspose(const PermutationBase<OtherDerived>& other) { + EIGEN_DEVICE_FUNC void assignTranspose(const PermutationBase<OtherDerived>& other) { for (Index i = 0; i < rows(); ++i) indices().coeffRef(Index(other.indices().coeff(i))) = StorageIndex(i); } template <typename Lhs, typename Rhs> @@ -405,7 +405,7 @@ #ifndef EIGEN_PARSED_BY_DOXYGEN template <typename Other> - PermutationMatrix(const InverseImpl<Other, PermutationStorage>& other) + EIGEN_DEVICE_FUNC PermutationMatrix(const InverseImpl<Other, PermutationStorage>& other) : m_indices(other.derived().nestedExpression().size()) { internal::convert_index<StorageIndex>(m_indices.size()); // asserts that size fits StorageIndex Base::assignTranspose(other.derived().nestedExpression());
diff --git a/Eigen/src/Core/SelfAdjointView.h b/Eigen/src/Core/SelfAdjointView.h index 18850c3..21cffb3 100644 --- a/Eigen/src/Core/SelfAdjointView.h +++ b/Eigen/src/Core/SelfAdjointView.h
@@ -34,6 +34,42 @@ namespace internal { +// Requires 0 <= i < j < matrix.rows(). The row prefixes may contain already-computed factors. +template <typename MatrixType> +EIGEN_DEVICE_FUNC void selfadjoint_swap_lower(MatrixType& matrix, Index i, Index j) { + using Scalar = typename MatrixType::Scalar; + const Index size = matrix.rows(); + eigen_internal_assert(matrix.cols() == size && i >= 0 && i < j && j < size); + matrix.template block<1, Dynamic>(i, 0, 1, i).swap(matrix.template block<1, Dynamic>(j, 0, 1, i)); + for (Index k = i + 1; k < j; ++k) { + const Scalar value = matrix.coeff(k, i); + matrix.coeffRef(k, i) = numext::conj(matrix.coeff(j, k)); + matrix.coeffRef(j, k) = numext::conj(value); + } + matrix.template block<Dynamic, 1>(j + 1, i, size - j - 1, 1) + .swap(matrix.template block<Dynamic, 1>(j + 1, j, size - j - 1, 1)); + EIGEN_IF_CONSTEXPR (NumTraits<Scalar>::IsComplex) { + matrix.coeffRef(j, i) = numext::conj(matrix.coeff(j, i)); + } + numext::swap(matrix.coeffRef(i, i), matrix.coeffRef(j, j)); +} + +// Consumes an owning copy of the permutation, disjoint from matrix. +template <typename MatrixType, typename PermutationType> +EIGEN_DEVICE_FUNC void selfadjoint_apply_permutation_lower(MatrixType& matrix, PermutationType& remaining) { + const Index size = matrix.rows(); + for (Index i = 0; i < size; ++i) { + while (Index(remaining.indices().coeff(i)) != i) { + const Index j = convert_index<Index>(remaining.indices().coeff(i)); + // Earlier indices are fixed, so j > i. Each swap removes one edge from the permutation cycle. + // If P(j) == j as well, index j repeats and this loop would not terminate. + eigen_assert(i < j && j < size && Index(remaining.indices().coeff(j)) != j && "indices are not a permutation"); + selfadjoint_swap_lower(matrix, i, j); + numext::swap(remaining.indices().coeffRef(i), remaining.indices().coeffRef(j)); + } + } +} + template <typename MatrixType, unsigned int UpLo> struct traits<SelfAdjointView<MatrixType, UpLo>> : traits<MatrixType> { using MatrixTypeNested = typename ref_selector<MatrixType>::non_const_type; @@ -130,6 +166,69 @@ return *this; } + /** Replaces this view's matrix \f$ S \f$ by \f$ P S P^T \f$ in place, where \a permutation represents \f$ P \f$. + * Only the referenced triangle is read or written; the opposite triangle is unchanged. Complex coefficients + * crossing the diagonal are conjugated. The matrix must be square and the permutation must have matching size. + * The operation takes \f$ O(n^2) \f$ time and \f$ O(n) \f$ auxiliary storage for a copy of the permutation. + * + * \returns a reference to \c *this, retaining its self-adjoint structure. + * + * \code + * A.selfadjointView<Lower>().applyPermutation(P); + * \endcode + */ + template <typename PermutationDerived> + EIGEN_DEVICE_FUNC SelfAdjointView& applyPermutation(const PermutationBase<PermutationDerived>& permutation) { + typename PermutationDerived::PlainPermutationType remaining(permutation); + return applyPermutationInPlace(remaining); + } + + /** Applies an inverse (or transposed) permutation in place. \sa applyPermutation() */ + template <typename PermutationDerived> + EIGEN_DEVICE_FUNC SelfAdjointView& applyPermutation( + const InverseImpl<PermutationDerived, PermutationStorage>& permutation) { + typename PermutationDerived::PlainPermutationType remaining(permutation); + return applyPermutationInPlace(remaining); + } + + /** Applies the permutation represented by a sequence of transpositions in place. + * The indices are copied before any matrix writes, including when they alias the matrix. + * The operation takes \f$ O(n^2) \f$ time and \f$ O(n) \f$ auxiliary storage. + * \sa applyPermutation(), applyTransposition() + */ + template <typename TranspositionsDerived> + EIGEN_DEVICE_FUNC SelfAdjointView& applyPermutation(const TranspositionsBase<TranspositionsDerived>& transpositions) { + return applyTranspositions<false>(transpositions); + } + + /** Applies an inverse (or transposed) sequence of transpositions in place. \sa applyPermutation() */ + template <typename TranspositionsDerived> + EIGEN_DEVICE_FUNC SelfAdjointView& applyPermutation( + const Transpose<TranspositionsBase<TranspositionsDerived>>& transpositions) { + return applyTranspositions<true>(transpositions.nestedExpression()); + } + + /** Interchanges rows and columns \a i and \a j of this self-adjoint view in place. + * Only the referenced triangle is read or written; complex coefficients crossing the diagonal are conjugated. + * The matrix must be square and both indices must lie in [0, rows()). Equal indices leave the view unchanged. + * The operation takes \f$ O(n) \f$ time and \f$ O(1) \f$ auxiliary storage, without heap allocation. + * \returns a reference to \c *this. + * \sa applyPermutation() + */ + EIGEN_DEVICE_FUNC SelfAdjointView& applyTransposition(Index i, Index j) { + EIGEN_STATIC_ASSERT_LVALUE(MatrixType); + eigen_assert(this->rows() == this->cols() && i >= 0 && j >= 0 && i < this->rows() && j < this->rows()); + if (i == j) return *this; + if (i > j) numext::swap(i, j); + EIGEN_IF_CONSTEXPR (UpLo == Lower) { + internal::selfadjoint_swap_lower(m_matrix, i, j); + } else { + auto transposed = m_matrix.transpose(); + internal::selfadjoint_swap_lower(transposed, i, j); + } + return *this; + } + /** \internal */ EIGEN_DEVICE_FUNC constexpr const MatrixTypeNestedCleaned& _expression() const noexcept { return m_matrix; } @@ -234,6 +333,34 @@ EIGEN_DEVICE_FUNC EigenvaluesReturnType eigenvalues() const; EIGEN_DEVICE_FUNC RealScalar operatorNorm() const; + private: + template <typename PermutationType> + EIGEN_DEVICE_FUNC SelfAdjointView& applyPermutationInPlace(PermutationType& remaining) { + EIGEN_STATIC_ASSERT_LVALUE(MatrixType); + eigen_assert(this->rows() == this->cols() && this->rows() == remaining.size()); + EIGEN_IF_CONSTEXPR (UpLo == Lower) { + internal::selfadjoint_apply_permutation_lower(m_matrix, remaining); + } else { + auto transposed = m_matrix.transpose(); + internal::selfadjoint_apply_permutation_lower(transposed, remaining); + } + return *this; + } + + template <bool Reverse, typename TranspositionsDerived> + EIGEN_DEVICE_FUNC SelfAdjointView& applyTranspositions( + const TranspositionsBase<TranspositionsDerived>& transpositions) { + EIGEN_STATIC_ASSERT_LVALUE(MatrixType); + eigen_assert(this->rows() == this->cols() && this->rows() == transpositions.size()); + const typename TranspositionsDerived::IndicesType::PlainObject indices(transpositions.indices()); + const Index size = this->rows(); + for (Index k = 0; k < size; ++k) { + const Index i = Reverse ? size - 1 - k : k; + applyTransposition(i, internal::convert_index<Index>(indices.coeff(i))); + } + return *this; + } + protected: MatrixTypeNested m_matrix; };
diff --git a/benchmarks/Cholesky/CMakeLists.txt b/benchmarks/Cholesky/CMakeLists.txt index ec05343..134f305 100644 --- a/benchmarks/Cholesky/CMakeLists.txt +++ b/benchmarks/Cholesky/CMakeLists.txt
@@ -9,3 +9,5 @@ eigen_add_benchmark(bench_llt_inverse_cplx bench_llt_inverse.cpp DEFINITIONS "SCALAR=std::complex<double>") eigen_add_benchmark(bench_llt_inverse_float bench_llt_inverse.cpp DEFINITIONS SCALAR=float) eigen_add_benchmark(bench_llt_inverse_cplxf bench_llt_inverse.cpp DEFINITIONS "SCALAR=std::complex<float>") + +eigen_add_benchmark(bench_cholesky_reconstruction bench_cholesky_reconstruction.cpp)
diff --git a/benchmarks/Cholesky/bench_cholesky_reconstruction.cpp b/benchmarks/Cholesky/bench_cholesky_reconstruction.cpp new file mode 100644 index 0000000..138e9ae --- /dev/null +++ b/benchmarks/Cholesky/bench_cholesky_reconstruction.cpp
@@ -0,0 +1,54 @@ +// SPDX-FileCopyrightText: The Eigen Authors +// SPDX-License-Identifier: MPL-2.0 + +#include <benchmark/benchmark.h> +#include <Eigen/Cholesky> +#include <cstdlib> +#include <type_traits> + +namespace { +using namespace Eigen; + +template <typename Scalar, int Order, int UpLo, bool Indefinite> +void BM_CholeskyReconstruction(benchmark::State& state) { + using MatrixType = Matrix<Scalar, Dynamic, Dynamic, Order>; + using RealScalar = typename MatrixType::RealScalar; + using Factorization = + typename std::conditional<Indefinite, BunchKaufman<MatrixType, UpLo>, LDLT<MatrixType, UpLo>>::type; + const Index size = state.range(0); + std::srand(3395); + const MatrixType random = MatrixType::Random(size, size); + MatrixType input(size, size); + if (Indefinite) { + input = random + random.adjoint(); + } else { + input.noalias() = random * random.adjoint(); + input.diagonal().array() += Scalar(1); + } + const Factorization factorization(input); + const MatrixType reconstructed = factorization.reconstructedMatrix(); + const RealScalar tolerance = RealScalar(64) * RealScalar(size) * NumTraits<RealScalar>::epsilon() * input.norm(); + if (factorization.info() != Success || !reconstructed.allFinite() || !((reconstructed - input).norm() <= tolerance)) { + state.SkipWithError("Incorrect Cholesky reconstruction"); + return; + } + for (auto _ : state) { + const MatrixType result = factorization.reconstructedMatrix(); + benchmark::DoNotOptimize(result.data()); + benchmark::ClobberMemory(); + } +} + +// Indefinite=false uses LDLT; true uses BunchKaufman. Factorization is outside timing. +#define EIGEN_BENCH_RECONSTRUCTION(Scalar, Order, UpLo) \ + BENCHMARK_TEMPLATE(BM_CholeskyReconstruction, Scalar, Order, UpLo, false) \ + ->ArgsProduct({{4, 16, 64, 128, 256, 1024}}); \ + BENCHMARK_TEMPLATE(BM_CholeskyReconstruction, Scalar, Order, UpLo, true)->ArgsProduct({{4, 16, 64, 128, 256, 1024}}) + +EIGEN_BENCH_RECONSTRUCTION(double, ColMajor, Lower); +EIGEN_BENCH_RECONSTRUCTION(double, RowMajor, Upper); +EIGEN_BENCH_RECONSTRUCTION(std::complex<double>, ColMajor, Upper); +EIGEN_BENCH_RECONSTRUCTION(std::complex<double>, RowMajor, Lower); + +#undef EIGEN_BENCH_RECONSTRUCTION +} // namespace
diff --git a/benchmarks/Core/CMakeLists.txt b/benchmarks/Core/CMakeLists.txt index 173c797..b8d0edb 100644 --- a/benchmarks/Core/CMakeLists.txt +++ b/benchmarks/Core/CMakeLists.txt
@@ -79,6 +79,7 @@ eigen_add_benchmark(bench_stable_norm bench_stable_norm.cpp) eigen_add_benchmark(bench_approx_comparisons bench_approx_comparisons.cpp) eigen_add_benchmark(bench_selfadjoint_l1norm bench_selfadjoint_l1norm.cpp) +eigen_add_benchmark(bench_selfadjoint_permutation bench_selfadjoint_permutation.cpp) eigen_add_benchmark(bench_cwise_math bench_cwise_math.cpp) eigen_add_benchmark(bench_complex_sqrt bench_complex_sqrt.cpp) eigen_add_benchmark(bench_cast bench_cast.cpp)
diff --git a/benchmarks/Core/bench_selfadjoint_permutation.cpp b/benchmarks/Core/bench_selfadjoint_permutation.cpp new file mode 100644 index 0000000..292d42c --- /dev/null +++ b/benchmarks/Core/bench_selfadjoint_permutation.cpp
@@ -0,0 +1,83 @@ +// SPDX-FileCopyrightText: The Eigen Authors +// SPDX-License-Identifier: MPL-2.0 + +#include <benchmark/benchmark.h> +#include <Eigen/Core> +#include <algorithm> +#include <cstdlib> +#include <random> + +namespace { +using namespace Eigen; + +template <bool InPlace, unsigned int UpLo, typename MatrixType, typename PermutationType> +void permute(MatrixType& matrix, const PermutationType& permutation, bool inverse) { + EIGEN_IF_CONSTEXPR (InPlace) { + if (inverse) + matrix.template selfadjointView<UpLo>().applyPermutation(permutation.inverse()); + else + matrix.template selfadjointView<UpLo>().applyPermutation(permutation); + } else { + if (inverse) + matrix.template selfadjointView<UpLo>() = + (permutation.inverse() * matrix.template selfadjointView<UpLo>() * permutation).eval(); + else + matrix.template selfadjointView<UpLo>() = + (permutation * matrix.template selfadjointView<UpLo>() * permutation.transpose()).eval(); + } +} + +template <typename Scalar, int Order, unsigned int UpLo, bool InPlace> +void BM_SelfAdjointPermutation(benchmark::State& state) { + using MatrixType = Matrix<Scalar, Dynamic, Dynamic, Order>; + const Index size = state.range(0); + std::srand(42); + MatrixType matrix = MatrixType::Random(size, size); + matrix.diagonal() = matrix.diagonal().real().template cast<Scalar>(); + PermutationMatrix<Dynamic> permutation(size); + permutation.setIdentity(); + std::mt19937 random(42); + if (state.range(1) == 0) { + std::shuffle(permutation.indices().data(), permutation.indices().data() + size, random); + } else { + permutation.indices().reverseInPlace(); + } + const MatrixType full = matrix.template selfadjointView<UpLo>(); + const bool inverse = state.range(2) != 0; + MatrixType reference(size, size); + if (inverse) + reference = permutation.inverse() * full * permutation; + else + reference = permutation * full * permutation.transpose(); + MatrixType expected = matrix; + expected.template triangularView<UpLo>() = reference; + permute<InPlace, UpLo>(matrix, permutation, inverse); + if (matrix != expected) { + state.SkipWithError("Incorrect self-adjoint permutation or modified opposite triangle"); + return; + } + for (auto _ : state) { + permute<InPlace, UpLo>(matrix, permutation, inverse); + benchmark::DoNotOptimize(matrix.data()); + benchmark::ClobberMemory(); + } +} + +// Permutation kind: 0 = shuffled (seed 42), 1 = reversal. Last argument selects the inverse. +#define EIGEN_BENCH_SELFADJOINT_PERMUTATION(Scalar, Order, UpLo) \ + BENCHMARK_TEMPLATE(BM_SelfAdjointPermutation, Scalar, Order, UpLo, false) \ + ->ArgsProduct({{4, 16, 64, 256, 1024}, {0, 1}, {0, 1}}); \ + BENCHMARK_TEMPLATE(BM_SelfAdjointPermutation, Scalar, Order, UpLo, true) \ + ->ArgsProduct({{4, 16, 64, 256, 1024}, {0, 1}, {0, 1}}) + +EIGEN_BENCH_SELFADJOINT_PERMUTATION(float, ColMajor, Lower); +EIGEN_BENCH_SELFADJOINT_PERMUTATION(float, RowMajor, Upper); +EIGEN_BENCH_SELFADJOINT_PERMUTATION(double, ColMajor, Lower); +EIGEN_BENCH_SELFADJOINT_PERMUTATION(double, RowMajor, Upper); +EIGEN_BENCH_SELFADJOINT_PERMUTATION(std::complex<float>, ColMajor, Upper); +EIGEN_BENCH_SELFADJOINT_PERMUTATION(std::complex<float>, RowMajor, Lower); +EIGEN_BENCH_SELFADJOINT_PERMUTATION(std::complex<double>, ColMajor, Upper); +EIGEN_BENCH_SELFADJOINT_PERMUTATION(std::complex<double>, RowMajor, Lower); + +#undef EIGEN_BENCH_SELFADJOINT_PERMUTATION +} // namespace
diff --git a/doc/QuickReference.dox b/doc/QuickReference.dox index 91f006a..7f5ad2d 100644 --- a/doc/QuickReference.dox +++ b/doc/QuickReference.dox
@@ -732,6 +732,15 @@ m3 -= s1 * m3.adjoint() * m1.selfadjointView<Eigen::Lower>();\endcode </td></tr> <tr><td> +Symmetric permutation in place: \n +\f$ S := P S P^T \f$ (only the selected triangle is accessed) +</td><td>\code +M.selfadjointView<Eigen::Lower>().applyPermutation(P); +M.selfadjointView<Eigen::Upper>().applyPermutation(P.inverse()); +// P can also be a Transpositions sequence. +M.selfadjointView<Eigen::Lower>().applyTransposition(i,j);\endcode +</td></tr> +<tr><td> Rank 1 and rank K update: \n \f$ upper(M_1) \mathrel{{+}{=}} s_1 M_2 M_2^* \f$ \n \f$ lower(M_1) \mathbin{{-}{=}} M_2^* M_2 \f$
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 323dd1c..e4c5998 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt
@@ -442,6 +442,7 @@ ei_add_test(trsm_packet_cancellation) ei_add_test(triangular_in_place) ei_add_test(selfadjoint) +ei_add_test(selfadjoint_permutation) ei_add_test(product_selfadjoint) ei_add_test(product_symm) ei_add_test(product_syrk) @@ -458,6 +459,7 @@ ei_add_test(bandmatrix) ei_add_test(cholesky) ei_add_test(bunchkaufman) +ei_add_test(cholesky_reconstruction) ei_add_test(condition_estimator) ei_add_test(lu) ei_add_test(determinant)
diff --git a/test/cholesky_reconstruction.cpp b/test/cholesky_reconstruction.cpp new file mode 100644 index 0000000..4a2f769 --- /dev/null +++ b/test/cholesky_reconstruction.cpp
@@ -0,0 +1,78 @@ +// SPDX-FileCopyrightText: The Eigen Authors +// SPDX-License-Identifier: MPL-2.0 + +#define EIGEN_RUNTIME_NO_MALLOC +#include "main.h" +#include <Eigen/Cholesky> + +template <typename MatrixType, typename Factorization> +void check_cholesky_reconstruction(const MatrixType& input, const Factorization& factorization) { + using RealScalar = typename MatrixType::RealScalar; + VERIFY_IS_EQUAL(factorization.info(), Success); + const MatrixType reconstructed = factorization.reconstructedMatrix(); + // Two triangular products and the factorization each accumulate O(n) rounding errors per coefficient. + const RealScalar tolerance = + RealScalar(64) * RealScalar((std::max)(Index(1), input.rows())) * NumTraits<RealScalar>::epsilon() * input.norm(); + VERIFY(reconstructed.allFinite()); + VERIFY((reconstructed - input).norm() <= tolerance); + VERIFY_IS_EQUAL(MatrixType(reconstructed.template triangularView<StrictlyUpper>()), + MatrixType(reconstructed.adjoint().template triangularView<StrictlyUpper>())); +} + +template <typename Scalar, int Order, int UpLo> +void cholesky_reconstruction() { + using MatrixType = Matrix<Scalar, Dynamic, Dynamic, Order>; + for (Index size : {0, 1, 2, 5, 17, 63, 64, 65}) { + const MatrixType random = MatrixType::Random(size, size); + MatrixType input = random * random.adjoint(); + input.diagonal().array() += Scalar(1); + check_cholesky_reconstruction(input, LDLT<MatrixType, UpLo>(input)); + input = random + random.adjoint(); + check_cholesky_reconstruction(input, BunchKaufman<MatrixType, UpLo>(input)); + } +} + +template <typename Scalar, int Order, int UpLo> +void cholesky_reconstruction_fixed() { + using MatrixType = Matrix<Scalar, 4, 4, Order>; + MatrixType input = MatrixType::Zero(); + input.diagonal() << Scalar(1), Scalar(4), Scalar(2), Scalar(3); + LDLT<MatrixType, UpLo> ldlt(input); + VERIFY_IS_EQUAL(ldlt.transpositionsP().coeff(0), 1); + internal::set_is_malloc_allowed(false); + check_cholesky_reconstruction(input, ldlt); + internal::set_is_malloc_allowed(true); + + // The off-diagonal entry forces a 2x2 pivot and a non-adjacent symmetric interchange. + const Scalar offDiagonal = Scalar(2) + Scalar(0.5) * internal::random<Scalar>(); + input.setZero(); + input(1, 1) = Scalar(2); + input(2, 2) = Scalar(-3); + input(3, 0) = offDiagonal; + input(0, 3) = numext::conj(offDiagonal); + BunchKaufman<MatrixType, UpLo> bk(input); + VERIFY_IS_EQUAL(bk.transpositionsP().coeff(1), 3); + VERIFY(!numext::is_exactly_zero(bk.subDiagonal().coeff(0))); + internal::set_is_malloc_allowed(false); + check_cholesky_reconstruction(input, bk); + internal::set_is_malloc_allowed(true); +} + +EIGEN_DECLARE_TEST(cholesky_reconstruction) { + CALL_SUBTEST_1((cholesky_reconstruction<double, ColMajor, Lower>())); + CALL_SUBTEST_1((cholesky_reconstruction<double, ColMajor, Upper>())); + CALL_SUBTEST_2((cholesky_reconstruction<double, RowMajor, Lower>())); + CALL_SUBTEST_2((cholesky_reconstruction<double, RowMajor, Upper>())); + CALL_SUBTEST_3((cholesky_reconstruction<std::complex<double>, ColMajor, Lower>())); + CALL_SUBTEST_3((cholesky_reconstruction<std::complex<double>, ColMajor, Upper>())); + CALL_SUBTEST_4((cholesky_reconstruction<std::complex<double>, RowMajor, Lower>())); + CALL_SUBTEST_4((cholesky_reconstruction<std::complex<double>, RowMajor, Upper>())); + CALL_SUBTEST_5((cholesky_reconstruction<float, ColMajor, Lower>())); + CALL_SUBTEST_5((cholesky_reconstruction<std::complex<float>, RowMajor, Upper>())); + for (int i = 0; i < g_repeat; ++i) { + CALL_SUBTEST_1((cholesky_reconstruction_fixed<double, ColMajor, Lower>())); + CALL_SUBTEST_2((cholesky_reconstruction_fixed<double, RowMajor, Upper>())); + CALL_SUBTEST_3((cholesky_reconstruction_fixed<std::complex<double>, ColMajor, Upper>())); + CALL_SUBTEST_4((cholesky_reconstruction_fixed<std::complex<double>, RowMajor, Lower>())); + } +}
diff --git a/test/gpu_basic.cu b/test/gpu_basic.cu index a32c66e..889e3cd 100644 --- a/test/gpu_basic.cu +++ b/test/gpu_basic.cu
@@ -34,6 +34,23 @@ // } // }; +template <typename MatrixType, int UpLo> +struct selfadjoint_permutation { + EIGEN_DEVICE_FUNC void operator()(int i, const typename MatrixType::Scalar* in, + typename MatrixType::Scalar* out) const { + using namespace Eigen; + MatrixType matrix = Map<const MatrixType>(in + i); + PermutationMatrix<4> permutation; + permutation.indices() << 1, 3, 0, 2; + matrix.template selfadjointView<UpLo>().applyPermutation(permutation); + Map<MatrixType>(out + 48 * i) = matrix; + matrix.template selfadjointView<UpLo>().applyPermutation(permutation.inverse()); + Map<MatrixType>(out + 48 * i + 16) = matrix; + matrix.template selfadjointView<UpLo>().applyTransposition(3, 0); + Map<MatrixType>(out + 48 * i + 32) = matrix; + } +}; + template <typename Scalar> struct scaled_permutation_product { EIGEN_DEVICE_FUNC void operator()(int i, const float* in, float* out) const { @@ -1045,6 +1062,13 @@ CALL_SUBTEST(run_and_compare_to_gpu(selfadjoint_rank2_update<Matrix6f, Lower>(), nthreads, in, out)); CALL_SUBTEST(run_and_compare_to_gpu(selfadjoint_rank2_update<Matrix6f, Upper>(), nthreads, in, out)); + CALL_SUBTEST((run_and_compare_to_gpu(selfadjoint_permutation<Matrix4f, Lower>(), nthreads, in, out))); + CALL_SUBTEST( + (run_and_compare_to_gpu(selfadjoint_permutation<Matrix<float, 4, 4, RowMajor>, Upper>(), nthreads, in, out))); + CALL_SUBTEST((run_and_compare_to_gpu(selfadjoint_permutation<Matrix4cf, Upper>(), nthreads, cfin, cfout))); + CALL_SUBTEST((run_and_compare_to_gpu(selfadjoint_permutation<Matrix<std::complex<float>, 4, 4, RowMajor>, Lower>(), + nthreads, cfin, cfout))); + CALL_SUBTEST(run_and_compare_to_gpu(selfadjoint_l1_norm<Matrix4f, Lower>(), nthreads, in, out)); CALL_SUBTEST(run_and_compare_to_gpu(selfadjoint_l1_norm<Matrix4f, Upper>(), nthreads, in, out)); CALL_SUBTEST(run_and_compare_to_gpu(selfadjoint_l1_norm<Matrix6f, Lower>(), nthreads, in, out));
diff --git a/test/selfadjoint_permutation.cpp b/test/selfadjoint_permutation.cpp new file mode 100644 index 0000000..bfaf885 --- /dev/null +++ b/test/selfadjoint_permutation.cpp
@@ -0,0 +1,248 @@ +// SPDX-FileCopyrightText: The Eigen Authors +// SPDX-License-Identifier: MPL-2.0 + +#define TEST_ENABLE_TEMPORARY_TRACKING +#define EIGEN_RUNTIME_NO_MALLOC +#include "main.h" +#include <Eigen/Core> +#include <algorithm> +#include <cstdint> + +template <unsigned int UpLo, typename MatrixType, typename Transformation> +void check_selfadjoint_permutation(MatrixType& matrix, const Transformation& transformation) { + using IndicesType = typename Transformation::IndicesType; + using Permutation = PermutationMatrix<IndicesType::SizeAtCompileTime, IndicesType::MaxSizeAtCompileTime, + typename IndicesType::Scalar>; + // NOLINTNEXTLINE(performance-unnecessary-copy-initialization): Transpositions inputs need this conversion. + const Permutation permutation(transformation); + using PlainObject = typename MatrixType::PlainObject; + const PlainObject original = matrix; + const Index size = matrix.rows(); + PlainObject full(size, size); + for (Index j = 0; j < size; ++j) { + for (Index i = 0; i < size; ++i) { + full(Index(permutation.indices()[i]), Index(permutation.indices()[j])) = + (UpLo == Lower ? i >= j : i <= j) ? original(i, j) : numext::conj(original(j, i)); + } + } + PlainObject expected = original; + expected.template triangularView<UpLo>() = full; + auto view = matrix.template selfadjointView<UpLo>(); + VERIFY(&view.applyPermutation(transformation) == &view); + VERIFY_IS_EQUAL(matrix, expected); + VERIFY_IS_EQUAL(view.toDenseMatrix(), full); + view.applyPermutation(transformation.inverse()); + VERIFY_IS_EQUAL(matrix, original); + view.applyPermutation(transformation.transpose()).applyPermutation(transformation); + VERIFY_IS_EQUAL(matrix, original); +} + +template <unsigned int UpLo, typename MatrixType> +void check_selfadjoint_transposition(MatrixType& matrix, Index i, Index j) { + using PlainObject = typename MatrixType::PlainObject; + const PlainObject original = matrix; + PlainObject full = original.template selfadjointView<UpLo>(); + full.row(i).swap(full.row(j)); + full.col(i).swap(full.col(j)); + PlainObject expected = original; + expected.template triangularView<UpLo>() = full; + auto view = matrix.template selfadjointView<UpLo>(); + internal::set_is_malloc_allowed(false); + VERIFY(&view.applyTransposition(i, j) == &view); + internal::set_is_malloc_allowed(true); + VERIFY_IS_EQUAL(matrix, expected); + view.applyTransposition(j, i); + VERIFY_IS_EQUAL(matrix, original); +} + +template <typename Scalar, int Order> +void selfadjoint_permutation_cycles() { + using MatrixType = Matrix<Scalar, Dynamic, Dynamic, Order>; + for (Index size = 0; size <= 5; ++size) { + MatrixType matrix = MatrixType::Random(size, size); + matrix.diagonal() = matrix.diagonal().real().template cast<Scalar>(); + PermutationMatrix<Dynamic> permutation(size); + permutation.setIdentity(); + do { + check_selfadjoint_permutation<Lower>(matrix, permutation); + check_selfadjoint_permutation<Upper>(matrix, permutation); + } while (size > 1 && std::next_permutation(permutation.indices().data(), permutation.indices().data() + size)); + for (Index i = 0; i < size; ++i) { + for (Index j = 0; j < size; ++j) { + check_selfadjoint_transposition<Lower>(matrix, i, j); + check_selfadjoint_transposition<Upper>(matrix, i, j); + } + } + } +} + +template <typename Scalar, int Order> +void selfadjoint_permutation_views() { + using MatrixType = Matrix<Scalar, Dynamic, Dynamic, Order>; + const Index size = internal::random<Index>(2, 30); + MatrixType storage = MatrixType::Random(2 * size + 3, 2 * size + 3); + const MatrixType original = storage; + auto block = storage.block(1, 2, size, size); + PermutationMatrix<Dynamic, Dynamic, std::int64_t> permutation(size); + randomPermutationVector(permutation.indices(), size); + const Map<PermutationMatrix<Dynamic, Dynamic, std::int64_t>> mapped_permutation(permutation.indices().data(), size); + check_selfadjoint_permutation<Lower>(block, mapped_permutation); + check_selfadjoint_permutation<Upper>(block, permutation.indices().asPermutation()); + VERIFY_IS_EQUAL(storage, original); + + Map<MatrixType, Unaligned, Stride<Dynamic, Dynamic>> strided(storage.data(), size, size, + Stride<Dynamic, Dynamic>(2 * storage.outerStride(), 2)); + check_selfadjoint_permutation<Lower>(strided, permutation); + check_selfadjoint_permutation<Upper>(strided, permutation); + for (Index i = 0; i < size; ++i) { + check_selfadjoint_transposition<Lower>(block, i, size - 1 - i); + check_selfadjoint_transposition<Upper>(strided, i, size - 1 - i); + } + VERIFY_IS_EQUAL(storage, original); +} + +template <typename Scalar, int Size> +void selfadjoint_permutation_fixed() { + Matrix<Scalar, Size, Size> matrix = Matrix<Scalar, Size, Size>::Random(); + PermutationMatrix<Size> permutation; + randomPermutationVector(permutation.indices(), Size); + internal::set_is_malloc_allowed(false); + check_selfadjoint_permutation<Lower>(matrix, permutation); + check_selfadjoint_permutation<Upper>(matrix, permutation); + Transpositions<Size> transpositions; + for (Index i = 0; i < Size; ++i) transpositions.indices()[i] = int(Size - 1 - i); + check_selfadjoint_permutation<Lower>(matrix, transpositions); + check_selfadjoint_permutation<Upper>(matrix, transpositions); + internal::set_is_malloc_allowed(true); +} + +void selfadjoint_permutation_preconditions() { + MatrixXd matrix = MatrixXd::Random(3, 3); + PermutationMatrix<Dynamic> permutation(2); + permutation.setIdentity(); + VERIFY_RAISES_ASSERT(matrix.selfadjointView<Lower>().applyPermutation(permutation)); + VERIFY_RAISES_ASSERT(matrix.selfadjointView<Upper>().applyPermutation(permutation.inverse())); + VERIFY_RAISES_ASSERT(matrix.selfadjointView<Lower>().applyTransposition(-1, 0)); + VERIFY_RAISES_ASSERT(matrix.selfadjointView<Upper>().applyTransposition(0, 3)); + Transpositions<Dynamic> transpositions(2); + transpositions.setIdentity(); + VERIFY_RAISES_ASSERT(matrix.selfadjointView<Lower>().applyPermutation(transpositions)); + VERIFY_RAISES_ASSERT(matrix.selfadjointView<Upper>().applyPermutation(transpositions.inverse())); + // Indices that do not form a permutation: a repeat would otherwise loop forever. + PermutationMatrix<Dynamic> invalid(3); + invalid.indices() << 1, 1, 2; + VERIFY_RAISES_ASSERT(matrix.selfadjointView<Lower>().applyPermutation(invalid)); + VERIFY_RAISES_ASSERT(matrix.selfadjointView<Upper>().applyPermutation(invalid)); + invalid.indices() << 2, 0, 0; + VERIFY_RAISES_ASSERT(matrix.selfadjointView<Lower>().applyPermutation(invalid)); + invalid.indices() << 3, 0, 1; + VERIFY_RAISES_ASSERT(matrix.selfadjointView<Upper>().applyPermutation(invalid)); + matrix.resize(2, 3); + VERIFY_RAISES_ASSERT(matrix.selfadjointView<Lower>().applyPermutation(permutation)); + VERIFY_RAISES_ASSERT(matrix.selfadjointView<Upper>().applyPermutation(permutation)); + VERIFY_RAISES_ASSERT(matrix.selfadjointView<Lower>().applyTransposition(0, 1)); + VERIFY_RAISES_ASSERT(matrix.selfadjointView<Upper>().applyPermutation(transpositions)); + VERIFY_RAISES_ASSERT(matrix.selfadjointView<Lower>().applyPermutation(transpositions.transpose())); +} + +void selfadjoint_permutation_aliased_indices() { + Matrix3i matrix; + matrix << 1, 2, 0, 4, 5, 6, 7, 8, 9; + const PermutationMatrix<3> permutation(matrix.row(0).transpose()); + Matrix3i expected = matrix; + expected.selfadjointView<Upper>().applyPermutation(permutation); + matrix.selfadjointView<Upper>().applyPermutation(matrix.row(0).transpose().asPermutation()); + VERIFY_IS_EQUAL(matrix, expected); + matrix << 1, 2, 0, 4, 5, 6, 7, 8, 9; + expected = matrix; + expected.selfadjointView<Upper>().applyPermutation(permutation.inverse()); + matrix.selfadjointView<Upper>().applyPermutation(matrix.row(0).transpose().asPermutation().inverse()); + VERIFY_IS_EQUAL(matrix, expected); + + for (bool inverse : {false, true}) { + matrix << 1, 2, 0, 4, 5, 6, 7, 8, 9; + auto indices = matrix.row(0).transpose(); + TranspositionsWrapper<decltype(indices)> wrapped(indices); + const Transpositions<3> transpositions(wrapped); + expected = matrix; + if (inverse) { + expected.selfadjointView<Upper>().applyPermutation(transpositions.inverse()); + matrix.selfadjointView<Upper>().applyPermutation(wrapped.inverse()); + } else { + expected.selfadjointView<Upper>().applyPermutation(transpositions); + matrix.selfadjointView<Upper>().applyPermutation(wrapped); + } + VERIFY_IS_EQUAL(matrix, expected); + } +} + +template <typename Scalar, int Order> +void selfadjoint_transpositions() { + using MatrixType = Matrix<Scalar, Dynamic, Dynamic, Order>; + for (Index size = 0; size <= 4; ++size) { + MatrixType matrix = MatrixType::Random(size, size); + matrix.diagonal() = matrix.diagonal().real().template cast<Scalar>(); + Transpositions<Dynamic, Dynamic, std::int64_t> transpositions(size); + transpositions.indices().setZero(); + // All n^n sequences, including swaps with earlier indices and repeated targets. + bool finished = false; + do { + check_selfadjoint_permutation<Lower>(matrix, transpositions); + check_selfadjoint_permutation<Upper>(matrix, transpositions); + Index i = 0; + while (i < size && ++transpositions.indices()[i] == size) transpositions.indices()[i++] = 0; + finished = i == size; + } while (!finished); + transpositions.setIdentity(); + const Map<Transpositions<Dynamic, Dynamic, std::int64_t>> mapped(transpositions.indices().data(), size); + check_selfadjoint_permutation<Lower>(matrix, mapped); + check_selfadjoint_permutation<Upper>(matrix, mapped); + } +} + +void selfadjoint_permutation_temporary_count() { + MatrixXd matrix = MatrixXd::Random(8, 8); + PermutationMatrix<Dynamic> permutation(8); + randomPermutationVector(permutation.indices(), 8); + VERIFY_EVALUATION_COUNT(matrix.selfadjointView<Lower>().applyPermutation(permutation), 1); + VERIFY_EVALUATION_COUNT(matrix.selfadjointView<Upper>().applyPermutation(permutation.inverse()), 1); + VERIFY_EVALUATION_COUNT(matrix.selfadjointView<Lower>().applyPermutation(permutation.transpose()), 1); + Transpositions<Dynamic> transpositions(8); + transpositions.setIdentity(); + VERIFY_EVALUATION_COUNT(matrix.selfadjointView<Lower>().applyPermutation(transpositions), 1); + VERIFY_EVALUATION_COUNT(matrix.selfadjointView<Upper>().applyPermutation(transpositions.transpose()), 1); + VERIFY_EVALUATION_COUNT(matrix.selfadjointView<Lower>().applyTransposition(1, 7), 0); +} + +EIGEN_DECLARE_TEST(selfadjoint_permutation) { + CALL_SUBTEST_1((selfadjoint_permutation_cycles<float, ColMajor>())); + CALL_SUBTEST_1((selfadjoint_permutation_cycles<float, RowMajor>())); + CALL_SUBTEST_2((selfadjoint_permutation_cycles<double, ColMajor>())); + CALL_SUBTEST_2((selfadjoint_permutation_cycles<double, RowMajor>())); + CALL_SUBTEST_3((selfadjoint_permutation_cycles<std::complex<float>, ColMajor>())); + CALL_SUBTEST_3((selfadjoint_permutation_cycles<std::complex<float>, RowMajor>())); + CALL_SUBTEST_4((selfadjoint_permutation_cycles<std::complex<double>, ColMajor>())); + CALL_SUBTEST_4((selfadjoint_permutation_cycles<std::complex<double>, RowMajor>())); + CALL_SUBTEST_5((selfadjoint_permutation_cycles<int, ColMajor>())); + CALL_SUBTEST_5((selfadjoint_permutation_cycles<short, RowMajor>())); + CALL_SUBTEST_5((selfadjoint_permutation_cycles<half, ColMajor>())); + CALL_SUBTEST_5((selfadjoint_permutation_cycles<bfloat16, RowMajor>())); + for (int i = 0; i < g_repeat; ++i) { + CALL_SUBTEST_1((selfadjoint_permutation_views<float, ColMajor>())); + CALL_SUBTEST_2((selfadjoint_permutation_views<double, RowMajor>())); + CALL_SUBTEST_3((selfadjoint_permutation_views<std::complex<float>, RowMajor>())); + CALL_SUBTEST_4((selfadjoint_permutation_views<std::complex<double>, ColMajor>())); + CALL_SUBTEST_1((selfadjoint_permutation_fixed<float, 0>())); + CALL_SUBTEST_1((selfadjoint_permutation_fixed<float, 1>())); + CALL_SUBTEST_2((selfadjoint_permutation_fixed<double, 4>())); + CALL_SUBTEST_3((selfadjoint_permutation_fixed<std::complex<float>, 5>())); + } + CALL_SUBTEST_6((selfadjoint_transpositions<double, ColMajor>())); + CALL_SUBTEST_6((selfadjoint_transpositions<double, RowMajor>())); + CALL_SUBTEST_6((selfadjoint_transpositions<std::complex<double>, ColMajor>())); + CALL_SUBTEST_6((selfadjoint_transpositions<std::complex<double>, RowMajor>())); + CALL_SUBTEST_6((selfadjoint_transpositions<int, RowMajor>())); + CALL_SUBTEST_6(selfadjoint_permutation_temporary_count()); + CALL_SUBTEST_1(selfadjoint_permutation_preconditions()); + CALL_SUBTEST_5(selfadjoint_permutation_aliased_indices()); +}