Core: Work around a GCC 15.1/15.2 miscompile of mixed real/complex TRMV libeigen/eigen!3418 Closes #3214 Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>
diff --git a/Eigen/src/Core/products/TriangularMatrixVector.h b/Eigen/src/Core/products/TriangularMatrixVector.h index d7bba5c..cf69b23 100644 --- a/Eigen/src/Core/products/TriangularMatrixVector.h +++ b/Eigen/src/Core/products/TriangularMatrixVector.h
@@ -18,6 +18,16 @@ namespace internal { +// GCC 15.1 and 15.2 at -O3 with AVX2 miscompile run() for a real matrix times a complex vector: the SLP vectorizer +// breaks copies of it specialized for a constant size (#3214). GCC 15.3 is correct. MinGW's GCC 15.2 fails the same +// way, but EIGEN_COMP_GNUC_STRICT leaves MinGW out. +#if (EIGEN_COMP_GNUC_STRICT || (EIGEN_COMP_MINGW && !EIGEN_COMP_CLANG)) && EIGEN_COMP_GNUC >= 1500 && \ + EIGEN_COMP_GNUC < 1530 +#define EIGEN_TRMV_DONT_INLINE EIGEN_DONT_INLINE __attribute__((optimize("no-tree-slp-vectorize"))) +#else +#define EIGEN_TRMV_DONT_INLINE EIGEN_DONT_INLINE +#endif + template <typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar, bool ConjRhs, int StorageOrder, int Version = Specialized> struct triangular_matrix_vector_product; @@ -28,9 +38,9 @@ static constexpr bool IsLower = (Mode & Lower) == Lower; static constexpr bool HasUnitDiag = (Mode & UnitDiag) == UnitDiag; static constexpr bool HasZeroDiag = (Mode & ZeroDiag) == ZeroDiag; - static EIGEN_DONT_INLINE void run(Index _rows, Index _cols, const LhsScalar* lhs_, Index lhsStride, - const RhsScalar* rhs_, Index rhsIncr, ResScalar* res_, Index resIncr, - const RhsScalar& alpha); + static EIGEN_TRMV_DONT_INLINE void run(Index _rows, Index _cols, const LhsScalar* lhs_, Index lhsStride, + const RhsScalar* rhs_, Index rhsIncr, ResScalar* res_, Index resIncr, + const RhsScalar& alpha); }; template <typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar, bool ConjRhs, int Version> @@ -151,9 +161,9 @@ static constexpr bool IsLower = (Mode & Lower) == Lower; static constexpr bool HasUnitDiag = (Mode & UnitDiag) == UnitDiag; static constexpr bool HasZeroDiag = (Mode & ZeroDiag) == ZeroDiag; - static EIGEN_DONT_INLINE void run(Index _rows, Index _cols, const LhsScalar* lhs_, Index lhsStride, - const RhsScalar* rhs_, Index rhsIncr, ResScalar* res_, Index resIncr, - const ResScalar& alpha); + static EIGEN_TRMV_DONT_INLINE void run(Index _rows, Index _cols, const LhsScalar* lhs_, Index lhsStride, + const RhsScalar* rhs_, Index rhsIncr, ResScalar* res_, Index resIncr, + const ResScalar& alpha); }; template <typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar, bool ConjRhs, int Version> @@ -214,6 +224,8 @@ } } +#undef EIGEN_TRMV_DONT_INLINE + /*************************************************************************** * Wrapper to product_triangular_vector ***************************************************************************/
diff --git a/benchmarks/Core/bench_trmv.cpp b/benchmarks/Core/bench_trmv.cpp index ecc0676..8713995 100644 --- a/benchmarks/Core/bench_trmv.cpp +++ b/benchmarks/Core/bench_trmv.cpp
@@ -1,6 +1,7 @@ // Benchmarks for triangular matrix-vector product (TRMV). // -// Tests y += triangularView(A) * x for various modes and sizes. +// Tests y += triangularView(A) * x for various modes and sizes, with real, complex and mixed (real matrix, complex +// vector) scalars, for column-major and row-major A. // Exercises TriangularMatrixVector.h kernel. // SPDX-FileCopyrightText: The Eigen Authors // SPDX-License-Identifier: MPL-2.0 @@ -10,39 +11,58 @@ using namespace Eigen; -template <typename Scalar> +template <typename LhsScalar, typename RhsScalar> double trmvFlops(Index n) { - // TRMV: ~n^2 multiply-adds - return (NumTraits<Scalar>::IsComplex ? 8.0 : 2.0) * n * n; + // TRMV: ~n^2 multiply-adds, of 2 flops (real), 4 (real times complex) or 8 (complex) each + const int complex_operands = int(NumTraits<LhsScalar>::IsComplex) + int(NumTraits<RhsScalar>::IsComplex); + return double(2 << complex_operands) * n * n; } // y = triangularView<Mode>(A) * x -template <typename Scalar, unsigned int Mode> +template <typename LhsScalar, typename RhsScalar, unsigned int Mode, int StorageOrder = ColMajor> static void BM_TRMV(benchmark::State& state) { const Index n = state.range(0); - using Mat = Matrix<Scalar, Dynamic, Dynamic>; - using Vec = Matrix<Scalar, Dynamic, 1>; + using Mat = Matrix<LhsScalar, Dynamic, Dynamic, StorageOrder>; + using Vec = Matrix<RhsScalar, Dynamic, 1>; + using ResVec = Matrix<typename ScalarBinaryOpTraits<LhsScalar, RhsScalar>::ReturnType, Dynamic, 1>; Mat A = Mat::Random(n, n); Vec x = Vec::Random(n); - Vec y(n); + ResVec y(n); for (auto _ : state) { y.noalias() = A.template triangularView<Mode>() * x; benchmark::DoNotOptimize(y.data()); benchmark::ClobberMemory(); } - state.counters["GFLOPS"] = benchmark::Counter(trmvFlops<Scalar>(n), benchmark::Counter::kIsIterationInvariantRate, - benchmark::Counter::kIs1000); + state.counters["GFLOPS"] = benchmark::Counter( + trmvFlops<LhsScalar, RhsScalar>(n), benchmark::Counter::kIsIterationInvariantRate, benchmark::Counter::kIs1000); } +using cfloat = std::complex<float>; +using cdouble = std::complex<double>; + // clang-format off #define TRMV_SIZES ->Arg(8)->Arg(16)->Arg(32)->Arg(64)->Arg(128)->Arg(256)->Arg(512)->Arg(1024)->Arg(2048) -BENCHMARK(BM_TRMV<float, Lower>) TRMV_SIZES ->Name("TRMV_float_Lower"); -BENCHMARK(BM_TRMV<float, Upper>) TRMV_SIZES ->Name("TRMV_float_Upper"); -BENCHMARK(BM_TRMV<float, UnitLower>) TRMV_SIZES ->Name("TRMV_float_UnitLower"); -BENCHMARK(BM_TRMV<float, UnitUpper>) TRMV_SIZES ->Name("TRMV_float_UnitUpper"); -BENCHMARK(BM_TRMV<double, Lower>) TRMV_SIZES ->Name("TRMV_double_Lower"); -BENCHMARK(BM_TRMV<double, Upper>) TRMV_SIZES ->Name("TRMV_double_Upper"); -BENCHMARK(BM_TRMV<double, UnitLower>) TRMV_SIZES ->Name("TRMV_double_UnitLower"); -BENCHMARK(BM_TRMV<double, UnitUpper>) TRMV_SIZES ->Name("TRMV_double_UnitUpper"); +BENCHMARK(BM_TRMV<float, float, Lower>) TRMV_SIZES ->Name("TRMV_float_Lower"); +BENCHMARK(BM_TRMV<float, float, Upper>) TRMV_SIZES ->Name("TRMV_float_Upper"); +BENCHMARK(BM_TRMV<float, float, UnitLower>) TRMV_SIZES ->Name("TRMV_float_UnitLower"); +BENCHMARK(BM_TRMV<float, float, UnitUpper>) TRMV_SIZES ->Name("TRMV_float_UnitUpper"); +BENCHMARK(BM_TRMV<double, double, Lower>) TRMV_SIZES ->Name("TRMV_double_Lower"); +BENCHMARK(BM_TRMV<double, double, Upper>) TRMV_SIZES ->Name("TRMV_double_Upper"); +BENCHMARK(BM_TRMV<double, double, UnitLower>) TRMV_SIZES ->Name("TRMV_double_UnitLower"); +BENCHMARK(BM_TRMV<double, double, UnitUpper>) TRMV_SIZES ->Name("TRMV_double_UnitUpper"); +BENCHMARK(BM_TRMV<cfloat, cfloat, Lower>) TRMV_SIZES ->Name("TRMV_cfloat_Lower"); +BENCHMARK(BM_TRMV<cfloat, cfloat, Upper>) TRMV_SIZES ->Name("TRMV_cfloat_Upper"); +BENCHMARK(BM_TRMV<cdouble, cdouble, Lower>) TRMV_SIZES ->Name("TRMV_cdouble_Lower"); +BENCHMARK(BM_TRMV<cdouble, cdouble, Upper>) TRMV_SIZES ->Name("TRMV_cdouble_Upper"); +BENCHMARK(BM_TRMV<float, cfloat, Lower>) TRMV_SIZES ->Name("TRMV_float_cfloat_Lower"); +BENCHMARK(BM_TRMV<float, cfloat, Upper>) TRMV_SIZES ->Name("TRMV_float_cfloat_Upper"); +BENCHMARK(BM_TRMV<double, cdouble, Lower>) TRMV_SIZES ->Name("TRMV_double_cdouble_Lower"); +BENCHMARK(BM_TRMV<double, cdouble, Upper>) TRMV_SIZES ->Name("TRMV_double_cdouble_Upper"); +BENCHMARK(BM_TRMV<double, double, Lower, RowMajor>) TRMV_SIZES ->Name("TRMV_double_Lower_RowMajor"); +BENCHMARK(BM_TRMV<double, double, Upper, RowMajor>) TRMV_SIZES ->Name("TRMV_double_Upper_RowMajor"); +BENCHMARK(BM_TRMV<cdouble, cdouble, Lower, RowMajor>) TRMV_SIZES ->Name("TRMV_cdouble_Lower_RowMajor"); +BENCHMARK(BM_TRMV<cdouble, cdouble, Upper, RowMajor>) TRMV_SIZES ->Name("TRMV_cdouble_Upper_RowMajor"); +BENCHMARK(BM_TRMV<double, cdouble, Lower, RowMajor>) TRMV_SIZES ->Name("TRMV_double_cdouble_Lower_RowMajor"); +BENCHMARK(BM_TRMV<double, cdouble, Upper, RowMajor>) TRMV_SIZES ->Name("TRMV_double_cdouble_Upper_RowMajor"); #undef TRMV_SIZES // clang-format on
diff --git a/test/product_trmv.cpp b/test/product_trmv.cpp index 067df4b..b7618b9 100644 --- a/test/product_trmv.cpp +++ b/test/product_trmv.cpp
@@ -117,6 +117,43 @@ trmv_negative_stride_destination_shapes<Scalar, RowMajor>(); } +// A real triangular matrix times a complex vector at a fixed size, through the column-major kernel (directly and as the +// transpose of a row-vector product) and the row-major one. The unevaluated s * v.adjoint() reaches the kernel as a +// conjugated vector and a scale factor, as in #3214. Without the workaround, GCC 15.2 at -O3 with AVX2 gets sizes 3 and +// 4 wrong here; which sizes it miscompiles depends on every product in the translation unit. +template <int Size, int Mode> +void trmv_real_times_complex_fixed_size() { + using CD = std::complex<double>; + using RealMat = Matrix<double, Size, Size>; + using ComplexMat = Matrix<CD, Size, Size>; + using ComplexVec = Matrix<CD, Size, 1>; + using ComplexRowVec = Matrix<CD, 1, Size>; + const RealMat m = RealMat::Random(); + const ComplexVec v = ComplexVec::Random(); + const CD s = internal::random<CD>(); + const ComplexMat dense = m.template triangularView<Mode>().toDenseMatrix().template cast<CD>(); + const ComplexMat dense_of_transpose = + m.transpose().template triangularView<Mode>().toDenseMatrix().template cast<CD>(); + + const ComplexVec sv = s * v; + VERIFY(verifyProduct(ComplexVec(m.template triangularView<Mode>() * sv), ComplexVec(dense * sv), dense, sv)); + VERIFY(verifyProduct(ComplexVec(m.transpose().template triangularView<Mode>() * sv), + ComplexVec(dense_of_transpose * sv), dense_of_transpose, sv)); + const ComplexRowVec svh = s * v.adjoint(); + VERIFY(verifyProduct(ComplexRowVec(svh * m.transpose().template triangularView<Mode>()), + ComplexRowVec(svh * dense_of_transpose), svh, dense_of_transpose)); + VERIFY(verifyProduct(ComplexRowVec(s * v.adjoint() * m.transpose().template triangularView<Mode>()), + ComplexRowVec(svh * dense_of_transpose), svh, dense_of_transpose)); +} + +template <int Size> +void trmv_real_times_complex_fixed_size_all_modes() { + trmv_real_times_complex_fixed_size<Size, Lower>(); + trmv_real_times_complex_fixed_size<Size, Upper>(); + trmv_real_times_complex_fixed_size<Size, UnitLower>(); + trmv_real_times_complex_fixed_size<Size, UnitUpper>(); +} + EIGEN_DECLARE_TEST(product_trmv) { int s = 0; for (int i = 0; i < g_repeat; i++) { @@ -132,6 +169,11 @@ s = internal::random<int>(1, EIGEN_TEST_MAX_SIZE); CALL_SUBTEST_6(trmv(Matrix<float, Dynamic, Dynamic, RowMajor>(s, s))); TEST_SET_BUT_UNUSED_VARIABLE(s); + + CALL_SUBTEST_8(trmv_real_times_complex_fixed_size_all_modes<2>()); + CALL_SUBTEST_8(trmv_real_times_complex_fixed_size_all_modes<3>()); + CALL_SUBTEST_8(trmv_real_times_complex_fixed_size_all_modes<4>()); + CALL_SUBTEST_8(trmv_real_times_complex_fixed_size_all_modes<5>()); } CALL_SUBTEST_7(trmv_negative_stride_destination_all<float>()); CALL_SUBTEST_7(trmv_negative_stride_destination_all<double>());