Core: Fix first_aligned and unaligned complex vectorization libeigen/eigen!3397 Co-authored-by: Antonio Sanchez <cantonios@google.com>
diff --git a/Eigen/src/Core/Redux.h b/Eigen/src/Core/Redux.h index 524dd92..861be0a 100644 --- a/Eigen/src/Core/Redux.h +++ b/Eigen/src/Core/Redux.h
@@ -451,7 +451,7 @@ ? int(packetAlignment) : int(Unaligned); constexpr int alignment = plain_enum_max(alignment0, Evaluator::Alignment); - const Index alignedStart = internal::first_default_aligned(xpr); + const Index alignedStart = internal::first_aligned<alignment>(xpr); const Index alignedSize4 = numext::round_down(size - alignedStart, 4 * packetSize); const Index alignedSize = numext::round_down(size - alignedStart, packetSize); const Index alignedEnd4 = alignedStart + alignedSize4;
diff --git a/Eigen/src/Core/util/Memory.h b/Eigen/src/Core/util/Memory.h index de32b35..75c9f7c 100644 --- a/Eigen/src/Core/util/Memory.h +++ b/Eigen/src/Core/util/Memory.h
@@ -537,13 +537,13 @@ * \param array the address of the start of the array * \param size the size of the array * - * \note If no element of the array is well aligned or the requested alignment is not a multiple of a scalar, - * the size of the array is returned. For example with SSE, the requested alignment is typically 16-bytes. If - * packet size for the given scalar type is 1, then everything is considered well-aligned. + * \note Returns 0 if \a Alignment <= alignof(Scalar). If alignof(Scalar) < \a Alignment <= sizeof(Scalar), all + * elements share the alignment of the first one, so 0 or \a size is returned. Otherwise \a size is returned if no + * element is well aligned or \a Alignment is not a multiple of sizeof(Scalar). * - * \note Otherwise, if the Alignment is larger that the scalar size, we rely on the assumptions that sizeof(Scalar) is a - * power of 2. On the other hand, we do not assume that the array address is a multiple of sizeof(Scalar), as that fails - * for example with Scalar=double on certain 32-bit platforms, see bug #79. + * \note We rely on the assumption that sizeof(Scalar) is a power of 2. On the other hand, we do not assume that the + * array address is a multiple of sizeof(Scalar), as that fails for example with Scalar=double on certain 32-bit + * platforms (see bug #79) or with Scalar=std::complex<double> when alignof(Scalar) < sizeof(Scalar). * * There is also the variant first_aligned(const MatrixBase&) defined in DenseCoeffsBase.h. * \sa first_default_aligned() @@ -554,10 +554,12 @@ constexpr Index AlignmentSize = Alignment / ScalarSize; constexpr Index AlignmentMask = AlignmentSize - 1; - EIGEN_IF_CONSTEXPR (AlignmentSize <= 1) { - // Either the requested alignment if smaller than a scalar, or it exactly match a 1 scalar - // so that all elements of the array have the same alignment. + EIGEN_IF_CONSTEXPR (Alignment <= alignof(Scalar)) { + // A valid Scalar* is already aligned to at least alignof(Scalar). return 0; + } else EIGEN_IF_CONSTEXPR (Alignment <= ScalarSize) { + // Stepping by sizeof(Scalar) preserves address modulo Alignment, so either element 0 is aligned or none is. + return (std::uintptr_t(array) & (Alignment - 1)) ? size : 0; } else EIGEN_IF_CONSTEXPR ((Alignment % ScalarSize) != 0) { // The requested alignment is not a multiple of the scalar size. Consequently, no element of the array is well // aligned.
diff --git a/Eigen/src/ThreadPool/CoreThreadPoolDevice.h b/Eigen/src/ThreadPool/CoreThreadPoolDevice.h index 30be008..28e381a 100644 --- a/Eigen/src/ThreadPool/CoreThreadPoolDevice.h +++ b/Eigen/src/ThreadPool/CoreThreadPoolDevice.h
@@ -296,28 +296,31 @@ struct dense_assignment_loop_with_device<Kernel, CoreThreadPoolDevice, LinearVectorizedTraversal, NoUnrolling> { using Scalar = typename Kernel::Scalar; using PacketType = typename Kernel::PacketType; - static constexpr Index XprEvaluationCost = cost_helper<Kernel>::Cost, - RequestedAlignment = Kernel::AssignmentTraits::LinearRequiredAlignment, - PacketSize = unpacket_traits<PacketType>::size, - DstIsAligned = Kernel::AssignmentTraits::DstAlignment >= RequestedAlignment, - DstAlignment = packet_traits<Scalar>::AlignedOnScalar ? RequestedAlignment - : Kernel::AssignmentTraits::DstAlignment, - SrcAlignment = Kernel::AssignmentTraits::JointAlignment; + static constexpr int PacketSize = unpacket_traits<PacketType>::size; + static constexpr int SrcAlignment = Kernel::AssignmentTraits::JointAlignment; + static constexpr int DstAlignment = plain_enum_max(Kernel::AssignmentTraits::DstAlignment, alignof(Scalar)); + static constexpr int RequestedAlignment = unpacket_traits<PacketType>::alignment; + static constexpr bool Alignable = (DstAlignment >= RequestedAlignment) || + (static_cast<std::size_t>(RequestedAlignment - DstAlignment) % sizeof(Scalar) == 0); + static constexpr int Alignment = Alignable ? RequestedAlignment : DstAlignment; + static constexpr bool DstIsAligned = DstAlignment >= Alignment; + static constexpr bool UsePacketSegment = Kernel::AssignmentTraits::UsePacketSegment; + static constexpr Index XprEvaluationCost = cost_helper<Kernel>::Cost; + struct AssignmentFunctor : public Kernel { EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE AssignmentFunctor(Kernel& kernel) : Kernel(kernel) {} EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void operator()(Index index) { - this->template assignPacket<DstAlignment, SrcAlignment, PacketType>(index); + this->template assignPacket<Alignment, SrcAlignment, PacketType>(index); } }; - static constexpr bool UsePacketSegment = Kernel::AssignmentTraits::UsePacketSegment; + using head_loop = unaligned_dense_assignment_loop<PacketType, DstAlignment, SrcAlignment, UsePacketSegment, DstIsAligned>; - using tail_loop = unaligned_dense_assignment_loop<PacketType, DstAlignment, SrcAlignment, UsePacketSegment, false>; + using tail_loop = unaligned_dense_assignment_loop<PacketType, Alignment, SrcAlignment, UsePacketSegment, false>; static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void run(Kernel& kernel, CoreThreadPoolDevice& device) { const Index size = kernel.size(); - const Index alignedStart = - DstIsAligned ? 0 : internal::first_aligned<RequestedAlignment>(kernel.dstDataPtr(), size); + const Index alignedStart = DstIsAligned ? 0 : internal::first_aligned<Alignment>(kernel.dstDataPtr(), size); const Index alignedEnd = alignedStart + numext::round_down(size - alignedStart, PacketSize); head_loop::run(kernel, 0, alignedStart);
diff --git a/test/assignment_threaded.cpp b/test/assignment_threaded.cpp index b7ffdff..b702616 100644 --- a/test/assignment_threaded.cpp +++ b/test/assignment_threaded.cpp
@@ -44,7 +44,7 @@ PlainObject dst(rows, cols), ref(rows, cols), rhs(rows, cols); rhs.setRandom(); - const auto rhs_xpr = rhs.cwiseAbs2(); + const auto rhs_xpr = rhs + rhs; // linear access dst.setRandom(); @@ -57,6 +57,15 @@ dst.device(threadPoolDevice) = rhs_xpr.unaryExpr(VectorizationOn()); VERIFY_IS_CWISE_EQUAL(ref, dst); + EIGEN_IF_CONSTEXPR (PlainObject::SizeAtCompileTime == Dynamic) { + using RealScalar = typename NumTraits<Scalar>::Real; + Matrix<RealScalar, Dynamic, 1> raw_buf(2 * rows * cols + 2); + Map<PlainObject, Unaligned> unaligned_dst(reinterpret_cast<Scalar*>(raw_buf.data() + 1), rows, cols); + unaligned_dst.setZero(); + unaligned_dst.device(threadPoolDevice) = rhs_xpr.unaryExpr(VectorizationOn()); + VERIFY_IS_CWISE_EQUAL(ref, unaligned_dst); + } + // outer-inner access Index blockRows = numext::maxi(Index(1), rows - 1); Index blockCols = numext::maxi(Index(1), cols - 1); @@ -80,6 +89,8 @@ EIGEN_DECLARE_TEST(assignment_threaded) { for (int i = 0; i < g_repeat; i++) { CALL_SUBTEST(test_threaded_assignment(MatrixXd(), 123, 123)); + CALL_SUBTEST(test_threaded_assignment(MatrixXcf(), 123, 123)); + CALL_SUBTEST(test_threaded_assignment(MatrixXcd(), 123, 123)); CALL_SUBTEST(test_threaded_assignment(Matrix<float, 16, 16>())); } }
diff --git a/test/first_aligned.cpp b/test/first_aligned.cpp index b3b545e..1f8744f 100644 --- a/test/first_aligned.cpp +++ b/test/first_aligned.cpp
@@ -20,7 +20,7 @@ void test_none_aligned_helper(Scalar *array, int size) { EIGEN_UNUSED_VARIABLE(array); EIGEN_UNUSED_VARIABLE(size); - VERIFY(internal::packet_traits<Scalar>::size == 1 || internal::first_default_aligned(array, size) == size); + VERIFY(!internal::packet_traits<Scalar>::Vectorizable || internal::first_default_aligned(array, size) == size); } struct some_non_vectorizable_type { @@ -45,6 +45,18 @@ test_none_aligned_helper(array_double_plus_4_bytes, 50); test_none_aligned_helper(array_double_plus_4_bytes + 1, 50); + // alignof(std::complex<double>) == 8 while sizeof(std::complex<double>) == 16; an 8-byte-offset buffer + // cannot reach >= 16-byte alignment at any element stride. + EIGEN_ALIGN16 std::complex<double> array_cdouble[100]; + test_first_aligned_helper(array_cdouble, 50); + test_first_aligned_helper(array_cdouble + 1, 50); + std::complex<double> *array_cdouble_plus_8_bytes = + reinterpret_cast<std::complex<double> *>(std::uintptr_t(array_cdouble) + 8); + test_none_aligned_helper(array_cdouble_plus_8_bytes, 50); + test_none_aligned_helper(array_cdouble_plus_8_bytes + 1, 50); + VERIFY_IS_EQUAL(internal::first_aligned<16>(array_cdouble_plus_8_bytes, 50), 50); + VERIFY_IS_EQUAL(internal::first_aligned<8>(array_cdouble_plus_8_bytes, 50), 0); + some_non_vectorizable_type array_nonvec[100]; test_first_aligned_helper(array_nonvec, 100); test_none_aligned_helper(array_nonvec, 100);
diff --git a/test/jacobi.cpp b/test/jacobi.cpp index 2616120..0f3e806 100644 --- a/test/jacobi.cpp +++ b/test/jacobi.cpp
@@ -53,6 +53,27 @@ b.applyOnTheRight(p, q, rot); VERIFY_IS_APPROX(b.col(p), c * a.col(p) - s * a.col(q)); VERIFY_IS_APPROX(b.col(q), numext::conj(s) * a.col(p) + numext::conj(c) * a.col(q)); + + EIGEN_IF_CONSTEXPR (RowsAtCompileTime == Dynamic) { + using Scalar = typename MatrixType::Scalar; + using RealScalar = typename NumTraits<Scalar>::Real; + using Vec = Matrix<Scalar, Dynamic, 1>; + Matrix<RealScalar, Dynamic, 1> raw_x(2 * rows + 1); + Matrix<RealScalar, Dynamic, 1> raw_y(2 * rows + 1); + Map<Vec, Unaligned> ux(reinterpret_cast<Scalar*>(raw_x.data() + 1), rows); + Map<Vec, Unaligned> uy(reinterpret_cast<Scalar*>(raw_y.data() + 1), rows); + ux = a.col(p); + uy = a.col(q); + internal::apply_rotation_in_the_plane(ux, uy, rot.transpose()); + VERIFY_IS_APPROX(ux, b.col(p)); + VERIFY_IS_APPROX(uy, b.col(q)); + + ux = a.col(p); + Vec aligned_y = a.col(q); + internal::apply_rotation_in_the_plane(ux, aligned_y, rot.transpose()); + VERIFY_IS_APPROX(ux, b.col(p)); + VERIFY_IS_APPROX(aligned_y, b.col(q)); + } } }
diff --git a/test/product_selfadjoint.cpp b/test/product_selfadjoint.cpp index a7fd13e..465fac3 100644 --- a/test/product_selfadjoint.cpp +++ b/test/product_selfadjoint.cpp
@@ -215,29 +215,44 @@ } } -// Same test for complex type (tests conjugation logic). -template <int> +// Same test for complex types, including buffers aligned only to alignof(RealScalar) < sizeof(Scalar). +template <typename Scalar> void product_selfadjoint_boundary_complex() { - typedef std::complex<float> Scalar; - typedef Matrix<Scalar, Dynamic, Dynamic> Mat; - typedef Matrix<Scalar, Dynamic, 1> Vec; + using RealScalar = typename NumTraits<Scalar>::Real; + using Mat = Matrix<Scalar, Dynamic, Dynamic>; + using Vec = Matrix<Scalar, Dynamic, 1>; + using RealVec = Matrix<RealScalar, Dynamic, 1>; - const int sizes[] = {1, 8, 47, 48, 49, 64}; - for (int si = 0; si < 6; ++si) { - int n = sizes[si]; + const int sizes[] = {1, 2, 3, 4, 8, 9, 16, 47, 48, 49, 64}; + for (int n : sizes) { Mat m1 = Mat::Random(n, n); m1 = (m1 + m1.adjoint()).eval(); // make Hermitian m1.diagonal() = m1.diagonal().real().template cast<Scalar>(); // real diagonal Vec v1 = Vec::Random(n); Mat rhs = Mat::Random(n, 3); + Vec ref = m1 * v1; - Mat m2 = m1.triangularView<Lower>(); - VERIFY_IS_APPROX(m2.selfadjointView<Lower>() * v1, m1 * v1); - VERIFY_IS_APPROX(m2.selfadjointView<Lower>() * rhs, m1 * rhs); + Mat m2 = m1.template triangularView<Lower>(); + VERIFY_IS_APPROX(m2.template selfadjointView<Lower>() * v1, ref); + VERIFY_IS_APPROX(m2.template selfadjointView<Lower>() * rhs, m1 * rhs); - m2 = m1.triangularView<Upper>(); - VERIFY_IS_APPROX(m2.selfadjointView<Upper>() * v1, m1 * v1); + m2 = m1.template triangularView<Upper>(); + VERIFY_IS_APPROX(m2.template selfadjointView<Upper>() * v1, ref); + + RealVec unaligned_res_storage(2 * n + 1); + RealVec unaligned_rhs_storage(2 * n + 1); + Map<Vec, Unaligned> unaligned_res(reinterpret_cast<Scalar*>(unaligned_res_storage.data() + 1), n); + Map<Vec, Unaligned> unaligned_rhs(reinterpret_cast<Scalar*>(unaligned_rhs_storage.data() + 1), n); + unaligned_rhs = v1; + + unaligned_res.setZero(); + unaligned_res.noalias() += m1.template selfadjointView<Lower>() * unaligned_rhs; + VERIFY_IS_APPROX(unaligned_res, ref); + + unaligned_res.setZero(); + unaligned_res.noalias() += m1.template selfadjointView<Upper>() * unaligned_rhs; + VERIFY_IS_APPROX(unaligned_res, ref); } } @@ -267,5 +282,6 @@ // Deterministic blocking boundary tests (outside g_repeat). CALL_SUBTEST_8(product_selfadjoint_boundary<0>()); - CALL_SUBTEST_9(product_selfadjoint_boundary_complex<0>()); + CALL_SUBTEST_9(product_selfadjoint_boundary_complex<std::complex<float>>()); + CALL_SUBTEST_9(product_selfadjoint_boundary_complex<std::complex<double>>()); }
diff --git a/test/redux.cpp b/test/redux.cpp index 7cbb914..81c7f87 100644 --- a/test/redux.cpp +++ b/test/redux.cpp
@@ -256,6 +256,16 @@ Map<const Vec, Unaligned> mapped_for_prod(buf_for_prod.data() + offset, n); redux_vec_boundary_check(mapped, mapped_for_prod); } + { + using RealScalar = typename NumTraits<Scalar>::Real; + Matrix<RealScalar, Dynamic, 1> raw_buf(2 * n + 1); + Matrix<RealScalar, Dynamic, 1> raw_buf_for_prod(2 * n + 1); + Map<Vec, Unaligned> mapped(reinterpret_cast<Scalar*>(raw_buf.data() + 1), n); + Map<Vec, Unaligned> mapped_for_prod(reinterpret_cast<Scalar*>(raw_buf_for_prod.data() + 1), n); + mapped = v; + mapped_for_prod = v_for_prod; + redux_vec_boundary_check(mapped, mapped_for_prod); + } } } @@ -715,6 +725,8 @@ // Integer reductions are already tested by matrixRedux/vectorRedux with clamped values. CALL_SUBTEST_12(redux_vec_boundary<float>()); CALL_SUBTEST_12(redux_vec_boundary<double>()); + CALL_SUBTEST_12(redux_vec_boundary<std::complex<float>>()); + CALL_SUBTEST_12(redux_vec_boundary<std::complex<double>>()); CALL_SUBTEST_12(redux_block_boundary<float>()); CALL_SUBTEST_12(redux_block_boundary<double>()); CALL_SUBTEST_12(redux_block_boundary<std::complex<float>>());