devtools/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-04-10 11:34:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

Add blocking and vectorization boundary tests for LU and Cholesky

Browse Source 
libeigen/eigen!2317
This commit is contained in:
Pavel Guzenfeld
2026-03-20 20:27:49 +00:00
committed by Rasmus Munk Larsen
parent 30128de0e3
commit a72172e563
2 changed files with 140 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

70
test/cholesky.cpp
View File

@@ -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 <typename Scalar>
void cholesky_blocking_boundary() {
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef Matrix<Scalar, Dynamic, Dynamic> MatrixType;
typedef Matrix<Scalar, Dynamic, 1> VectorType;
const Index PS = internal::packet_traits<Scalar>::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<MatrixType> 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<MatrixType> 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 <typename Scalar>
void cholesky_rowmajor_boundary() {
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef Matrix<Scalar, Dynamic, Dynamic, RowMajor> RowMatrixType;
typedef Matrix<Scalar, Dynamic, 1> 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<RowMatrixType> llt(m);
VERIFY(llt.info() == Success);
VERIFY_IS_APPROX(m, llt.reconstructedMatrix());
LDLT<RowMatrixType> 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<void>());
// Blocking and vectorization boundary tests (deterministic, outside g_repeat).
CALL_SUBTEST_2(cholesky_blocking_boundary<double>());
CALL_SUBTEST_8(cholesky_blocking_boundary<float>());
CALL_SUBTEST_6(cholesky_blocking_boundary<std::complex<double> >());
CALL_SUBTEST_2(cholesky_rowmajor_boundary<double>());
CALL_SUBTEST_8(cholesky_rowmajor_boundary<float>());
TEST_SET_BUT_UNUSED_VARIABLE(nb_temporaries)
}

70
test/lu.cpp
View File

@@ -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 <typename Scalar>
void lu_blocking_boundary() {
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef Matrix<Scalar, Dynamic, Dynamic> MatrixType;
const Index PS = internal::packet_traits<Scalar>::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<MatrixType> 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<MatrixType> 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<MatrixType> flu(m);
VERIFY_IS_APPROX(m, flu.reconstructedMatrix());
}
}
// Test PartialPivLU with RowMajor storage order at blocking boundaries.
template <typename Scalar>
void lu_rowmajor_boundary() {
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef Matrix<Scalar, Dynamic, Dynamic, RowMajor> 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<RowMatrixType> 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<Matrix3f>());
@@ -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<float>());
CALL_SUBTEST_4(lu_blocking_boundary<double>());
CALL_SUBTEST_5(lu_blocking_boundary<std::complex<float> >());
CALL_SUBTEST_6(lu_blocking_boundary<std::complex<double> >());
CALL_SUBTEST_4(lu_rowmajor_boundary<double>());
CALL_SUBTEST_5(lu_rowmajor_boundary<std::complex<float> >());
}