diff --git a/test/linearstructure.cpp b/test/linearstructure.cpp index 404b4319c..082fa4739 100644 --- a/test/linearstructure.cpp +++ b/test/linearstructure.cpp @@ -109,6 +109,56 @@ void real_complex(DenseIndex rows = MatrixType::RowsAtCompileTime, DenseIndex co VERIFY(g_called && "matrix - real not properly optimized"); } +// Test linear structure operations between matrices with different storage orders. +// When storage orders differ, vectorization is disabled (StorageOrdersAgree=false in +// AssignEvaluator.h), exercising the scalar fallback path. +template +void linearStructure_mixed_storage() { + const Index PS = internal::packet_traits::size; + // Sizes at vectorization boundaries to expose any mismatch in traversal + const Index sizes[] = {1, PS, PS + 1, 2 * PS, 2 * PS + 1, 4 * PS + 1, 16}; + typedef Matrix ColMat; + typedef Matrix RowMat; + + for (int si = 0; si < 7; ++si) { + Index n = sizes[si]; + if (n <= 0) continue; + ColMat mc = ColMat::Random(n, n); + RowMat mr = RowMat::Random(n, n); + + // ColMajor + RowMajor → ColMajor + ColMat sum_c = mc + mr; + for (Index j = 0; j < n; ++j) + for (Index i = 0; i < n; ++i) VERIFY_IS_APPROX(sum_c(i, j), mc(i, j) + mr(i, j)); + + // ColMajor - RowMajor → ColMajor + ColMat diff_c = mc - mr; + for (Index j = 0; j < n; ++j) + for (Index i = 0; i < n; ++i) VERIFY_IS_APPROX(diff_c(i, j), mc(i, j) - mr(i, j)); + + // RowMajor + ColMajor → RowMajor + RowMat sum_r = mr + mc; + for (Index j = 0; j < n; ++j) + for (Index i = 0; i < n; ++i) VERIFY_IS_APPROX(sum_r(i, j), mr(i, j) + mc(i, j)); + + // Assignment between storage orders + ColMat from_row = mr; + VERIFY_IS_APPROX(from_row, mr); + RowMat from_col = mc; + VERIFY_IS_APPROX(from_col, mc); + + // cwiseProduct with mixed storage + ColMat cwp = mc.cwiseProduct(mr); + for (Index j = 0; j < n; ++j) + for (Index i = 0; i < n; ++i) VERIFY_IS_APPROX(cwp(i, j), mc(i, j) * mr(i, j)); + + // += with mixed storage + ColMat mc2 = mc; + mc2 += mr; + VERIFY_IS_APPROX(mc2, sum_c); + } +} + template void linearstructure_overflow() { // make sure that /=scalar and /scalar do not overflow @@ -147,4 +197,9 @@ EIGEN_DECLARE_TEST(linearstructure) { CALL_SUBTEST_11(real_complex(10, 10)); } CALL_SUBTEST_4(linearstructure_overflow<0>()); + + // Mixed storage order tests (deterministic, outside g_repeat). + CALL_SUBTEST_12(linearStructure_mixed_storage()); + CALL_SUBTEST_12(linearStructure_mixed_storage()); + CALL_SUBTEST_12(linearStructure_mixed_storage>()); } diff --git a/test/mapped_matrix.cpp b/test/mapped_matrix.cpp index 78ae09cb0..71ecb3483 100644 --- a/test/mapped_matrix.cpp +++ b/test/mapped_matrix.cpp @@ -167,12 +167,83 @@ void check_const_correctness(const PlainObjectType&) { // verify that map-to-const don't have LvalueBit typedef std::add_const_t ConstPlainObjectType; - VERIFY(!(internal::traits >::Flags & LvalueBit)); - VERIFY(!(internal::traits >::Flags & LvalueBit)); + VERIFY(!(internal::traits>::Flags & LvalueBit)); + VERIFY(!(internal::traits>::Flags & LvalueBit)); VERIFY(!(Map::Flags & LvalueBit)); VERIFY(!(Map::Flags & LvalueBit)); } +// Test Map with InnerStride at vectorization boundary sizes. +// Strided Maps exercise different traversal paths (SliceVectorized or Default) +// in assignment and reductions. +template +void map_inner_stride_boundary() { + const Index PS = internal::packet_traits::size; + const Index sizes[] = {1, 2, 3, PS - 1, PS, PS + 1, 2 * PS, 2 * PS + 1, 4 * PS, 4 * PS + 1}; + for (int si = 0; si < 10; ++si) { + const Index n = sizes[si]; + if (n <= 0) continue; + typedef Matrix Vec; + // InnerStride<2>: every other element + Vec data = Vec::Random(2 * n); + Map> strided(data.data(), n); + + // Test assignment to/from strided map + Vec dense = strided; + for (Index k = 0; k < n; ++k) VERIFY_IS_APPROX(dense(k), data(2 * k)); + + // Test scalar operations on strided map + Vec result = Scalar(2) * strided; + for (Index k = 0; k < n; ++k) VERIFY_IS_APPROX(result(k), Scalar(2) * data(2 * k)); + + // Test strided map + dense vector + Vec other = Vec::Random(n); + Vec sum_result = strided + other; + for (Index k = 0; k < n; ++k) VERIFY_IS_APPROX(sum_result(k), data(2 * k) + other(k)); + + // Test writing to strided map + Map> strided_dst(data.data(), n); + strided_dst = other; + for (Index k = 0; k < n; ++k) VERIFY_IS_APPROX(data(2 * k), other(k)); + } +} + +// Test Map with OuterStride on matrices at boundary sizes. +template +void map_outer_stride_boundary() { + const Index PS = internal::packet_traits::size; + typedef Matrix Mat; + // Test various inner dimensions around packet size + const Index inner_sizes[] = {1, PS - 1, PS, PS + 1, 2 * PS, 2 * PS + 1}; + const Index outer_stride = 64; // large enough for any inner size + const Index cols = 4; + + for (int si = 0; si < 6; ++si) { + Index rows = inner_sizes[si]; + if (rows <= 0) continue; + typedef Matrix Vec; + Vec data = Vec::Random(outer_stride * cols); + Map> mapped(data.data(), rows, cols, OuterStride<>(outer_stride)); + + // Test that mapped values match expected layout + Mat dense = mapped; + for (Index j = 0; j < cols; ++j) + for (Index i = 0; i < rows; ++i) VERIFY_IS_APPROX(dense(i, j), data(j * outer_stride + i)); + + // Test reduction on mapped matrix + Scalar ref_sum(0); + for (Index j = 0; j < cols; ++j) + for (Index i = 0; i < rows; ++i) ref_sum += data(j * outer_stride + i); + VERIFY_IS_APPROX(mapped.sum(), ref_sum); + + // Test matrix product with mapped matrix + Vec x = Vec::Random(cols); + Vec y = mapped * x; + Vec y_ref = dense * x; + VERIFY_IS_APPROX(y, y_ref); + } +} + EIGEN_DECLARE_TEST(mapped_matrix) { for (int i = 0; i < g_repeat; i++) { CALL_SUBTEST_1(map_class_vector(Matrix())); @@ -197,4 +268,10 @@ EIGEN_DECLARE_TEST(mapped_matrix) { CALL_SUBTEST_9(map_static_methods(VectorXcd(8))); CALL_SUBTEST_10(map_static_methods(VectorXf(12))); } + + // Strided map tests at vectorization boundaries (deterministic, outside g_repeat). + CALL_SUBTEST_12(map_inner_stride_boundary()); + CALL_SUBTEST_12(map_inner_stride_boundary()); + CALL_SUBTEST_13(map_outer_stride_boundary()); + CALL_SUBTEST_13(map_outer_stride_boundary()); } diff --git a/test/triangular.cpp b/test/triangular.cpp index 1a3fbcbc1..210228d1d 100644 --- a/test/triangular.cpp +++ b/test/triangular.cpp @@ -258,6 +258,71 @@ void triangular_rect(const MatrixType& m) { VERIFY_IS_APPROX(m2, m3); } +// Test triangular solve and product at sizes that exercise GEBP blocking. +// The standard test caps at maxsize=20, which never triggers the blocked code paths +// in TriangularSolverMatrix.h (requires size >= 48 with EIGEN_DEBUG_SMALL_PRODUCT_BLOCKS). +template +void triangular_at_blocking_boundaries() { + typedef double Scalar; + typedef Matrix Mat; + typedef Matrix Vec; + + const int sizes[] = {47, 48, 49, 64, 96, 128}; + for (int si = 0; si < 6; ++si) { + int n = sizes[si]; + Mat m1 = Mat::Random(n, n); + // Make well-conditioned: dominant diagonal + for (int i = 0; i < n; ++i) m1(i, i) += Scalar(n); + + Vec v = Vec::Random(n); + Mat rhs = Mat::Random(n, 5); + + // Upper triangular solve with vector + Mat U = m1.triangularView(); + Vec x = m1.triangularView().solve(v); + VERIFY_IS_APPROX(U * x, v); + + // Lower triangular solve with vector + Mat L = m1.triangularView(); + x = m1.triangularView().solve(v); + VERIFY_IS_APPROX(L * x, v); + + // Upper triangular solve with matrix rhs + Mat X = m1.triangularView().solve(rhs); + VERIFY_IS_APPROX(U * X, rhs); + + // Lower triangular solve with matrix rhs + X = m1.triangularView().solve(rhs); + VERIFY_IS_APPROX(L * X, rhs); + + // Triangular product + Mat prod = m1.triangularView() * rhs; + VERIFY_IS_APPROX(prod, U * rhs); + prod = rhs.transpose() * m1.triangularView(); + VERIFY_IS_APPROX(prod, rhs.transpose() * U); + } + + // Also test with float and RowMajor + { + typedef Matrix RMat; + typedef Matrix FVec; + for (int si = 0; si < 6; ++si) { + int n = sizes[si]; + RMat m1 = RMat::Random(n, n); + for (int i = 0; i < n; ++i) m1(i, i) += float(n); + + FVec v = FVec::Random(n); + RMat U = m1.triangularView(); + FVec x = m1.triangularView().solve(v); + VERIFY_IS_APPROX(U * x, v); + + RMat L = m1.triangularView(); + x = m1.triangularView().solve(v); + VERIFY_IS_APPROX(L * x, v); + } + } +} + void bug_159() { Matrix3d m = Matrix3d::Random().triangularView(); EIGEN_UNUSED_VARIABLE(m) @@ -289,4 +354,7 @@ EIGEN_DECLARE_TEST(triangular) { } CALL_SUBTEST_1(bug_159()); + + // Triangular solve/product at blocking boundaries (deterministic, outside g_repeat). + CALL_SUBTEST_11(triangular_at_blocking_boundaries<0>()); }