From a72172e56314feb76f9ffff5a11950c8abdd5ea3 Mon Sep 17 00:00:00 2001 From: Pavel Guzenfeld Date: Fri, 20 Mar 2026 20:27:49 +0000 Subject: [PATCH] Add blocking and vectorization boundary tests for LU and Cholesky libeigen/eigen!2317 --- test/cholesky.cpp | 70 +++++++++++++++++++++++++++++++++++++++++++++++ test/lu.cpp | 70 +++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 140 insertions(+) diff --git a/test/cholesky.cpp b/test/cholesky.cpp index 1c1e18ec9..1c57600e3 100644 --- a/test/cholesky.cpp +++ b/test/cholesky.cpp @@ -461,6 +461,69 @@ void cholesky_verify_assert() { VERIFY_RAISES_ASSERT(ldlt.solveInPlace(tmp)) } +// Test Cholesky decomposition at blocking and vectorization boundaries. +// LLT uses blocks of size max(8, min(size/8 rounded to 16, 128)). +// Sizes near these boundaries exercise the transition between full +// blocked and unblocked paths, including triangular solve boundaries. +template +void cholesky_blocking_boundary() { + typedef typename NumTraits::Real RealScalar; + typedef Matrix MatrixType; + typedef Matrix VectorType; + + const Index PS = internal::packet_traits::size; + const Index sizes[] = {1, 2, 3, PS - 1, PS, PS + 1, 2 * PS - 1, 2 * PS, 2 * PS + 1, 4 * PS, 4 * PS + 1, 7, + 8, 9, 15, 16, 17, 31, 32, 33, 63, 64, 65}; + for (Index si = 0; si < Index(sizeof(sizes) / sizeof(sizes[0])); ++si) { + Index n = sizes[si]; + if (n < 1) continue; + + // Create a symmetric positive definite matrix: A = R'*R + n*I + MatrixType R = MatrixType::Random(n, n); + MatrixType m = R.adjoint() * R; + m.diagonal().array() += RealScalar(n); + + // LLT + LLT llt(m); + VERIFY(llt.info() == Success); + VERIFY_IS_APPROX(m, llt.reconstructedMatrix()); + VectorType rhs = VectorType::Random(n); + VectorType x = llt.solve(rhs); + VERIFY_IS_APPROX(m * x, rhs); + + // LDLT + LDLT ldlt(m); + VERIFY(ldlt.info() == Success); + VERIFY_IS_APPROX(m, ldlt.reconstructedMatrix()); + x = ldlt.solve(rhs); + VERIFY_IS_APPROX(m * x, rhs); + } +} + +// Test Cholesky with RowMajor storage at blocking boundaries. +template +void cholesky_rowmajor_boundary() { + typedef typename NumTraits::Real RealScalar; + typedef Matrix RowMatrixType; + typedef Matrix VectorType; + + const Index sizes[] = {7, 8, 9, 15, 16, 17, 31, 32, 33}; + for (Index si = 0; si < Index(sizeof(sizes) / sizeof(sizes[0])); ++si) { + Index n = sizes[si]; + RowMatrixType R = RowMatrixType::Random(n, n); + RowMatrixType m = R.adjoint() * R; + m.diagonal().array() += RealScalar(n); + + LLT llt(m); + VERIFY(llt.info() == Success); + VERIFY_IS_APPROX(m, llt.reconstructedMatrix()); + + LDLT ldlt(m); + VERIFY(ldlt.info() == Success); + VERIFY_IS_APPROX(m, ldlt.reconstructedMatrix()); + } +} + EIGEN_DECLARE_TEST(cholesky) { int s = 0; for (int i = 0; i < g_repeat; i++) { @@ -496,5 +559,12 @@ EIGEN_DECLARE_TEST(cholesky) { CALL_SUBTEST_2(cholesky_faillure_cases()); + // Blocking and vectorization boundary tests (deterministic, outside g_repeat). + CALL_SUBTEST_2(cholesky_blocking_boundary()); + CALL_SUBTEST_8(cholesky_blocking_boundary()); + CALL_SUBTEST_6(cholesky_blocking_boundary >()); + CALL_SUBTEST_2(cholesky_rowmajor_boundary()); + CALL_SUBTEST_8(cholesky_rowmajor_boundary()); + TEST_SET_BUT_UNUSED_VARIABLE(nb_temporaries) } diff --git a/test/lu.cpp b/test/lu.cpp index 7946c0b4d..a77c86d6a 100644 --- a/test/lu.cpp +++ b/test/lu.cpp @@ -218,6 +218,68 @@ void test_2889() { VERIFY_IS_EQUAL(rcond, 0.0); } +// Test LU decomposition at blocking and vectorization boundaries. +// PartialPivLU uses blocks of size max(8, min(size/8 rounded to 16, 256)). +// Sizes near these boundaries exercise transitions between full blocks +// and remainder tails, including pivot propagation across block edges. +template +void lu_blocking_boundary() { + typedef typename NumTraits::Real RealScalar; + typedef Matrix MatrixType; + + const Index PS = internal::packet_traits::size; + const Index sizes[] = {1, 2, 3, PS - 1, PS, PS + 1, 2 * PS - 1, 2 * PS, 2 * PS + 1, 4 * PS, 4 * PS + 1, 7, + 8, 9, 15, 16, 17, 31, 32, 33, 63, 64, 65}; + for (Index si = 0; si < Index(sizeof(sizes) / sizeof(sizes[0])); ++si) { + Index n = sizes[si]; + if (n < 1) continue; + + // Create a diagonally dominant (invertible) matrix. + MatrixType m = MatrixType::Random(n, n); + m.diagonal().array() += RealScalar(2 * n); + + // PartialPivLU + PartialPivLU plu(m); + VERIFY_IS_APPROX(m, plu.reconstructedMatrix()); + MatrixType rhs = MatrixType::Random(n, 3); + MatrixType x = plu.solve(rhs); + VERIFY_IS_APPROX(m * x, rhs); + + // FullPivLU + FullPivLU flu(m); + VERIFY_IS_APPROX(m, flu.reconstructedMatrix()); + VERIFY(flu.isInvertible()); + x = flu.solve(rhs); + VERIFY_IS_APPROX(m * x, rhs); + } + + // Non-square matrices at boundary sizes for FullPivLU. + const Index rect_sizes[][2] = {{PS, 2 * PS}, {2 * PS, PS}, {15, 33}, {33, 15}, {1, 5}, {5, 1}}; + for (Index si = 0; si < Index(sizeof(rect_sizes) / sizeof(rect_sizes[0])); ++si) { + Index rows = rect_sizes[si][0]; + Index cols = rect_sizes[si][1]; + MatrixType m = MatrixType::Random(rows, cols); + FullPivLU flu(m); + VERIFY_IS_APPROX(m, flu.reconstructedMatrix()); + } +} + +// Test PartialPivLU with RowMajor storage order at blocking boundaries. +template +void lu_rowmajor_boundary() { + typedef typename NumTraits::Real RealScalar; + typedef Matrix RowMatrixType; + + const Index sizes[] = {7, 8, 9, 15, 16, 17, 31, 32, 33}; + for (Index si = 0; si < Index(sizeof(sizes) / sizeof(sizes[0])); ++si) { + Index n = sizes[si]; + RowMatrixType m = RowMatrixType::Random(n, n); + m.diagonal().array() += RealScalar(2 * n); + PartialPivLU plu(m); + VERIFY_IS_APPROX(m, plu.reconstructedMatrix()); + } +} + EIGEN_DECLARE_TEST(lu) { for (int i = 0; i < g_repeat; i++) { CALL_SUBTEST_1(lu_non_invertible()); @@ -257,4 +319,12 @@ EIGEN_DECLARE_TEST(lu) { CALL_SUBTEST_9(test_2889()); } + + // Blocking and vectorization boundary tests (deterministic, outside g_repeat). + CALL_SUBTEST_3(lu_blocking_boundary()); + CALL_SUBTEST_4(lu_blocking_boundary()); + CALL_SUBTEST_5(lu_blocking_boundary >()); + CALL_SUBTEST_6(lu_blocking_boundary >()); + CALL_SUBTEST_4(lu_rowmajor_boundary()); + CALL_SUBTEST_5(lu_rowmajor_boundary >()); }