Decompositions: Reduce cache thrashing in Schur and LU libeigen/eigen!3051 Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>
diff --git a/Eigen/src/Householder/Householder.h b/Eigen/src/Householder/Householder.h index 5271643..f977a29 100644 --- a/Eigen/src/Householder/Householder.h +++ b/Eigen/src/Householder/Householder.h
@@ -236,6 +236,47 @@ tau = conj((beta - c0) / beta); } +namespace internal { + +template <typename Derived, typename EssentialPart, + bool Fused = !Derived::IsRowMajor && EssentialPart::ColsAtCompileTime == 1 && + (EssentialPart::RowsAtCompileTime == 1 || EssentialPart::RowsAtCompileTime == 2)> +struct householder_apply_left_impl { + using Scalar = typename Derived::Scalar; + static EIGEN_DEVICE_FUNC void run(MatrixBase<Derived>& mat, const EssentialPart& essential, const Scalar& tau, + Scalar* workspace) { + Map<typename plain_row_type<typename Derived::PlainObject>::type> tmp(workspace, mat.cols()); + Block<Derived, EssentialPart::SizeAtCompileTime, Derived::ColsAtCompileTime> bottom(mat.derived(), 1, 0, + mat.rows() - 1, mat.cols()); + tmp.noalias() = essential.adjoint() * bottom.unwind(); + tmp = tau * (tmp + mat.row(0)); + mat.row(0) -= tmp; + bottom.unwind().noalias() -= essential * tmp; + } +}; + +// Two- and three-element reflectors on column-major storage: finishing each column before advancing replaces +// four strided row passes, which thrash the cache at large outer strides (issue #3160). +template <typename Derived, typename EssentialPart> +struct householder_apply_left_impl<Derived, EssentialPart, true> { + using Scalar = typename Derived::Scalar; + static EIGEN_DEVICE_FUNC void run(MatrixBase<Derived>& mat, const EssentialPart& essential, const Scalar& tau, + Scalar*) { + // Evaluated once so that the column loop reads plain coefficients; tau may reference a coefficient of mat. + const Matrix<Scalar, EssentialPart::RowsAtCompileTime, 1> v = essential; + const Scalar tauValue = tau; + Block<Derived, EssentialPart::RowsAtCompileTime, Derived::ColsAtCompileTime> bottom(mat.derived(), 1, 0, + mat.rows() - 1, mat.cols()); + for (Index j = 0; j < mat.cols(); ++j) { + const Scalar tmp = tauValue * (v.dot(bottom.col(j)) + mat.coeff(0, j)); + mat.coeffRef(0, j) -= tmp; + bottom.col(j) -= v * tmp; + } + } +}; + +} // namespace internal + /** Apply the elementary reflector H given by * \f$ H = I - tau v v^*\f$ * with @@ -258,13 +299,7 @@ if (rows() == 1) { *this *= Scalar(1) - tau; } else if (!numext::is_exactly_zero(tau)) { - Map<typename internal::plain_row_type<PlainObject>::type> tmp(workspace, cols()); - Block<Derived, EssentialPart::SizeAtCompileTime, Derived::ColsAtCompileTime> bottom(derived(), 1, 0, rows() - 1, - cols()); - tmp.noalias() = essential.adjoint() * bottom.unwind(); - tmp = tau * (tmp + this->row(0)); - this->row(0) = this->row(0) - tmp; - bottom.unwind().noalias() -= essential * tmp; + internal::householder_apply_left_impl<Derived, EssentialPart>::run(*this, essential, tau, workspace); } }
diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h index ff92663..06b5978 100644 --- a/Eigen/src/LU/PartialPivLU.h +++ b/Eigen/src/LU/PartialPivLU.h
@@ -363,6 +363,17 @@ using BlockType = Ref<Matrix<Scalar, Dynamic, Dynamic, StorageOrder>>; using RealScalar = typename MatrixType::RealScalar; + static void apply_row_transpositions(BlockType& matrix, Index first, Index count, const PivIndex* transpositions) { + EIGEN_IF_CONSTEXPR (StorageOrder == ColMajor) { + // Keep the pivot rows of one column in cache, even when the outer stride maps every column to the same set. + for (Index j = 0; j < matrix.cols(); ++j) + for (Index i = first; i < first + count; ++i) + numext::swap(matrix.coeffRef(i, j), matrix.coeffRef(transpositions[i], j)); + } else { + for (Index i = first; i < first + count; ++i) matrix.row(i).swap(matrix.row(transpositions[i])); + } + } + /** \internal performs the LU decomposition in-place of the matrix \a lu * using an unblocked algorithm. * @@ -489,15 +500,13 @@ // update permutations and apply them to A_0 if (k > 0) { BlockType A_0 = lu.block(0, 0, rows, k); - for (Index i = k; i < k + bs; ++i) { - Index piv = (row_transpositions[i] += internal::convert_index<PivIndex>(k)); - A_0.row(i).swap(A_0.row(piv)); - } + for (Index i = k; i < k + bs; ++i) row_transpositions[i] += internal::convert_index<PivIndex>(k); + apply_row_transpositions(A_0, k, bs, row_transpositions); } if (trows) { // apply permutations to A_2 - for (Index i = k; i < k + bs; ++i) A_2.row(i).swap(A_2.row(row_transpositions[i])); + apply_row_transpositions(A_2, k, bs, row_transpositions); // A12 = A11^-1 A12 A11.template triangularView<UnitLower>().solveInPlace(A12);
diff --git a/benchmarks/Eigenvalues/CMakeLists.txt b/benchmarks/Eigenvalues/CMakeLists.txt index 986d701..f59b368 100644 --- a/benchmarks/Eigenvalues/CMakeLists.txt +++ b/benchmarks/Eigenvalues/CMakeLists.txt
@@ -3,6 +3,7 @@ eigen_add_benchmark(bench_eigensolver bench_eigensolver.cpp) eigen_add_benchmark(bench_eigensolver_double bench_eigensolver.cpp DEFINITIONS SCALAR=double) +eigen_add_benchmark(bench_eigensolver_stride bench_eigensolver_stride.cpp) eigen_add_benchmark(bench_eig33 bench_eig33.cpp) eigen_add_benchmark(bench_tridiagonal_bisection bench_tridiagonal_bisection.cpp) eigen_add_benchmark(bench_tridiagonal_inverse_iteration bench_tridiagonal_inverse_iteration.cpp)
diff --git a/benchmarks/Eigenvalues/bench_eigensolver_stride.cpp b/benchmarks/Eigenvalues/bench_eigensolver_stride.cpp new file mode 100644 index 0000000..a73f8df --- /dev/null +++ b/benchmarks/Eigenvalues/bench_eigensolver_stride.cpp
@@ -0,0 +1,38 @@ +// SPDX-FileCopyrightText: The Eigen Authors +// SPDX-License-Identifier: MPL-2.0 + +#include <benchmark/benchmark.h> +#include <Eigen/Eigenvalues> +#include <cstdlib> + +using namespace Eigen; + +template <typename Scalar> +static void BM_EigenSolverStride(benchmark::State& state) { + using Mat = Matrix<Scalar, Dynamic, Dynamic>; + const Index n = state.range(0); + const bool computeVectors = state.range(1) != 0; + std::srand(1); + const Mat a = Mat::Random(n, n); + EigenSolver<Mat> solver(a, true); + if (solver.info() != Success) { + state.SkipWithError("EigenSolver did not converge"); + return; + } + const Mat vectors = solver.pseudoEigenvectors(); + const Mat values = solver.pseudoEigenvalueMatrix(); + const Scalar residual = (a * vectors - vectors * values).norm(); + const Scalar bound = Scalar(64 * n) * NumTraits<Scalar>::epsilon() * a.norm() * vectors.norm(); + if (!(residual <= bound)) { + state.SkipWithError("Eigenpair residual failed"); + return; + } + for (auto _ : state) { + solver.compute(a, computeVectors); + benchmark::DoNotOptimize(solver.eigenvalues().data()); + benchmark::ClobberMemory(); + } +} + +BENCHMARK_TEMPLATE(BM_EigenSolverStride, float)->ArgsProduct({{32, 128, 500, 512, 768, 1000, 1001, 1024}, {0, 1}}); +BENCHMARK_TEMPLATE(BM_EigenSolverStride, double)->ArgsProduct({{32, 128, 500, 512, 768, 1000, 1001, 1024}, {0, 1}});
diff --git a/benchmarks/Householder/CMakeLists.txt b/benchmarks/Householder/CMakeLists.txt index 8310c92..bb3f283 100644 --- a/benchmarks/Householder/CMakeLists.txt +++ b/benchmarks/Householder/CMakeLists.txt
@@ -2,3 +2,4 @@ # SPDX-License-Identifier: MPL-2.0 eigen_add_benchmark(bench_householder bench_householder.cpp) +eigen_add_benchmark(bench_householder_short bench_householder_short.cpp)
diff --git a/benchmarks/Householder/bench_householder_short.cpp b/benchmarks/Householder/bench_householder_short.cpp new file mode 100644 index 0000000..37bef4c --- /dev/null +++ b/benchmarks/Householder/bench_householder_short.cpp
@@ -0,0 +1,53 @@ +// SPDX-FileCopyrightText: The Eigen Authors +// SPDX-License-Identifier: MPL-2.0 + +#include <benchmark/benchmark.h> +#include <Eigen/Householder> +#include <cstdlib> + +using namespace Eigen; + +template <typename Scalar, int Size, int Side> +static void BM_ShortHouseholder(benchmark::State& state) { + using Mat = Matrix<Scalar, Dynamic, Dynamic>; + const Index n = state.range(0); + const Index stride = n + state.range(1); + std::srand(1); + Mat storage = Mat::Random(stride, n); + Matrix<Scalar, Dynamic, 1> workspace(n); + Matrix<Scalar, Size, 1> vector = Matrix<Scalar, Size, 1>::Random(); + Matrix<Scalar, Size - 1, 1> essential; + Scalar tau, beta; + vector.makeHouseholder(essential, tau, beta); + vector << Scalar(1), essential; + const Matrix<Scalar, Size, Size> h = Matrix<Scalar, Size, Size>::Identity() - tau * vector * vector.adjoint(); + auto block = storage.block(1, 1, Side == OnTheLeft ? Size : n - 1, Side == OnTheLeft ? n - 1 : Size); + const Mat original = block; + Mat expected; + if (Side == OnTheLeft) { + expected = h * original; + block.applyHouseholderOnTheLeft(essential, tau, workspace.data()); + } else { + expected = original * h; + block.applyHouseholderOnTheRight(essential, tau, workspace.data()); + } + if (!((block - expected).norm() <= Scalar(16 * Size) * NumTraits<Scalar>::epsilon() * original.norm())) { + state.SkipWithError("Householder application failed"); + return; + } + for (auto _ : state) { + if (Side == OnTheLeft) + block.applyHouseholderOnTheLeft(essential, tau, workspace.data()); + else + block.applyHouseholderOnTheRight(essential, tau, workspace.data()); + benchmark::DoNotOptimize(block.data()); + benchmark::ClobberMemory(); + } +} + +BENCHMARK_TEMPLATE(BM_ShortHouseholder, float, 2, OnTheLeft)->ArgsProduct({{32, 500, 512, 768, 1000, 1024}, {0, 8}}); +BENCHMARK_TEMPLATE(BM_ShortHouseholder, float, 3, OnTheLeft)->ArgsProduct({{32, 500, 512, 768, 1000, 1024}, {0, 8}}); +BENCHMARK_TEMPLATE(BM_ShortHouseholder, double, 2, OnTheLeft)->ArgsProduct({{32, 500, 512, 768, 1000, 1024}, {0, 8}}); +BENCHMARK_TEMPLATE(BM_ShortHouseholder, double, 3, OnTheLeft)->ArgsProduct({{32, 500, 512, 768, 1000, 1024}, {0, 8}}); +BENCHMARK_TEMPLATE(BM_ShortHouseholder, float, 3, OnTheRight)->ArgsProduct({{32, 500, 512, 768, 1000, 1024}, {0, 8}}); +BENCHMARK_TEMPLATE(BM_ShortHouseholder, double, 3, OnTheRight)->ArgsProduct({{32, 500, 512, 768, 1000, 1024}, {0, 8}});
diff --git a/benchmarks/LU/CMakeLists.txt b/benchmarks/LU/CMakeLists.txt index b6923e8..45aa594 100644 --- a/benchmarks/LU/CMakeLists.txt +++ b/benchmarks/LU/CMakeLists.txt
@@ -2,4 +2,5 @@ # SPDX-License-Identifier: MPL-2.0 eigen_add_benchmark(bench_lu bench_lu.cpp) +eigen_add_benchmark(bench_lu_stride bench_lu_stride.cpp) eigen_add_benchmark(bench_rcond bench_rcond.cpp)
diff --git a/benchmarks/LU/bench_lu_stride.cpp b/benchmarks/LU/bench_lu_stride.cpp new file mode 100644 index 0000000..b72ef59 --- /dev/null +++ b/benchmarks/LU/bench_lu_stride.cpp
@@ -0,0 +1,44 @@ +// SPDX-FileCopyrightText: The Eigen Authors +// SPDX-License-Identifier: MPL-2.0 + +#include <benchmark/benchmark.h> +#include <Eigen/LU> +#include <cstdlib> + +using namespace Eigen; + +template <typename Scalar, int StorageOrder> +static void BM_PartialPivLUStride(benchmark::State& state) { + using Mat = Matrix<Scalar, Dynamic, Dynamic, StorageOrder>; + const Index n = state.range(0); + const Index stride = n + state.range(1); + std::srand(1); + const Mat a = Mat::Random(n, n); + Matrix<Scalar, Dynamic, 1> storage(stride * n); + Map<Mat, 0, OuterStride<>> work(storage.data(), n, n, OuterStride<>(stride)); + work = a; + PartialPivLU<Ref<Mat>> solver(work); + const Mat lower = work.template triangularView<UnitLower>(); + const Mat upper = work.template triangularView<Upper>(); + const Scalar residual = (solver.permutationP() * a - lower * upper).norm(); + const Scalar bound = Scalar(64 * n) * NumTraits<Scalar>::epsilon() * a.norm(); + if (!(residual <= bound)) { + state.SkipWithError("LU reconstruction failed"); + return; + } + for (auto _ : state) { + state.PauseTiming(); + work = a; + state.ResumeTiming(); + solver.compute(work); + benchmark::DoNotOptimize(work.data()); + benchmark::ClobberMemory(); + } +} + +BENCHMARK_TEMPLATE(BM_PartialPivLUStride, float, ColMajor) + ->ArgsProduct({{32, 128, 500, 512, 768, 1000, 1001, 1024}, {0, 8}}); +BENCHMARK_TEMPLATE(BM_PartialPivLUStride, double, ColMajor) + ->ArgsProduct({{32, 128, 500, 512, 768, 1000, 1001, 1024}, {0, 8}}); +BENCHMARK_TEMPLATE(BM_PartialPivLUStride, float, RowMajor)->ArgsProduct({{32, 128, 500, 512, 1000, 1024}, {0, 8}}); +BENCHMARK_TEMPLATE(BM_PartialPivLUStride, double, RowMajor)->ArgsProduct({{32, 128, 500, 512, 1000, 1024}, {0, 8}});
diff --git a/test/householder.cpp b/test/householder.cpp index e3f9f6f..6c7280b 100644 --- a/test/householder.cpp +++ b/test/householder.cpp
@@ -783,6 +783,222 @@ } } +template <typename Scalar, int Size, int StorageOrder> +void householder_short_strided() { + using Mat = Matrix<Scalar, Dynamic, Dynamic, StorageOrder>; + using RealScalar = typename NumTraits<Scalar>::Real; + using Vector = Matrix<Scalar, Size, 1>; + using Essential = Matrix<Scalar, Size - 1, 1>; + const Vector vector = Vector::Random(); + Essential essential; + Scalar tau; + RealScalar beta; + vector.makeHouseholder(essential, tau, beta); + Vector v; + v << Scalar(1), essential; + const Matrix<Scalar, Size, Size> h = Matrix<Scalar, Size, Size>::Identity() - tau * v * v.adjoint(); + + for (Index cols : {0, 1, 2, 17, 32, 33, 500, 512}) { + for (Index innerStride : {1, 2}) { + for (Index outerStride : {2048, 2055}) { + Matrix<Scalar, Dynamic, 1> storage = + Matrix<Scalar, Dynamic, 1>::Random(outerStride * (StorageOrder == ColMajor ? cols + 2 : Size + 2)); + Map<Mat, 0, Stride<Dynamic, Dynamic>> mapped(storage.data(), Size + 2, cols + 2, + Stride<Dynamic, Dynamic>(outerStride, innerStride)); + const auto originalStorage = storage.eval(); + const Mat original = mapped; + auto block = mapped.block(1, 1, Size, cols); + const Mat expected = h * original.block(1, 1, Size, cols); + Matrix<Scalar, Dynamic, 1> workspace(cols + 2); + workspace.setConstant(Scalar(7)); + // Fixed essential length selects the fused column-major path even with dynamic block dimensions. + block.applyHouseholderOnTheLeft(essential, tau, workspace.data() + 1); + const RealScalar bound = RealScalar(16 * Size) * NumTraits<RealScalar>::epsilon() * original.norm(); + VERIFY((block - expected).norm() <= bound); + VERIFY_IS_EQUAL(workspace[0], Scalar(7)); + VERIFY_IS_EQUAL(workspace[cols + 1], Scalar(7)); + Mat expectedFull = original; + expectedFull.block(1, 1, Size, cols) = block; + const auto resultStorage = storage.eval(); + storage = originalStorage; + mapped = expectedFull; + VERIFY_IS_EQUAL(storage, resultStorage); + + block.applyHouseholderOnTheLeft(essential, Scalar(0), workspace.data() + 1); + VERIFY_IS_EQUAL(storage, resultStorage); + + storage = originalStorage; + const Matrix<Scalar, Dynamic, 1> dynamicEssential = essential; + block.applyHouseholderOnTheLeft(dynamicEssential, tau, workspace.data() + 1); + VERIFY((block - expected).norm() <= bound); + + if (cols > 0) { + storage = originalStorage; + const Scalar aliasedTau = block.coeff(0, 0); + const Matrix<Scalar, Size, Size> aliasedH = + Matrix<Scalar, Size, Size>::Identity() - aliasedTau * v * v.adjoint(); + const Mat aliasedExpected = aliasedH * original.block(1, 1, Size, cols); + block.applyHouseholderOnTheLeft(essential, block.coeffRef(0, 0), workspace.data() + 1); + VERIFY((block - aliasedExpected).norm() <= bound); + } + } + } + } +} + +// Reference for H * a with v = [1; essential] and the operation order of applyHouseholderOnTheLeft, which a +// noncommutative scalar distinguishes: w_j = tau * (a_0j + sum_i conj(e_i) a_(i+1)j), a_0j -= w_j, a_(i+1)j -= e_i w_j. +template <typename MatrixType, typename EssentialType> +MatrixType householder_left_reference(const MatrixType& a, const EssentialType& essential, + const typename MatrixType::Scalar& tau) { + using Scalar = typename MatrixType::Scalar; + MatrixType result = a; + for (Index j = 0; j < a.cols(); ++j) { + Scalar dot = a(0, j); + for (Index i = 0; i < essential.size(); ++i) dot += numext::conj(essential(i)) * a(i + 1, j); + const Scalar w = tau * dot; + result(0, j) -= w; + for (Index i = 0; i < essential.size(); ++i) result(i + 1, j) -= essential(i) * w; + } + return result; +} + +// Whether applyHouseholderOnTheLeft on Derived with this essential type selects the fused column loop. +template <typename Derived, typename EssentialPart> +constexpr bool householder_fused() { + return std::is_same<internal::householder_apply_left_impl<Derived, EssentialPart>, + internal::householder_apply_left_impl<Derived, EssentialPart, true>>::value; +} + +template <bool ExpectFused, typename MatrixType, typename EssentialType> +void verify_householder_left(MatrixType a, const EssentialType& essential, const typename MatrixType::Scalar& tau) { + using RealScalar = typename NumTraits<typename MatrixType::Scalar>::Real; + STATIC_CHECK((householder_fused<MatrixType, EssentialType>() == ExpectFused)); + const MatrixType expected = householder_left_reference(a, essential.eval(), tau); + // Applying a reflector is backward stable; both evaluations carry O(rows * eps * |a|) rounding. + const RealScalar bound = RealScalar(16 * a.rows()) * NumTraits<RealScalar>::epsilon() * a.norm(); + Matrix<typename MatrixType::Scalar, Dynamic, 1> workspace(a.cols()); + a.applyHouseholderOnTheLeft(essential, tau, workspace.data()); + VERIFY((a - expected).norm() <= bound); +} + +// The essential part may be any column-vector expression, not only a plain vector: a runtime one-column block +// (no compile-time vector shape), the adjoint of a row tail as CompleteOrthogonalDecomposition passes, or a +// fixed-size segment of either. Each must compile in C++14, where no branch of the dispatch is discarded, and +// only the fixed-size ones on column-major storage select the fused loop. +template <typename Scalar, int StorageOrder> +void householder_essential_expressions() { + using Mat = Matrix<Scalar, Dynamic, Dynamic, StorageOrder>; + using RowMajorMat = Matrix<Scalar, Dynamic, Dynamic, RowMajor>; + using VectorType = Matrix<Scalar, Dynamic, 1>; + constexpr bool kFusedIfFixed = StorageOrder == ColMajor; + const Index cols = internal::random<Index>(1, 20); + for (Index rows : {Index(2), Index(3), internal::random<Index>(4, 20)}) { + const Mat a = Mat::Random(rows, cols); + Scalar tau; + typename NumTraits<Scalar>::Real beta; + VectorType essential(rows - 1); + VectorType::Random(rows).makeHouseholder(essential, tau, beta); + + Mat column(rows + 1, 3); + column.col(1).tail(rows - 1) = essential; + verify_householder_left<false>(a, column.block(2, 1, rows - 1, 1), tau); + + Mat row(3, rows + 1); + row.row(1).tail(rows - 1) = essential.adjoint(); + verify_householder_left<false>(a, row.row(1).tail(rows - 1).adjoint(), tau); + + RowMajorMat rowMajorColumn(rows + 1, 3); + rowMajorColumn.col(1).tail(rows - 1) = essential; + verify_householder_left<false>(a, rowMajorColumn.col(1).tail(rows - 1), tau); + + if (rows == 2) { + verify_householder_left<kFusedIfFixed>(a, row.row(1).template segment<1>(2).adjoint(), tau); + verify_householder_left<kFusedIfFixed>(a, rowMajorColumn.col(1).template tail<1>(), tau); + } + if (rows == 3) { + verify_householder_left<kFusedIfFixed>(a, row.row(1).template segment<2>(2).adjoint(), tau); + verify_householder_left<kFusedIfFixed>(a, rowMajorColumn.col(1).template tail<2>(), tau); + } + } +} + +namespace noncommutative_scalar { + +// Quaternions: multiplication is associative and conjugation reverses products, but ab != ba in general. +struct Quaternion { + double r, i, j, k; + Quaternion(double real = 0, double x = 0, double y = 0, double z = 0) : r(real), i(x), j(y), k(z) {} + Quaternion operator+(const Quaternion& b) const { return Quaternion(r + b.r, i + b.i, j + b.j, k + b.k); } + Quaternion operator-(const Quaternion& b) const { return Quaternion(r - b.r, i - b.i, j - b.j, k - b.k); } + Quaternion operator-() const { return Quaternion(-r, -i, -j, -k); } + Quaternion operator*(const Quaternion& b) const { + return Quaternion(r * b.r - i * b.i - j * b.j - k * b.k, r * b.i + i * b.r + j * b.k - k * b.j, + r * b.j - i * b.k + j * b.r + k * b.i, r * b.k + i * b.j - j * b.i + k * b.r); + } + Quaternion& operator+=(const Quaternion& b) { return *this = *this + b; } + Quaternion& operator-=(const Quaternion& b) { return *this = *this - b; } + Quaternion& operator*=(const Quaternion& b) { return *this = *this * b; } + bool operator==(const Quaternion& b) const { return r == b.r && i == b.i && j == b.j && k == b.k; } + bool operator!=(const Quaternion& b) const { return !(*this == b); } +}; + +inline Quaternion conj(const Quaternion& a) { return Quaternion(a.r, -a.i, -a.j, -a.k); } +inline double real(const Quaternion& a) { return a.r; } +inline double imag(const Quaternion& a) { return a.i; } + +} // namespace noncommutative_scalar + +namespace Eigen { +template <> +struct NumTraits<noncommutative_scalar::Quaternion> : GenericNumTraits<noncommutative_scalar::Quaternion> { + using Real = double; + using Literal = double; + static constexpr bool IsComplex = true; + static constexpr bool RequireInitialization = true; +}; +} // namespace Eigen + +// The fused loop must multiply in the order the general update's expressions spell out, essential * tmp with +// tmp = tau * (essential.adjoint() * column + top), which only a noncommutative scalar observes. Small integer +// components keep every operation exact. The general path is not a reference here: its product kernels reorder +// operands (the gemv selector transposes essential.adjoint() * bottom, the column-major outer-product selector +// forms tmp_j * essential). +template <int Size> +void householder_noncommutative_scalar() { + using Scalar = noncommutative_scalar::Quaternion; + using Mat = Matrix<Scalar, Dynamic, Dynamic>; + using EssentialType = Matrix<Scalar, Size - 1, 1>; + STATIC_CHECK((householder_fused<Mat, EssentialType>())); + const Scalar i(0, 1), j(0, 0, 1), k(0, 0, 0, 1); + { + // With essential = [i, 0...]^T, tau = 1 and input [j, 0...]^T, the second entry becomes -i * j = -k. + Mat a = Mat::Constant(Size, 1, Scalar(0)); + a(0, 0) = j; + EssentialType essential = EssentialType::Constant(Scalar(0)); + essential(0) = i; + Scalar workspace[1]; + a.applyHouseholderOnTheLeft(essential, Scalar(1), workspace); + VERIFY(a(0, 0) == Scalar(0)); + VERIFY(a(1, 0) == -k); + } + const Index cols = internal::random<Index>(1, 7); + Mat a(Size, cols); + for (Index c = 0; c < cols; ++c) + for (Index r = 0; r < Size; ++r) + a(r, c) = Scalar(internal::random<int>(-2, 2), internal::random<int>(-2, 2), internal::random<int>(-2, 2), + internal::random<int>(-2, 2)); + EssentialType essential; + for (Index r = 0; r + 1 < Size; ++r) + essential(r) = Scalar(internal::random<int>(-2, 2), internal::random<int>(-2, 2), internal::random<int>(-2, 2), + internal::random<int>(-2, 2)); + const Scalar tau(internal::random<int>(-2, 2), internal::random<int>(-2, 2), internal::random<int>(-2, 2), 1); + const Mat expected = householder_left_reference(a, essential, tau); + Matrix<Scalar, Dynamic, 1> workspace(cols); + a.applyHouseholderOnTheLeft(essential, tau, workspace.data()); + VERIFY(a == expected); +} + EIGEN_DECLARE_TEST(householder) { for (int i = 0; i < g_repeat; i++) { CALL_SUBTEST_1(householder(Matrix<double, 2, 2>())); @@ -807,4 +1023,23 @@ CALL_SUBTEST_11(householder_blocked_right_regression<std::complex<double>>()); CALL_SUBTEST_12(householder_small_tail()); CALL_SUBTEST_13(householder_large_components()); + CALL_SUBTEST_14((householder_short_strided<float, 2, ColMajor>())); + CALL_SUBTEST_14((householder_short_strided<float, 3, ColMajor>())); + CALL_SUBTEST_15((householder_short_strided<double, 2, ColMajor>())); + CALL_SUBTEST_15((householder_short_strided<double, 3, ColMajor>())); + CALL_SUBTEST_16((householder_short_strided<std::complex<float>, 2, ColMajor>())); + CALL_SUBTEST_16((householder_short_strided<std::complex<float>, 3, ColMajor>())); + CALL_SUBTEST_17((householder_short_strided<std::complex<double>, 2, ColMajor>())); + CALL_SUBTEST_17((householder_short_strided<std::complex<double>, 3, ColMajor>())); + CALL_SUBTEST_18((householder_short_strided<double, 2, RowMajor>())); + CALL_SUBTEST_18((householder_short_strided<double, 3, RowMajor>())); + CALL_SUBTEST_19((householder_short_strided<std::complex<double>, 3, RowMajor>())); + for (int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_20((householder_essential_expressions<double, ColMajor>())); + CALL_SUBTEST_20((householder_essential_expressions<double, RowMajor>())); + CALL_SUBTEST_21((householder_essential_expressions<std::complex<double>, ColMajor>())); + CALL_SUBTEST_21((householder_essential_expressions<std::complex<double>, RowMajor>())); + CALL_SUBTEST_22(householder_noncommutative_scalar<2>()); + CALL_SUBTEST_22(householder_noncommutative_scalar<3>()); + } }
diff --git a/test/lu.cpp b/test/lu.cpp index b5e17f0..72131ac 100644 --- a/test/lu.cpp +++ b/test/lu.cpp
@@ -411,6 +411,36 @@ } } +template <typename Scalar, int StorageOrder> +void lu_strided_pivots() { + using Mat = Matrix<Scalar, Dynamic, Dynamic, StorageOrder>; + using RealScalar = typename NumTraits<Scalar>::Real; + for (Index n : {17, 32, 33, 127, 128, 129, 256, 512}) { + Mat a = Mat::Random(n, n); + // A cyclic permutation forces overlapping row swaps across panel boundaries. + const Index shift = n / 2 + 1; + a.bottomLeftCorner(n - shift, n - shift).diagonal().array() += RealScalar(2 * n); + a.topRightCorner(shift, shift).diagonal().array() += RealScalar(2 * n); + for (Index padding : {0, 7}) { + const Index stride = n + padding; + Matrix<Scalar, Dynamic, 1> storage = Matrix<Scalar, Dynamic, 1>::Constant(stride * n + 2, Scalar(7)); + Map<Mat, 0, OuterStride<>> work(storage.data() + 1, n, n, OuterStride<>(stride)); + work = a; + PartialPivLU<Ref<Mat>> lu(work); + VERIFY_IS_EQUAL(lu.matrixLU().data(), work.data()); + const RealScalar bound = RealScalar(32 * n) * NumTraits<RealScalar>::epsilon() * a.norm(); + const Mat lower = work.template triangularView<UnitLower>(); + const Mat upper = work.template triangularView<Upper>(); + VERIFY((lu.permutationP() * a - lower * upper).norm() <= bound); + VERIFY_IS_EQUAL(lu.permutationP().indices()[0], n - shift); + VERIFY_IS_EQUAL(storage[0], Scalar(7)); + VERIFY_IS_EQUAL(storage[storage.size() - 1], Scalar(7)); + for (Index j = 0; j < n; ++j) + for (Index i = n; i < stride; ++i) VERIFY_IS_EQUAL(storage[1 + j * stride + i], Scalar(7)); + } + } +} + EIGEN_DECLARE_TEST(lu) { for (int i = 0; i < g_repeat; i++) { CALL_SUBTEST_1(lu_non_invertible<Matrix3f>()); @@ -460,6 +490,13 @@ CALL_SUBTEST_9(test_2889()); } + CALL_SUBTEST_10((lu_strided_pivots<float, ColMajor>())); + CALL_SUBTEST_11((lu_strided_pivots<double, ColMajor>())); + CALL_SUBTEST_12((lu_strided_pivots<std::complex<float>, ColMajor>())); + CALL_SUBTEST_13((lu_strided_pivots<std::complex<double>, ColMajor>())); + CALL_SUBTEST_14((lu_strided_pivots<double, RowMajor>())); + CALL_SUBTEST_15((lu_strided_pivots<std::complex<double>, RowMajor>())); + // Blocking and vectorization boundary tests (deterministic, outside g_repeat). CALL_SUBTEST_3(lu_blocking_boundary<float>()); CALL_SUBTEST_4(lu_blocking_boundary<double>());