Add boundary test coverage: stableNorm, LinSpaced, complex GEMV, triangular solve

libeigen/eigen!2291

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

test/nullary.cpp

@@ -293,6 +293,130 @@ void nullary_internal_logic() {
   }
 }
 
+// Test LinSpaced at vectorization boundary sizes.
+// The packetOp in linspaced_op_impl uses mask/select logic to handle
+// the last partial packet (when vector size is not a multiple of PacketSize).
+// This exercises those boundaries with element-by-element verification.
+template <typename Scalar>
+void linspaced_boundary() {
+  typedef typename NumTraits<Scalar>::Real RealScalar;
+  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};
+  typedef Matrix<Scalar, Dynamic, 1> Vec;
+
+  for (int si = 0; si < 11; ++si) {
+    Index n = sizes[si];
+    if (n <= 0) continue;
+
+    Scalar low(1), high(100);
+    Vec v = Vec::LinSpaced(n, low, high);
+
+    // With n==1, LinSpaced returns [high] by design.
+    if (n == 1) {
+      VERIFY_IS_EQUAL(v(0), high);
+    } else {
+      VERIFY_IS_EQUAL(v(0), low);
+      VERIFY_IS_EQUAL(v(n - 1), high);
+
+      // Verify monotonicity.
+      for (Index k = 1; k < n; ++k) {
+        VERIFY(numext::real(v(k)) >= numext::real(v(k - 1)));
+      }
+
+      // Verify against scalar reference computation.
+      for (Index k = 0; k < n; ++k) {
+        Scalar ref = Scalar(low + (high - low) * RealScalar(k) / RealScalar(n - 1));
+        VERIFY_IS_APPROX(v(k), ref);
+      }
+    }
+  }
+
+  // Test the "flip" path: when |high| < |low|, the implementation uses
+  // a reversed computation for better precision. Verify at packet boundaries.
+  for (int si = 0; si < 11; ++si) {
+    Index n = sizes[si];
+    if (n <= 0 || n == 1) continue;  // skip n=1, flip irrelevant for single element
+
+    Scalar low(1000), high(1);  // |high| < |low| triggers flip
+    Vec v = Vec::LinSpaced(n, low, high);
+
+    VERIFY_IS_EQUAL(v(0), low);
+    VERIFY_IS_EQUAL(v(n - 1), high);
+
+    // Verify monotonicity (decreasing).
+    for (Index k = 1; k < n; ++k) {
+      VERIFY(numext::real(v(k)) <= numext::real(v(k - 1)));
+    }
+
+    // Verify against scalar reference.
+    for (Index k = 0; k < n; ++k) {
+      Scalar ref = Scalar(low + (high - low) * RealScalar(k) / RealScalar(n - 1));
+      VERIFY_IS_APPROX(v(k), ref);
+    }
+  }
+}
+
+// Test integer LinSpaced divisor path.
+// When (abs(high - low) + 1) < num_steps, the integer LinSpaced uses
+// a divisor-based formula instead of multiplication. This path is
+// barely covered by existing tests which use random ranges.
+template <int>
+void linspaced_integer_divisor() {
+  typedef Matrix<int, Dynamic, 1> VecI;
+
+  // Case: num_steps much larger than range → triggers divisor path.
+  // LinSpaced(12, 0, 5): 12 steps over range [0,5], so range+1=6, 6 < 12.
+  {
+    VecI v = VecI::LinSpaced(12, 0, 5);
+    VERIFY_IS_EQUAL(v(0), 0);
+    // All values must be in [0, 5].
+    for (Index k = 0; k < 12; ++k) {
+      VERIFY(v(k) >= 0 && v(k) <= 5);
+    }
+    // Must be non-decreasing.
+    for (Index k = 1; k < 12; ++k) {
+      VERIFY(v(k) >= v(k - 1));
+    }
+    // Each integer 0-5 should appear at least once.
+    for (int val = 0; val <= 5; ++val) {
+      bool found = false;
+      for (Index k = 0; k < 12; ++k) {
+        if (v(k) == val) {
+          found = true;
+          break;
+        }
+      }
+      VERIFY(found);
+    }
+  }
+
+  // Case: range exactly divides steps → each value should appear equally.
+  // LinSpaced(20, 0, 3): range+1=4, 20%4==0, so each of 0,1,2,3 appears 5 times.
+  {
+    VecI v = VecI::LinSpaced(20, 0, 3);
+    VERIFY_IS_EQUAL(v(0), 0);
+    for (Index k = 0; k < 20; ++k) {
+      VERIFY(v(k) >= 0 && v(k) <= 3);
+    }
+    for (Index k = 1; k < 20; ++k) {
+      VERIFY(v(k) >= v(k - 1));
+    }
+  }
+
+  // Reverse: LinSpaced(12, 5, 0) should be reverse of LinSpaced(12, 0, 5).
+  {
+    VecI fwd = VecI::LinSpaced(12, 0, 5);
+    VecI rev = VecI::LinSpaced(12, 5, 0);
+    VERIFY_IS_APPROX(fwd, rev.reverse());
+  }
+
+  // Single step: always returns high.
+  {
+    VecI v = VecI::LinSpaced(1, 3, 7);
+    VERIFY_IS_EQUAL(v(0), 7);
+  }
+}
+
 EIGEN_DECLARE_TEST(nullary) {
   CALL_SUBTEST_1(testMatrixType(Matrix2d()));
   CALL_SUBTEST_2(testMatrixType(MatrixXcf(internal::random<int>(1, 300), internal::random<int>(1, 300))));
@@ -318,4 +442,11 @@ EIGEN_DECLARE_TEST(nullary) {
   CALL_SUBTEST_6(bug1630<0>());
   CALL_SUBTEST_9(nullary_overflow<0>());
   CALL_SUBTEST_10(nullary_internal_logic<0>());
+
+  // LinSpaced at vectorization boundaries (deterministic, outside g_repeat).
+  CALL_SUBTEST_11(linspaced_boundary<float>());
+  CALL_SUBTEST_11(linspaced_boundary<double>());
+
+  // Integer LinSpaced divisor path tests.
+  CALL_SUBTEST_12(linspaced_integer_divisor<0>());
 }

test/product_extra.cpp

@@ -569,6 +569,88 @@ void product_custom_scalar_types() {
   }
 }
 
+// Test complex GEMV with all conjugation combinations at sizes that
+// exercise full, half, and quarter packet code paths.
+// The GEMV kernels in GeneralMatrixVector.h use conj_helper at three
+// packet levels. The existing product_extra tests cover conjugation
+// but only at random sizes, never systematically at packet boundaries.
+template <int>
+void gemv_complex_conjugate() {
+  typedef std::complex<float> Scf;
+  typedef std::complex<double> Scd;
+  const Index PS_f = internal::packet_traits<Scf>::size;
+  const Index PS_d = internal::packet_traits<Scd>::size;
+
+  // Sizes chosen to exercise packet boundaries for both float and double.
+  const Index sizes[] = {1, 2, 3, 4, 5, 7, 8, 9, 15, 16, 17, 31, 32, 33};
+
+  for (int si = 0; si < 14; ++si) {
+    Index m = sizes[si];
+    // Test complex<float> GEMV with all conjugation combos.
+    {
+      typedef Matrix<Scf, Dynamic, Dynamic> Mat;
+      typedef Matrix<Scf, Dynamic, 1> Vec;
+      Mat A = Mat::Random(m, m);
+      Vec v = Vec::Random(m);
+      Vec res(m);
+
+      // A * v (no conjugation)
+      res.noalias() = A * v;
+      VERIFY_IS_APPROX(res, (A.eval() * v.eval()).eval());
+
+      // A.conjugate() * v
+      res.noalias() = A.conjugate() * v;
+      VERIFY_IS_APPROX(res, (A.conjugate().eval() * v.eval()).eval());
+
+      // A * v.conjugate()
+      res.noalias() = A * v.conjugate();
+      VERIFY_IS_APPROX(res, (A.eval() * v.conjugate().eval()).eval());
+
+      // A.conjugate() * v.conjugate()
+      res.noalias() = A.conjugate() * v.conjugate();
+      VERIFY_IS_APPROX(res, (A.conjugate().eval() * v.conjugate().eval()).eval());
+
+      // A.adjoint() * v (transpose + conjugate of lhs)
+      Vec res2(m);
+      res2.noalias() = A.adjoint() * v;
+      VERIFY_IS_APPROX(res2, (A.adjoint().eval() * v.eval()).eval());
+
+      // Row-major complex GEMV
+      typedef Matrix<Scf, Dynamic, Dynamic, RowMajor> RMat;
+      RMat B = A;
+      res.noalias() = B * v;
+      VERIFY_IS_APPROX(res, (A.eval() * v.eval()).eval());
+
+      res.noalias() = B.conjugate() * v;
+      VERIFY_IS_APPROX(res, (A.conjugate().eval() * v.eval()).eval());
+    }
+
+    // Test complex<double> GEMV with conjugation.
+    {
+      typedef Matrix<Scd, Dynamic, Dynamic> Mat;
+      typedef Matrix<Scd, Dynamic, 1> Vec;
+      Mat A = Mat::Random(m, m);
+      Vec v = Vec::Random(m);
+      Vec res(m);
+
+      res.noalias() = A.conjugate() * v;
+      VERIFY_IS_APPROX(res, (A.conjugate().eval() * v.eval()).eval());
+
+      res.noalias() = A * v.conjugate();
+      VERIFY_IS_APPROX(res, (A.eval() * v.conjugate().eval()).eval());
+
+      // Non-square: wide matrix × vector (exercises different cols path).
+      Mat C = Mat::Random(m, m + 3);
+      Vec w = Vec::Random(m + 3);
+      Vec res3(m);
+      res3.noalias() = C.conjugate() * w;
+      VERIFY_IS_APPROX(res3, (C.conjugate().eval() * w.eval()).eval());
+    }
+  }
+  (void)PS_f;
+  (void)PS_d;
+}
+
 EIGEN_DECLARE_TEST(product_extra) {
   for (int i = 0; i < g_repeat; i++) {
     CALL_SUBTEST_1(product_extra(
@@ -593,4 +675,7 @@ EIGEN_DECLARE_TEST(product_extra) {
   CALL_SUBTEST_8(aliasing_with_resize<void>());
   CALL_SUBTEST_9(product_custom_scalar_types<0>());
   CALL_SUBTEST_10(test_small_block_correctness<0>());
+
+  // Complex GEMV conjugation at varied sizes (deterministic, outside g_repeat).
+  CALL_SUBTEST_11(gemv_complex_conjugate<0>());
 }

test/product_trsolve.cpp

@@ -108,6 +108,119 @@ void trsolve(int size = Size, int cols = Cols) {
   }
 }
 
+// Test triangular solve with non-unit inner stride at blocking boundary sizes.
+// The scalar fallback path in trsmKernelR (TriangularSolverMatrix.h lines 156-166)
+// is used when OtherInnerStride != 1. The existing bug 1741 test only uses
+// InnerStride=2 at random sizes. This exercises the scalar path at sizes that
+// trigger blocking transitions and tests additional configurations.
+template <int>
+void trsolve_strided_boundary() {
+  typedef double Scalar;
+  typedef Matrix<Scalar, Dynamic, Dynamic> MatrixX;
+
+  const int sizes[] = {1, 2, 3, 4, 8, 12, 16, 24, 32, 47, 48, 49, 64};
+  for (int si = 0; si < 13; ++si) {
+    int n = sizes[si];
+
+    MatrixX lhs = MatrixX::Random(n, n);
+    lhs *= 0.1;
+    lhs.diagonal().array() += 1.0;
+
+    // InnerStride = 2: ColMajor RHS, OnTheLeft, Lower
+    {
+      int cols = 5;
+      MatrixX buffer(2 * n, 2 * cols);
+      Map<MatrixX, 0, Stride<Dynamic, 2> > map(buffer.data(), n, cols, Stride<Dynamic, 2>(2 * n, 2));
+      MatrixX ref(n, cols);
+      buffer.setZero();
+      map.setRandom();
+      ref = map;
+      lhs.triangularView<Lower>().solveInPlace(map);
+      VERIFY_IS_APPROX(lhs.triangularView<Lower>().toDenseMatrix() * MatrixX(map), ref);
+    }
+
+    // InnerStride = 2: Upper triangular
+    {
+      int cols = 5;
+      MatrixX buffer(2 * n, 2 * cols);
+      Map<MatrixX, 0, Stride<Dynamic, 2> > map(buffer.data(), n, cols, Stride<Dynamic, 2>(2 * n, 2));
+      MatrixX ref(n, cols);
+      buffer.setZero();
+      map.setRandom();
+      ref = map;
+      lhs.triangularView<Upper>().solveInPlace(map);
+      VERIFY_IS_APPROX(lhs.triangularView<Upper>().toDenseMatrix() * MatrixX(map), ref);
+    }
+
+    // InnerStride = 2: UnitLower (tests the UnitDiag path without diagonal scaling)
+    {
+      int cols = 3;
+      MatrixX buffer(2 * n, 2 * cols);
+      Map<MatrixX, 0, Stride<Dynamic, 2> > map(buffer.data(), n, cols, Stride<Dynamic, 2>(2 * n, 2));
+      MatrixX ref(n, cols);
+      buffer.setZero();
+      map.setRandom();
+      ref = map;
+      lhs.triangularView<UnitLower>().solveInPlace(map);
+      VERIFY_IS_APPROX(lhs.triangularView<UnitLower>().toDenseMatrix() * MatrixX(map), ref);
+    }
+
+    // InnerStride = 3: Less common stride to exercise the scalar path more thoroughly
+    {
+      int cols = 4;
+      MatrixX buffer(3 * n, 3 * cols);
+      Map<MatrixX, 0, Stride<Dynamic, 3> > map(buffer.data(), n, cols, Stride<Dynamic, 3>(3 * n, 3));
+      MatrixX ref(n, cols);
+      buffer.setZero();
+      map.setRandom();
+      ref = map;
+      lhs.triangularView<Lower>().solveInPlace(map);
+      VERIFY_IS_APPROX(lhs.triangularView<Lower>().toDenseMatrix() * MatrixX(map), ref);
+    }
+
+    // Vector RHS with InnerStride = 2
+    {
+      typedef Matrix<Scalar, Dynamic, 1> VecX;
+      VecX buffer(2 * n);
+      Map<VecX, 0, InnerStride<2> > map(buffer.data(), n, InnerStride<2>(2));
+      buffer.setZero();
+      map.setRandom();
+      VecX ref = map;
+      lhs.triangularView<Lower>().solveInPlace(map);
+      VERIFY_IS_APPROX(lhs.triangularView<Lower>().toDenseMatrix() * VecX(map), ref);
+    }
+  }
+
+  // Complex with non-unit stride: tests conjugation in the scalar fallback path.
+  {
+    typedef std::complex<double> CScalar;
+    typedef Matrix<CScalar, Dynamic, Dynamic> CMatrixX;
+    int n = 32;
+    CMatrixX lhs = CMatrixX::Random(n, n);
+    lhs *= CScalar(0.1);
+    lhs.diagonal().array() += CScalar(1.0);
+
+    int cols = 4;
+    CMatrixX buffer(2 * n, 2 * cols);
+    Map<CMatrixX, 0, Stride<Dynamic, 2> > map(buffer.data(), n, cols, Stride<Dynamic, 2>(2 * n, 2));
+    CMatrixX ref(n, cols);
+
+    // Conjugate Lower
+    buffer.setZero();
+    map.setRandom();
+    ref = map;
+    lhs.conjugate().triangularView<Lower>().solveInPlace(map);
+    VERIFY_IS_APPROX(lhs.conjugate().triangularView<Lower>().toDenseMatrix() * CMatrixX(map), ref);
+
+    // Adjoint Upper
+    buffer.setZero();
+    map.setRandom();
+    ref = map;
+    lhs.adjoint().triangularView<Lower>().solveInPlace(map);
+    VERIFY_IS_APPROX(lhs.adjoint().triangularView<Lower>().toDenseMatrix() * CMatrixX(map), ref);
+  }
+}
+
 EIGEN_DECLARE_TEST(product_trsolve) {
   for (int i = 0; i < g_repeat; i++) {
     // matrices
@@ -134,4 +247,7 @@ EIGEN_DECLARE_TEST(product_trsolve) {
     CALL_SUBTEST_13((trsolve<float, 1, 2>()));
     CALL_SUBTEST_14((trsolve<float, 3, 1>()));
   }
+
+  // Strided solve at blocking boundaries (deterministic, outside g_repeat).
+  CALL_SUBTEST_15(trsolve_strided_boundary<0>());
 }

test/stable_norm.cpp

@@ -236,6 +236,71 @@ void test_hypot() {
   VERIFY((numext::isnan)(numext::hypot(a, nan)));
 }
 
+// Test stableNorm at the 4096-element block boundary.
+// stable_norm_impl_inner_step processes vectors in blocks of 4096.
+// Sizes near this boundary exercise the transition between full blocks
+// and the remainder tail, including scale propagation across blocks.
+template <typename Scalar>
+void stable_norm_block_boundary() {
+  using std::abs;
+  using std::sqrt;
+  typedef typename NumTraits<Scalar>::Real RealScalar;
+  typedef Matrix<Scalar, Dynamic, 1> VecType;
+
+  // Test sizes around the 4096 block boundary.
+  const Index sizes[] = {4095, 4096, 4097, 8191, 8192, 8193, 12288};
+  for (int si = 0; si < 7; ++si) {
+    Index n = sizes[si];
+    VecType v = VecType::Random(n);
+    VERIFY_IS_APPROX(v.stableNorm(), v.norm());
+    VERIFY_IS_APPROX(v.blueNorm(), v.norm());
+  }
+
+  // Test scale transitions across blocks: first block has tiny values,
+  // second block has huge values. This exercises the scale/invScale
+  // update logic when maxCoeff > scale in stable_norm_kernel.
+  {
+    RealScalar tiny = (std::numeric_limits<RealScalar>::min)() * RealScalar(1e4);
+    RealScalar huge_val = (std::numeric_limits<RealScalar>::max)() * RealScalar(1e-4);
+    Index n = 8192;
+    VecType v(n);
+    // First 4096 elements: tiny. Second 4096 elements: huge.
+    v.head(4096).setConstant(Scalar(tiny));
+    v.tail(4096).setConstant(Scalar(huge_val));
+    // The huge part dominates, so the expected norm is sqrt(4096)*huge_val.
+    RealScalar expected = sqrt(RealScalar(4096)) * abs(huge_val);
+    VERIFY_IS_APPROX(v.stableNorm(), expected);
+    VERIFY_IS_APPROX(v.blueNorm(), expected);
+  }
+
+  // Reverse: first block huge, second block tiny.
+  {
+    RealScalar tiny = (std::numeric_limits<RealScalar>::min)() * RealScalar(1e4);
+    RealScalar huge_val = (std::numeric_limits<RealScalar>::max)() * RealScalar(1e-4);
+    Index n = 8192;
+    VecType v(n);
+    v.head(4096).setConstant(Scalar(huge_val));
+    v.tail(4096).setConstant(Scalar(tiny));
+    RealScalar expected = sqrt(RealScalar(4096)) * abs(huge_val);
+    VERIFY_IS_APPROX(v.stableNorm(), expected);
+    VERIFY_IS_APPROX(v.blueNorm(), expected);
+  }
+
+  // Matrix version: columns with different magnitudes.
+  // Scale must propagate correctly across columns.
+  {
+    RealScalar tiny = (std::numeric_limits<RealScalar>::min)() * RealScalar(1e4);
+    RealScalar huge_val = (std::numeric_limits<RealScalar>::max)() * RealScalar(1e-4);
+    typedef Matrix<Scalar, Dynamic, Dynamic> MatType;
+    MatType m(100, 2);
+    m.col(0).setConstant(Scalar(tiny));
+    m.col(1).setConstant(Scalar(huge_val));
+    RealScalar expected = sqrt(RealScalar(100)) * abs(huge_val);
+    VERIFY_IS_APPROX(m.stableNorm(), expected);
+    VERIFY_IS_APPROX(m.blueNorm(), expected);
+  }
+}
+
 EIGEN_DECLARE_TEST(stable_norm) {
   CALL_SUBTEST_1(test_empty());
 
@@ -253,4 +318,8 @@ EIGEN_DECLARE_TEST(stable_norm) {
     CALL_SUBTEST_5(stable_norm(VectorXcd(internal::random<int>(10, 2000))));
     CALL_SUBTEST_6(stable_norm(VectorXcf(internal::random<int>(10, 2000))));
   }
+
+  // Block boundary and scale transition tests (deterministic, outside g_repeat).
+  CALL_SUBTEST_7(stable_norm_block_boundary<float>());
+  CALL_SUBTEST_7(stable_norm_block_boundary<double>());
 }