diff --git a/test/nullary.cpp b/test/nullary.cpp index 1f78287c7..a54c61842 100644 --- a/test/nullary.cpp +++ b/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 +void linspaced_boundary() { + typedef typename NumTraits::Real RealScalar; + 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}; + typedef Matrix 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 +void linspaced_integer_divisor() { + typedef Matrix 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(1, 300), internal::random(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()); + CALL_SUBTEST_11(linspaced_boundary()); + + // Integer LinSpaced divisor path tests. + CALL_SUBTEST_12(linspaced_integer_divisor<0>()); } diff --git a/test/product_extra.cpp b/test/product_extra.cpp index 5cad11a12..1e3c6654e 100644 --- a/test/product_extra.cpp +++ b/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 +void gemv_complex_conjugate() { + typedef std::complex Scf; + typedef std::complex Scd; + const Index PS_f = internal::packet_traits::size; + const Index PS_d = internal::packet_traits::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 GEMV with all conjugation combos. + { + typedef Matrix Mat; + typedef Matrix 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 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 GEMV with conjugation. + { + typedef Matrix Mat; + typedef Matrix 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()); 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>()); } diff --git a/test/product_trsolve.cpp b/test/product_trsolve.cpp index 033aa86ed..c7dfb2581 100644 --- a/test/product_trsolve.cpp +++ b/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 +void trsolve_strided_boundary() { + typedef double Scalar; + typedef Matrix 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 > map(buffer.data(), n, cols, Stride(2 * n, 2)); + MatrixX ref(n, cols); + buffer.setZero(); + map.setRandom(); + ref = map; + lhs.triangularView().solveInPlace(map); + VERIFY_IS_APPROX(lhs.triangularView().toDenseMatrix() * MatrixX(map), ref); + } + + // InnerStride = 2: Upper triangular + { + int cols = 5; + MatrixX buffer(2 * n, 2 * cols); + Map > map(buffer.data(), n, cols, Stride(2 * n, 2)); + MatrixX ref(n, cols); + buffer.setZero(); + map.setRandom(); + ref = map; + lhs.triangularView().solveInPlace(map); + VERIFY_IS_APPROX(lhs.triangularView().toDenseMatrix() * MatrixX(map), ref); + } + + // InnerStride = 2: UnitLower (tests the UnitDiag path without diagonal scaling) + { + int cols = 3; + MatrixX buffer(2 * n, 2 * cols); + Map > map(buffer.data(), n, cols, Stride(2 * n, 2)); + MatrixX ref(n, cols); + buffer.setZero(); + map.setRandom(); + ref = map; + lhs.triangularView().solveInPlace(map); + VERIFY_IS_APPROX(lhs.triangularView().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 > map(buffer.data(), n, cols, Stride(3 * n, 3)); + MatrixX ref(n, cols); + buffer.setZero(); + map.setRandom(); + ref = map; + lhs.triangularView().solveInPlace(map); + VERIFY_IS_APPROX(lhs.triangularView().toDenseMatrix() * MatrixX(map), ref); + } + + // Vector RHS with InnerStride = 2 + { + typedef Matrix VecX; + VecX buffer(2 * n); + Map > map(buffer.data(), n, InnerStride<2>(2)); + buffer.setZero(); + map.setRandom(); + VecX ref = map; + lhs.triangularView().solveInPlace(map); + VERIFY_IS_APPROX(lhs.triangularView().toDenseMatrix() * VecX(map), ref); + } + } + + // Complex with non-unit stride: tests conjugation in the scalar fallback path. + { + typedef std::complex CScalar; + typedef Matrix 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 > map(buffer.data(), n, cols, Stride(2 * n, 2)); + CMatrixX ref(n, cols); + + // Conjugate Lower + buffer.setZero(); + map.setRandom(); + ref = map; + lhs.conjugate().triangularView().solveInPlace(map); + VERIFY_IS_APPROX(lhs.conjugate().triangularView().toDenseMatrix() * CMatrixX(map), ref); + + // Adjoint Upper + buffer.setZero(); + map.setRandom(); + ref = map; + lhs.adjoint().triangularView().solveInPlace(map); + VERIFY_IS_APPROX(lhs.adjoint().triangularView().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())); CALL_SUBTEST_14((trsolve())); } + + // Strided solve at blocking boundaries (deterministic, outside g_repeat). + CALL_SUBTEST_15(trsolve_strided_boundary<0>()); } diff --git a/test/stable_norm.cpp b/test/stable_norm.cpp index 567c862df..c0ef4b652 100644 --- a/test/stable_norm.cpp +++ b/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 +void stable_norm_block_boundary() { + using std::abs; + using std::sqrt; + typedef typename NumTraits::Real RealScalar; + typedef Matrix 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::min)() * RealScalar(1e4); + RealScalar huge_val = (std::numeric_limits::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::min)() * RealScalar(1e4); + RealScalar huge_val = (std::numeric_limits::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::min)() * RealScalar(1e4); + RealScalar huge_val = (std::numeric_limits::max)() * RealScalar(1e-4); + typedef Matrix 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(10, 2000)))); CALL_SUBTEST_6(stable_norm(VectorXcf(internal::random(10, 2000)))); } + + // Block boundary and scale transition tests (deterministic, outside g_repeat). + CALL_SUBTEST_7(stable_norm_block_boundary()); + CALL_SUBTEST_7(stable_norm_block_boundary()); }