Add test coverage for strided maps, triangular blocking, and mixed storage orders

libeigen/eigen!2288

Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>

test/linearstructure.cpp

@@ -109,6 +109,56 @@ void real_complex(DenseIndex rows = MatrixType::RowsAtCompileTime, DenseIndex co
   VERIFY(g_called && "matrix<complex> - 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 <typename Scalar>
+void linearStructure_mixed_storage() {
+  const Index PS = internal::packet_traits<Scalar>::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<Scalar, Dynamic, Dynamic, ColMajor> ColMat;
+  typedef Matrix<Scalar, Dynamic, Dynamic, RowMajor> 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 <int>
 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<ArrayXXcf>(10, 10));
   }
   CALL_SUBTEST_4(linearstructure_overflow<0>());
+
+  // Mixed storage order tests (deterministic, outside g_repeat).
+  CALL_SUBTEST_12(linearStructure_mixed_storage<float>());
+  CALL_SUBTEST_12(linearStructure_mixed_storage<double>());
+  CALL_SUBTEST_12(linearStructure_mixed_storage<std::complex<float>>());
 }

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<PlainObjectType> ConstPlainObjectType;
-  VERIFY(!(internal::traits<Map<ConstPlainObjectType> >::Flags & LvalueBit));
-  VERIFY(!(internal::traits<Map<ConstPlainObjectType, AlignedMax> >::Flags & LvalueBit));
+  VERIFY(!(internal::traits<Map<ConstPlainObjectType>>::Flags & LvalueBit));
+  VERIFY(!(internal::traits<Map<ConstPlainObjectType, AlignedMax>>::Flags & LvalueBit));
   VERIFY(!(Map<ConstPlainObjectType>::Flags & LvalueBit));
   VERIFY(!(Map<ConstPlainObjectType, AlignedMax>::Flags & LvalueBit));
 }
 
+// Test Map with InnerStride at vectorization boundary sizes.
+// Strided Maps exercise different traversal paths (SliceVectorized or Default)
+// in assignment and reductions.
+template <typename Scalar>
+void map_inner_stride_boundary() {
+  const Index PS = internal::packet_traits<Scalar>::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<Scalar, Dynamic, 1> Vec;
+    // InnerStride<2>: every other element
+    Vec data = Vec::Random(2 * n);
+    Map<Vec, 0, InnerStride<2>> 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<Vec, 0, InnerStride<2>> 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 <typename Scalar>
+void map_outer_stride_boundary() {
+  const Index PS = internal::packet_traits<Scalar>::size;
+  typedef Matrix<Scalar, Dynamic, Dynamic> 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<Scalar, Dynamic, 1> Vec;
+    Vec data = Vec::Random(outer_stride * cols);
+    Map<Mat, 0, OuterStride<>> 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<float, 1, 1>()));
@@ -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<float>());
+  CALL_SUBTEST_12(map_inner_stride_boundary<double>());
+  CALL_SUBTEST_13(map_outer_stride_boundary<float>());
+  CALL_SUBTEST_13(map_outer_stride_boundary<double>());
 }

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 <int>
+void triangular_at_blocking_boundaries() {
+  typedef double Scalar;
+  typedef Matrix<Scalar, Dynamic, Dynamic> Mat;
+  typedef Matrix<Scalar, Dynamic, 1> 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<Upper>();
+    Vec x = m1.triangularView<Upper>().solve(v);
+    VERIFY_IS_APPROX(U * x, v);
+
+    // Lower triangular solve with vector
+    Mat L = m1.triangularView<Lower>();
+    x = m1.triangularView<Lower>().solve(v);
+    VERIFY_IS_APPROX(L * x, v);
+
+    // Upper triangular solve with matrix rhs
+    Mat X = m1.triangularView<Upper>().solve(rhs);
+    VERIFY_IS_APPROX(U * X, rhs);
+
+    // Lower triangular solve with matrix rhs
+    X = m1.triangularView<Lower>().solve(rhs);
+    VERIFY_IS_APPROX(L * X, rhs);
+
+    // Triangular product
+    Mat prod = m1.triangularView<Upper>() * rhs;
+    VERIFY_IS_APPROX(prod, U * rhs);
+    prod = rhs.transpose() * m1.triangularView<Upper>();
+    VERIFY_IS_APPROX(prod, rhs.transpose() * U);
+  }
+
+  // Also test with float and RowMajor
+  {
+    typedef Matrix<float, Dynamic, Dynamic, RowMajor> RMat;
+    typedef Matrix<float, Dynamic, 1> 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<Upper>();
+      FVec x = m1.triangularView<Upper>().solve(v);
+      VERIFY_IS_APPROX(U * x, v);
+
+      RMat L = m1.triangularView<Lower>();
+      x = m1.triangularView<Lower>().solve(v);
+      VERIFY_IS_APPROX(L * x, v);
+    }
+  }
+}
+
 void bug_159() {
   Matrix3d m = Matrix3d::Random().triangularView<Lower>();
   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>());
 }