diff --git a/Eigen/src/Core/arch/Default/Half.h b/Eigen/src/Core/arch/Default/Half.h index fde72a84d..6034f9f1e 100644 --- a/Eigen/src/Core/arch/Default/Half.h +++ b/Eigen/src/Core/arch/Default/Half.h @@ -175,7 +175,13 @@ struct half_base : public __half_raw { #if defined(EIGEN_GPUCC) #if defined(EIGEN_HIPCC) - EIGEN_DEVICE_FUNC _EIGEN_MAYBE_CONSTEXPR half_base(const __half& h) { x = __half_as_ushort(h); } + // Delegate to raw_uint16_to_half, which reinterprets the raw bits for every storage type of + // __half_raw::x. In the host compile phase on platforms with a native fp16 type (e.g. __fp16 on + // arm64), a direct "x = __half_as_ushort(h)" would perform a numeric integer-to-float conversion. + // numext::bit_cast is used to extract the bits because hip_fp16.h defines __half_as_ushort for + // the device compile phase only (host references fail to link against ROCm 6.3). + EIGEN_DEVICE_FUNC _EIGEN_MAYBE_CONSTEXPR half_base(const __half& h) + : __half_raw(raw_uint16_to_half(numext::bit_cast(h))) {} #elif defined(EIGEN_CUDACC) EIGEN_DEVICE_FUNC _EIGEN_MAYBE_CONSTEXPR half_base(const __half& h) : __half_raw(*(__half_raw*)&h) {} #endif @@ -212,10 +218,12 @@ struct half : public half_impl::half_base { #endif #endif -#if EIGEN_HAS_ARM64_FP16 +// In the device compile phase __half_raw is the vendor type, which has no construct_from_rep_tag, +// so these constructors are restricted to the host compile phase and non-GPU builds. +#if EIGEN_HAS_ARM64_FP16 && !defined(EIGEN_GPU_COMPILE_PHASE) explicit EIGEN_DEVICE_FUNC _EIGEN_MAYBE_CONSTEXPR half(__fp16 b) : half(__half_raw(__half_raw::construct_from_rep_tag(), b)) {} -#elif defined(EIGEN_HAS_BUILTIN_FLOAT16) +#elif defined(EIGEN_HAS_BUILTIN_FLOAT16) && !defined(EIGEN_GPU_COMPILE_PHASE) explicit EIGEN_DEVICE_FUNC _EIGEN_MAYBE_CONSTEXPR half(_Float16 b) : half(__half_raw(__half_raw::construct_from_rep_tag(), b)) {} #endif @@ -240,7 +248,10 @@ struct half : public half_impl::half_base { #if defined(EIGEN_HAS_GPU_FP16) && !defined(EIGEN_GPU_COMPILE_PHASE) EIGEN_DEVICE_FUNC operator __half() const { ::__half_raw hr; - hr.x = x; + // raw_half_as_uint16 reinterprets the raw bits for every storage type of __half_raw::x. + // A direct "hr.x = x" would perform a numeric float-to-integer conversion when x has a + // native fp16 type (e.g. __fp16 on arm64), since the vendor ::__half_raw::x is an integer. + hr.x = half_impl::raw_half_as_uint16(*this); return __half(hr); } #endif @@ -616,13 +627,17 @@ EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half operator--(half& a, int) { // also possible to vectorize directly. EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC _EIGEN_MAYBE_CONSTEXPR __half_raw raw_uint16_to_half(numext::uint16_t x) { - // We cannot simply do a "return __half_raw(x)" here, because __half_raw is union type - // in the hip_fp16 header file, and that will trigger a compile error - // On the other hand, having anything but a return statement also triggers a compile error - // because this is constexpr function. + // In the device compile phase of a GPU build we cannot simply do a "return __half_raw(x)", + // because there __half_raw is the vendor type (a union in the hip_fp16 header file) that has + // no uint16 constructor, and that will trigger a compile error. There "h.x = x" assigns the + // raw bits, since the vendor member x is an integer. + // In the host compile phase (and in non-GPU builds) Eigen's own __half_raw is in effect and its + // member x may be a native fp16 type (__fp16 on arm64, _Float16 with AVX512FP16 or riscv-zfh), + // so "h.x = x" would perform a numeric integer-to-float conversion that corrupts the raw bits; + // the explicit uint16 constructor reinterprets the bits for every storage type instead. // Fortunately, since we need to disable EIGEN_CONSTEXPR for GPU anyway, we can get out - // of this catch22 by having separate bodies for GPU / non GPU -#if defined(EIGEN_GPUCC) + // of this catch22 by having separate bodies for the GPU device phase / everything else. +#if defined(EIGEN_GPUCC) && defined(EIGEN_GPU_COMPILE_PHASE) __half_raw h; h.x = x; return h; @@ -893,7 +908,8 @@ struct NumTraits : GenericNumTraits { enum { IsSigned = true, IsInteger = false, IsComplex = false, RequireInitialization = false }; EIGEN_DEVICE_FUNC _EIGEN_MAYBE_CONSTEXPR static EIGEN_STRONG_INLINE Eigen::half epsilon() { - return half_impl::raw_uint16_to_half(0x0800); + // 0x1400 is 2^-10, the fp16 machine epsilon, matching std::numeric_limits::epsilon(). + return half_impl::raw_uint16_to_half(0x1400); } EIGEN_DEVICE_FUNC _EIGEN_MAYBE_CONSTEXPR static EIGEN_STRONG_INLINE Eigen::half dummy_precision() { return half_impl::raw_uint16_to_half(0x211f); // Eigen::half(1e-2f); diff --git a/Eigen/src/Core/arch/RVV10/PacketMathBF16.h b/Eigen/src/Core/arch/RVV10/PacketMathBF16.h index 6d0c47734..14929b013 100644 --- a/Eigen/src/Core/arch/RVV10/PacketMathBF16.h +++ b/Eigen/src/Core/arch/RVV10/PacketMathBF16.h @@ -239,24 +239,40 @@ EIGEN_STRONG_INLINE Packet1Xbf pmax(const Packet1X return F32ToBf16(pmax(Bf16ToF32(a), Bf16ToF32(b))); } +// Comparisons are performed in float32 and the resulting vbool mask is expanded to all-ones/all-zeros +// 16-bit lanes with a vmerge. Narrowing the float32 comparison result with F32ToBf16 instead would go +// through an arithmetic conversion (vfncvtbf16) that canonicalizes the all-ones (NaN) lanes to 0x7fc0 +// and corrupts the mask. template <> EIGEN_STRONG_INLINE Packet1Xbf pcmp_le(const Packet1Xbf& a, const Packet1Xbf& b) { - return F32ToBf16(pcmp_le(Bf16ToF32(a), Bf16ToF32(b))); + PacketMask16 mask = __riscv_vmfle_vv_f32m2_b16(Bf16ToF32(a), Bf16ToF32(b), unpacket_traits::size); + return __riscv_vreinterpret_v_u16m1_bf16m1(__riscv_vmerge_vvm_u16m1( + __riscv_vreinterpret_v_bf16m1_u16m1(pzero(a)), + __riscv_vreinterpret_v_bf16m1_u16m1(ptrue(a)), mask, unpacket_traits::size)); } template <> EIGEN_STRONG_INLINE Packet1Xbf pcmp_lt(const Packet1Xbf& a, const Packet1Xbf& b) { - return F32ToBf16(pcmp_lt(Bf16ToF32(a), Bf16ToF32(b))); + PacketMask16 mask = __riscv_vmflt_vv_f32m2_b16(Bf16ToF32(a), Bf16ToF32(b), unpacket_traits::size); + return __riscv_vreinterpret_v_u16m1_bf16m1(__riscv_vmerge_vvm_u16m1( + __riscv_vreinterpret_v_bf16m1_u16m1(pzero(a)), + __riscv_vreinterpret_v_bf16m1_u16m1(ptrue(a)), mask, unpacket_traits::size)); } template <> EIGEN_STRONG_INLINE Packet1Xbf pcmp_eq(const Packet1Xbf& a, const Packet1Xbf& b) { - return F32ToBf16(pcmp_eq(Bf16ToF32(a), Bf16ToF32(b))); + PacketMask16 mask = __riscv_vmfeq_vv_f32m2_b16(Bf16ToF32(a), Bf16ToF32(b), unpacket_traits::size); + return __riscv_vreinterpret_v_u16m1_bf16m1(__riscv_vmerge_vvm_u16m1( + __riscv_vreinterpret_v_bf16m1_u16m1(pzero(a)), + __riscv_vreinterpret_v_bf16m1_u16m1(ptrue(a)), mask, unpacket_traits::size)); } template <> EIGEN_STRONG_INLINE Packet1Xbf pcmp_lt_or_nan(const Packet1Xbf& a, const Packet1Xbf& b) { - return F32ToBf16(pcmp_lt_or_nan(Bf16ToF32(a), Bf16ToF32(b))); + PacketMask16 mask = __riscv_vmfge_vv_f32m2_b16(Bf16ToF32(a), Bf16ToF32(b), unpacket_traits::size); + return __riscv_vreinterpret_v_u16m1_bf16m1( + __riscv_vmerge_vxm_u16m1(__riscv_vreinterpret_v_bf16m1_u16m1(ptrue(a)), + static_cast(0), mask, unpacket_traits::size)); } EIGEN_STRONG_INLINE Packet1Xbf pselect(const PacketMask16& mask, const Packet1Xbf& a, const Packet1Xbf& b) { @@ -575,24 +591,39 @@ EIGEN_STRONG_INLINE Packet2Xbf pmax(const Packet2X return F32ToBf16(pmax(Bf16ToF32(a), Bf16ToF32(b))); } +// See the Packet1Xbf comparisons above: the vbool mask is expanded with a vmerge because an arithmetic +// narrowing conversion (vfncvtbf16) would canonicalize the all-ones (NaN) lanes to 0x7fc0 and corrupt +// the mask. template <> EIGEN_STRONG_INLINE Packet2Xbf pcmp_le(const Packet2Xbf& a, const Packet2Xbf& b) { - return F32ToBf16(pcmp_le(Bf16ToF32(a), Bf16ToF32(b))); + PacketMask8 mask = __riscv_vmfle_vv_f32m4_b8(Bf16ToF32(a), Bf16ToF32(b), unpacket_traits::size); + return __riscv_vreinterpret_v_u16m2_bf16m2(__riscv_vmerge_vvm_u16m2( + __riscv_vreinterpret_v_bf16m2_u16m2(pzero(a)), + __riscv_vreinterpret_v_bf16m2_u16m2(ptrue(a)), mask, unpacket_traits::size)); } template <> EIGEN_STRONG_INLINE Packet2Xbf pcmp_lt(const Packet2Xbf& a, const Packet2Xbf& b) { - return F32ToBf16(pcmp_lt(Bf16ToF32(a), Bf16ToF32(b))); + PacketMask8 mask = __riscv_vmflt_vv_f32m4_b8(Bf16ToF32(a), Bf16ToF32(b), unpacket_traits::size); + return __riscv_vreinterpret_v_u16m2_bf16m2(__riscv_vmerge_vvm_u16m2( + __riscv_vreinterpret_v_bf16m2_u16m2(pzero(a)), + __riscv_vreinterpret_v_bf16m2_u16m2(ptrue(a)), mask, unpacket_traits::size)); } template <> EIGEN_STRONG_INLINE Packet2Xbf pcmp_eq(const Packet2Xbf& a, const Packet2Xbf& b) { - return F32ToBf16(pcmp_eq(Bf16ToF32(a), Bf16ToF32(b))); + PacketMask8 mask = __riscv_vmfeq_vv_f32m4_b8(Bf16ToF32(a), Bf16ToF32(b), unpacket_traits::size); + return __riscv_vreinterpret_v_u16m2_bf16m2(__riscv_vmerge_vvm_u16m2( + __riscv_vreinterpret_v_bf16m2_u16m2(pzero(a)), + __riscv_vreinterpret_v_bf16m2_u16m2(ptrue(a)), mask, unpacket_traits::size)); } template <> EIGEN_STRONG_INLINE Packet2Xbf pcmp_lt_or_nan(const Packet2Xbf& a, const Packet2Xbf& b) { - return F32ToBf16(pcmp_lt_or_nan(Bf16ToF32(a), Bf16ToF32(b))); + PacketMask8 mask = __riscv_vmfge_vv_f32m4_b8(Bf16ToF32(a), Bf16ToF32(b), unpacket_traits::size); + return __riscv_vreinterpret_v_u16m2_bf16m2( + __riscv_vmerge_vxm_u16m2(__riscv_vreinterpret_v_bf16m2_u16m2(ptrue(a)), + static_cast(0), mask, unpacket_traits::size)); } EIGEN_STRONG_INLINE Packet2Xbf pselect(const PacketMask8& mask, const Packet2Xbf& a, const Packet2Xbf& b) { diff --git a/Eigen/src/SparseQR/SparseQR.h b/Eigen/src/SparseQR/SparseQR.h index c25e93e72..f89533108 100644 --- a/Eigen/src/SparseQR/SparseQR.h +++ b/Eigen/src/SparseQR/SparseQR.h @@ -449,9 +449,30 @@ void SparseQR::factorize(const MatrixType& mat) { RealScalar pivotThreshold; RealScalar max2Norm = RealScalar(0.0); if (m_useDefaultThreshold) { - for (int j = 0; j < n; j++) max2Norm = numext::maxi(max2Norm, m_pmat.col(j).norm()); - if (max2Norm == RealScalar(0)) max2Norm = RealScalar(1); - pivotThreshold = RealScalar(20 * (m + n)) * max2Norm * NumTraits::epsilon(); + // For scalar types narrower than float (half, bfloat16), compute the default threshold in float: + // the squared column norms and the factor 20*(m+n) overflow a 16-bit type for moderate problem + // sizes, and 20*(m+n)*epsilon() exceeds 1 for m+n as small as 52 (half) or 7 (bfloat16), which + // would reject every pivot and collapse the reported rank to zero. + typedef typename std::conditional<(sizeof(RealScalar) < sizeof(float)), float, RealScalar>::type ThresholdReal; + if (EIGEN_CONST_CONDITIONAL((std::is_same::value))) { + // ThresholdReal == RealScalar (float and wider scalars): keep the historical computation bit-for-bit. + for (int j = 0; j < n; j++) max2Norm = numext::maxi(max2Norm, m_pmat.col(j).norm()); + if (max2Norm == RealScalar(0)) max2Norm = RealScalar(1); + pivotThreshold = RealScalar(20 * (m + n)) * max2Norm * NumTraits::epsilon(); + } else { + ThresholdReal maxColNorm = ThresholdReal(0); + for (int j = 0; j < n; j++) { + ThresholdReal colSquaredNorm = ThresholdReal(0); + for (typename QRMatrixType::InnerIterator it(m_pmat, j); it; ++it) + colSquaredNorm += numext::abs2(ThresholdReal(numext::real(it.value()))) + + numext::abs2(ThresholdReal(numext::imag(it.value()))); + maxColNorm = numext::maxi(maxColNorm, numext::sqrt(colSquaredNorm)); + } + if (maxColNorm == ThresholdReal(0)) maxColNorm = ThresholdReal(1); + // Convert back to the narrow type only once, after the whole threshold has been formed in float. + pivotThreshold = RealScalar(ThresholdReal(20 * (m + n)) * maxColNorm * NumTraits::epsilon()); + max2Norm = RealScalar(maxColNorm); + } } else { pivotThreshold = m_threshold; } diff --git a/ci/build.linux.gitlab-ci.yml b/ci/build.linux.gitlab-ci.yml index f06a09c54..4b735e4ab 100644 --- a/ci/build.linux.gitlab-ci.yml +++ b/ci/build.linux.gitlab-ci.yml @@ -526,6 +526,13 @@ build:linux:riscv64:gcc-15:default: EIGEN_CI_CROSS_INSTALL: g++-15-riscv64-linux-gnu EIGEN_CI_CROSS_C_COMPILER: riscv64-linux-gnu-gcc-15 EIGEN_CI_CROSS_CXX_COMPILER: riscv64-linux-gnu-g++-15 + # Zvfbfmin+Zvfbfwma additionally enable the bfloat16 packets + # (PacketMathBF16.h), which are otherwise never compiled by CI. gcc-only: + # clang 21 crashes mangling any function signature containing a + # fixed-length (riscv_rvv_vector_bits) RVV bfloat16 vector type; fixed in + # clang 22. + EIGEN_CI_ADDITIONAL_ARGS: >- + -DEIGEN_TEST_CUSTOM_CXX_FLAGS=-march=rv64gc_v_zvl256b_zfh_zvfh_zvfbfmin_zvfbfwma;-mrvv-vector-bits=zvl;-DEIGEN_RISCV64_USE_RVV10 build:linux:riscv64:clang-21:default: extends: .build:linux:riscv64 diff --git a/test/bfloat16_float.cpp b/test/bfloat16_float.cpp index a5d81d91c..7edf71aa4 100644 --- a/test/bfloat16_float.cpp +++ b/test/bfloat16_float.cpp @@ -181,6 +181,13 @@ void test_numtraits() { VERIFY(NumTraits::IsSigned); + // The bfloat16 machine epsilon is 2^-7 = 0x3c00, and NumTraits must agree + // with std::numeric_limits. + VERIFY_BFLOAT16_BITS_EQUAL(NumTraits::epsilon(), 0x3c00); + VERIFY_BFLOAT16_BITS_EQUAL(std::numeric_limits::epsilon(), 0x3c00); + VERIFY_IS_EQUAL(numext::bit_cast(NumTraits::epsilon()), + numext::bit_cast(std::numeric_limits::epsilon())); + VERIFY_IS_EQUAL(numext::bit_cast(std::numeric_limits::infinity()), numext::bit_cast(bfloat16(std::numeric_limits::infinity()))); // There is no guarantee that casting a 32-bit NaN to bfloat16 has a precise diff --git a/test/gpu_basic.cu b/test/gpu_basic.cu index 68d026b6b..79eba600f 100644 --- a/test/gpu_basic.cu +++ b/test/gpu_basic.cu @@ -438,6 +438,31 @@ bool verifyIsApproxWithInfsNans(const Type1& a, const Type2& b, return true; } +#if defined(EIGEN_HAS_GPU_FP16) && !defined(EIGEN_GPU_COMPILE_PHASE) +// Host-side check that converting between Eigen::half and the vendor __half type preserves the +// raw bits. This is a regression test for builds where Eigen::half stores a native fp16 type +// (e.g. __fp16 on arm64): the host phase used to perform numeric value conversions instead of +// bit reinterpretations, corrupting every raw-bit constant (NumTraits, numeric_limits, ...). +void test_half_raw_bit_interop() { + const numext::uint16_t raw_bits[] = {0x0000, 0x3c00 /*1*/, 0x7c00 /*inf*/, 0x7e00 /*qNaN*/, 0xfbff /*lowest*/}; + for (int i = 0; i < 5; ++i) { + const numext::uint16_t raw = raw_bits[i]; + const Eigen::half h = numext::bit_cast(raw); + // Eigen::half -> __half must preserve the bits (sizeof(__half) == 2 on both CUDA and HIP). + // Call the conversion operator explicitly: a static_cast would be ambiguous because + // Eigen::half also converts to __half via operator float() and __half(float). + const __half v = h.operator __half(); + VERIFY_IS_EQUAL(numext::bit_cast(v), raw); + // __half -> Eigen::half must preserve the bits as well. + const Eigen::half h2(v); + VERIFY_IS_EQUAL(numext::bit_cast(h2), raw); + } + // Raw-bit constants must survive the host phase of a GPU build. + VERIFY((numext::isinf)(NumTraits::infinity())); + VERIFY((numext::isnan)(NumTraits::quiet_NaN())); +} +#endif + template void test_with_infs_nans(const Kernel& ker, int n, const Input& in, Output& out) { Output out_ref, out_gpu; @@ -509,6 +534,12 @@ EIGEN_DECLARE_TEST(gpu_basic) { // numeric_limits CALL_SUBTEST(test_with_infs_nans(numeric_limits_test(), 1, in, out)); + + // Eigen::half <-> __half raw-bit interop on the host. +#if defined(EIGEN_HAS_GPU_FP16) && !defined(EIGEN_GPU_COMPILE_PHASE) + CALL_SUBTEST(test_half_raw_bit_interop()); +#endif + CALL_SUBTEST(test_custom_less_scalar_minmax()); CALL_SUBTEST(test_float_nan_minmax()); diff --git a/test/half_float.cpp b/test/half_float.cpp index 2f81855d2..39a97c90d 100644 --- a/test/half_float.cpp +++ b/test/half_float.cpp @@ -127,6 +127,14 @@ void test_numtraits() { VERIFY(NumTraits::IsSigned); + // The fp16 machine epsilon is 2^-10 = 0x1400, and NumTraits must agree with + // std::numeric_limits (NumTraits used to return 2^-13, making every + // epsilon-based threshold 8x too tight). + VERIFY_HALF_BITS_EQUAL(NumTraits::epsilon(), 0x1400); + VERIFY_HALF_BITS_EQUAL(std::numeric_limits::epsilon(), 0x1400); + VERIFY_IS_EQUAL(numext::bit_cast(NumTraits::epsilon()), + numext::bit_cast(std::numeric_limits::epsilon())); + VERIFY_IS_EQUAL(numext::bit_cast(std::numeric_limits::infinity()), numext::bit_cast(half(std::numeric_limits::infinity()))); // There is no guarantee that casting a 32-bit NaN to 16-bit has a precise diff --git a/test/packetmath_test_shared.h b/test/packetmath_test_shared.h index 64e7619cb..81f2915f0 100644 --- a/test/packetmath_test_shared.h +++ b/test/packetmath_test_shared.h @@ -98,6 +98,25 @@ bool areEqual(const Scalar* a, const Scalar* b, int size) { return true; } +// Checks that each lane of `data` is a full bit mask: all bits set when `zero_mask[i]` is false, all +// bits cleared when it is true. Packet comparisons must return all-ones masks on true lanes so that +// the result can be consumed bitwise (pselect, pand, pandnot, ...); a lane that is merely nonzero +// (e.g. a NaN produced by an arithmetic float-to-half conversion of the mask) is not enough. +template +bool areFullBitMasks(const Scalar* data, const bool* zero_mask, int size) { + for (int i = 0; i < size; ++i) { + const unsigned char expected = zero_mask[i] ? 0x00 : 0xff; + if (!(bits(data[i]) == expected).all()) { + std::cout << "Mask lane " << i << " should be all-" << (zero_mask[i] ? "zero" : "one") + << " bits, got bytes:" << std::hex; + for (std::size_t k = 0; k < sizeof(Scalar); ++k) std::cout << " 0x" << static_cast(bits(data[i])[k]); + std::cout << std::dec << std::endl; + return false; + } + } + return true; +} + template bool areApprox(const Scalar* a, const Scalar* b, int size, const typename NumTraits::Real& precision) { for (int i = 0; i < size; ++i) { @@ -141,6 +160,11 @@ bool areApprox(const Scalar* a, const Scalar* b, int size, const typename NumTra data_mask[i] = numext::is_exactly_zero(data2[i]); \ } \ VERIFY(test::areEqual(ref_mask, data_mask, PacketSize) && #POP); \ + /* Vectorized comparisons must produce all-ones bit masks on true lanes. Scalar fallbacks return */ \ + /* Scalar(1) and bool packets store `true` (loading 0xff into a bool is UB), so those are exempt. */ \ + if (internal::unpacket_traits::vectorizable && !std::is_same::value) { \ + VERIFY(test::areFullBitMasks(data2, ref_mask, PacketSize) && #POP); \ + } \ } // Checks component-wise for input of size N. All of data1, data2, and ref diff --git a/test/sparseqr.cpp b/test/sparseqr.cpp index 16e7663e6..dee7c346d 100644 --- a/test/sparseqr.cpp +++ b/test/sparseqr.cpp @@ -269,6 +269,40 @@ void test_sparseqr_explicit_threshold_disables_lookahead() { VERIFY_IS_EQUAL(explicit_threshold_solver.colsPermutation().indices()(2), 2); } +template +void test_sparseqr_16bit_default_threshold(int n) { + // Factorize a well-conditioned, full-rank, diagonally dominant tridiagonal + // (4 on the diagonal, -1 off-diagonal) system with the default pivot + // threshold. The default threshold 20*(m+n)*max2Norm*epsilon used to be + // computed entirely in the 16-bit scalar type, where 20*(m+n)*epsilon + // exceeds 1 for m+n as small as 52 (half) or 7 (bfloat16) -- so every pivot + // was rejected and the reported rank collapsed to 0 with info()==Success. + typedef SparseMatrix MatrixType; + typedef Matrix DenseVector; + + MatrixType A(n, n); + A.reserve(VectorXi::Constant(n, 3)); + for (int j = 0; j < n; ++j) { + if (j > 0) A.insert(j - 1, j) = Scalar(-1); + A.insert(j, j) = Scalar(4); + if (j + 1 < n) A.insert(j + 1, j) = Scalar(-1); + } + A.makeCompressed(); + + SparseQR > solver(A); + VERIFY_IS_EQUAL(solver.info(), Success); + VERIFY_IS_EQUAL(solver.rank(), Index(n)); + + const DenseVector b = A * DenseVector::Ones(n); + const DenseVector x = solver.solve(b); + VERIFY_IS_EQUAL(solver.info(), Success); + // Check the relative residual in float: 16-bit norms of length-n vectors are + // unreliable. A is diagonally dominant with cond(A) <= 3, so a small + // residual also implies an accurate solution. + const VectorXf residual = (A * x - b).template cast(); + VERIFY(residual.norm() <= 0.05f * b.template cast().norm()); +} + EIGEN_DECLARE_TEST(sparseqr) { for (int i = 0; i < g_repeat; ++i) { CALL_SUBTEST_1(test_sparseqr_scalar()); @@ -279,4 +313,6 @@ EIGEN_DECLARE_TEST(sparseqr) { CALL_SUBTEST_5(test_sparseqr_tiny_independent_column()); CALL_SUBTEST_6(test_sparseqr_explicit_threshold_disables_lookahead()); CALL_SUBTEST_7(test_sparseqr_lookahead_preserves_needed_weak_direction()); + CALL_SUBTEST_8(test_sparseqr_16bit_default_threshold(400)); + CALL_SUBTEST_8(test_sparseqr_16bit_default_threshold(64)); } diff --git a/unsupported/Eigen/src/Tensor/TensorRandom.h b/unsupported/Eigen/src/Tensor/TensorRandom.h index fc0cfb660..c786fda90 100644 --- a/unsupported/Eigen/src/Tensor/TensorRandom.h +++ b/unsupported/Eigen/src/Tensor/TensorRandom.h @@ -209,6 +209,21 @@ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T RandomToTypeNormal(uint64_t* state, uint return v / u; } +// For 16-bit types, compute the deviate in float and round once. Running the +// rejection algorithm above directly in 16-bit arithmetic truncates the tails +// (|v/u| is limited by the coarse uniform grid) and emits NaN/Inf: the 16-bit +// uniform draw is exactly 0 with probability 2^-10 (half) / 2^-7 (bfloat16), +// so log(u) = -inf poisons the acceptance test and v/u returns +/-inf. +template <> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Eigen::half RandomToTypeNormal(uint64_t* state, uint64_t stream) { + return Eigen::half(RandomToTypeNormal(state, stream)); +} +template <> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Eigen::bfloat16 RandomToTypeNormal(uint64_t* state, + uint64_t stream) { + return Eigen::bfloat16(RandomToTypeNormal(state, stream)); +} + template <> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE std::complex RandomToTypeNormal >(uint64_t* state, uint64_t stream) { diff --git a/unsupported/test/tensor_random.cpp b/unsupported/test/tensor_random.cpp index 0e50b535e..e1fca6dc8 100644 --- a/unsupported/test/tensor_random.cpp +++ b/unsupported/test/tensor_random.cpp @@ -55,6 +55,33 @@ static void test_normal() { VERIFY(all_distinct); } +template +static void test_normal_all_finite(Eigen::Index size) { + // Regression test: the 16-bit uniform draw is exactly 0 with probability + // 2^-10 (half) / 2^-7 (bfloat16). Running the ratio-of-uniforms rejection + // in 16-bit arithmetic let log(0) = -inf poison the acceptance test and + // returned v / 0 = +/-inf (or 0/0 = NaN) at measurable rates. + Tensor vec(size); + vec.template setRandom>(); + Eigen::Index num_not_finite = 0; + for (Eigen::Index i = 0; i < size; ++i) { + if (!(numext::isfinite)(vec(i))) ++num_not_finite; + } + VERIFY_IS_EQUAL(num_not_finite, Eigen::Index(0)); +} + +template +static void test_uniform_range(Eigen::Index size) { + // All uniform draws must lie in [0, 1). + Tensor vec(size); + vec.setRandom(); + Eigen::Index num_out_of_range = 0; + for (Eigen::Index i = 0; i < size; ++i) { + if (!(vec(i) >= Scalar(0.0f) && vec(i) < Scalar(1.0f))) ++num_out_of_range; + } + VERIFY_IS_EQUAL(num_out_of_range, Eigen::Index(0)); +} + struct MyGenerator { MyGenerator() {} MyGenerator(const MyGenerator&) {} @@ -96,5 +123,9 @@ EIGEN_DECLARE_TEST(tensor_random) { CALL_SUBTEST((test_normal())); CALL_SUBTEST((test_default())); CALL_SUBTEST((test_normal())); + CALL_SUBTEST((test_normal_all_finite(Eigen::Index(1) << 21))); + CALL_SUBTEST((test_normal_all_finite(Eigen::Index(1) << 18))); + CALL_SUBTEST((test_uniform_range(Eigen::Index(1) << 16))); + CALL_SUBTEST((test_uniform_range(Eigen::Index(1) << 16))); CALL_SUBTEST(test_custom()); } diff --git a/unsupported/test/tensor_random_gpu.cu b/unsupported/test/tensor_random_gpu.cu index 17bdc4474..e522b82f7 100644 --- a/unsupported/test/tensor_random_gpu.cu +++ b/unsupported/test/tensor_random_gpu.cu @@ -61,6 +61,64 @@ void test_gpu_random_normal() { assert(gpuStreamSynchronize(gpu_device.stream()) == gpuSuccess); } +template +void test_gpu_random_uniform_range(int rows, int cols) { + Tensor out(rows, cols); + out.setZero(); + + std::size_t out_bytes = out.size() * sizeof(Scalar); + + Scalar* d_out; + gpuMalloc((void**)(&d_out), out_bytes); + + Eigen::GpuStreamDevice stream; + Eigen::GpuDevice gpu_device(&stream); + + Eigen::TensorMap > gpu_out(d_out, rows, cols); + + gpu_out.device(gpu_device) = gpu_out.random(); + + assert(gpuMemcpyAsync(out.data(), d_out, out_bytes, gpuMemcpyDeviceToHost, gpu_device.stream()) == gpuSuccess); + assert(gpuStreamSynchronize(gpu_device.stream()) == gpuSuccess); + + // All uniform draws must lie in [0, 1). + int num_out_of_range = 0; + for (int i = 0; i < out.size(); ++i) { + if (!(out.data()[i] >= Scalar(0.0f) && out.data()[i] < Scalar(1.0f))) ++num_out_of_range; + } + VERIFY_IS_EQUAL(num_out_of_range, 0); +} + +template +void test_gpu_random_normal_all_finite(int rows, int cols) { + Tensor out(rows, cols); + out.setZero(); + + std::size_t out_bytes = out.size() * sizeof(Scalar); + + Scalar* d_out; + gpuMalloc((void**)(&d_out), out_bytes); + + Eigen::GpuStreamDevice stream; + Eigen::GpuDevice gpu_device(&stream); + + Eigen::TensorMap > gpu_out(d_out, rows, cols); + + Eigen::internal::NormalRandomGenerator gen(true); + gpu_out.device(gpu_device) = gpu_out.random(gen); + + assert(gpuMemcpyAsync(out.data(), d_out, out_bytes, gpuMemcpyDeviceToHost, gpu_device.stream()) == gpuSuccess); + assert(gpuStreamSynchronize(gpu_device.stream()) == gpuSuccess); + + // Regression test for 16-bit types: the deviate must be computed in float, + // otherwise log(0) in the rejection algorithm emits NaN/Inf. + int num_not_finite = 0; + for (int i = 0; i < out.size(); ++i) { + if (!(numext::isfinite)(out.data()[i])) ++num_not_finite; + } + VERIFY_IS_EQUAL(num_not_finite, 0); +} + static void test_complex() { Tensor, 1> vec(6); vec.setRandom(); @@ -75,5 +133,9 @@ static void test_complex() { EIGEN_DECLARE_TEST(tensor_random_gpu) { CALL_SUBTEST(test_gpu_random_uniform()); CALL_SUBTEST(test_gpu_random_normal()); + CALL_SUBTEST(test_gpu_random_uniform_range(256, 256)); + CALL_SUBTEST(test_gpu_random_uniform_range(256, 256)); + CALL_SUBTEST(test_gpu_random_normal_all_finite(1024, 2048)); + CALL_SUBTEST(test_gpu_random_normal_all_finite(512, 512)); CALL_SUBTEST(test_complex()); }