Core bug fixes: products, triangular and self-adjoint views libeigen/eigen!3343 Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>
diff --git a/Eigen/src/Core/Assign_AOCL.h b/Eigen/src/Core/Assign_AOCL.h index 5350e93..3f09ce6 100644 --- a/Eigen/src/Core/Assign_AOCL.h +++ b/Eigen/src/Core/Assign_AOCL.h
@@ -53,7 +53,7 @@ * 3. Dense2Dense assignment context with AOCL-compatible traits * * Dispatch conditions (all must be true): - * - Source and destination have DirectAccessBit (contiguous memory) + * - Source and destination are contiguous arrays (direct access, unit inner stride, linear access) * - Compatible storage orders (both row-major or both column-major) * - Vector size ≥ EIGEN_AOCL_VML_THRESHOLD or Dynamic size * - Supported data type (currently double precision for VRDA) @@ -118,13 +118,20 @@ namespace Eigen { namespace internal { +// The AOCL calls take one contiguous array per operand: unit inner stride and linear access, which a block of a +// larger matrix or a map with an outer stride lacks. +template <typename Xpr> +struct aocl_is_contiguous + : bool_constant<(int(Xpr::Flags) & DirectAccessBit) && int(Xpr::InnerStrideAtCompileTime) == 1 && + (int(evaluator<Xpr>::Flags) & LinearAccessBit)> {}; + // Traits for unary operations. template <typename Dst, typename Src> class aocl_assign_traits { private: enum { - DstHasDirectAccess = !!(Dst::Flags & DirectAccessBit), - SrcHasDirectAccess = !!(Src::Flags & DirectAccessBit), + DstIsContiguous = aocl_is_contiguous<Dst>::value, + SrcIsContiguous = aocl_is_contiguous<Src>::value, StorageOrdersAgree = (int(Dst::IsRowMajor) == int(Src::IsRowMajor)), InnerSize = Dst::IsVectorAtCompileTime ? int(Dst::SizeAtCompileTime) : (Dst::Flags & RowMajorBit) ? int(Dst::ColsAtCompileTime) @@ -134,7 +141,7 @@ public: enum { - EnableAoclVML = DstHasDirectAccess && SrcHasDirectAccess && StorageOrdersAgree && LargeEnough, + EnableAoclVML = DstIsContiguous && SrcIsContiguous && StorageOrdersAgree && LargeEnough, Traversal = LinearTraversal }; }; @@ -144,9 +151,9 @@ class aocl_assign_binary_traits { private: enum { - DstHasDirectAccess = !!(Dst::Flags & DirectAccessBit), - LhsHasDirectAccess = !!(Lhs::Flags & DirectAccessBit), - RhsHasDirectAccess = !!(Rhs::Flags & DirectAccessBit), + DstIsContiguous = aocl_is_contiguous<Dst>::value, + LhsIsContiguous = aocl_is_contiguous<Lhs>::value, + RhsIsContiguous = aocl_is_contiguous<Rhs>::value, StorageOrdersAgree = (int(Dst::IsRowMajor) == int(Lhs::IsRowMajor)) && (int(Dst::IsRowMajor) == int(Rhs::IsRowMajor)), InnerSize = Dst::IsVectorAtCompileTime ? int(Dst::SizeAtCompileTime) @@ -156,9 +163,7 @@ }; public: - enum { - EnableAoclVML = DstHasDirectAccess && LhsHasDirectAccess && RhsHasDirectAccess && StorageOrdersAgree && LargeEnough - }; + enum { EnableAoclVML = DstIsContiguous && LhsIsContiguous && RhsIsContiguous && StorageOrdersAgree && LargeEnough }; }; // Unary operation dispatch for float (scalar fallback). @@ -167,7 +172,8 @@ struct Assignment<DstXprType, CwiseUnaryOp<scalar_##EIGENOP##_op<float>, SrcXprNested>, assign_op<float, float>, \ Dense2Dense, std::enable_if_t<aocl_assign_traits<DstXprType, SrcXprNested>::EnableAoclVML>> { \ typedef CwiseUnaryOp<scalar_##EIGENOP##_op<float>, SrcXprNested> SrcXprType; \ - static void run(DstXprType &dst, const SrcXprType &src, const assign_op<float, float> &) { \ + static void run(DstXprType &dst, const SrcXprType &src, const assign_op<float, float> &func) { \ + resize_if_allowed(dst, src, func); \ eigen_assert(dst.rows() == src.rows() && dst.cols() == src.cols()); \ Eigen::Index n = dst.size(); \ if (n <= 0) return; \ @@ -185,7 +191,8 @@ struct Assignment<DstXprType, CwiseUnaryOp<scalar_##EIGENOP##_op<double>, SrcXprNested>, assign_op<double, double>, \ Dense2Dense, std::enable_if_t<aocl_assign_traits<DstXprType, SrcXprNested>::EnableAoclVML>> { \ typedef CwiseUnaryOp<scalar_##EIGENOP##_op<double>, SrcXprNested> SrcXprType; \ - static void run(DstXprType &dst, const SrcXprType &src, const assign_op<double, double> &) { \ + static void run(DstXprType &dst, const SrcXprType &src, const assign_op<double, double> &func) { \ + resize_if_allowed(dst, src, func); \ eigen_assert(dst.rows() == src.rows() && dst.cols() == src.cols()); \ Eigen::Index n = dst.size(); \ eigen_assert(n <= INT_MAX && "AOCL does not support arrays larger than INT_MAX"); \ @@ -220,7 +227,8 @@ assign_op<float, float>, Dense2Dense, \ std::enable_if_t<aocl_assign_binary_traits<DstXprType, LhsXprNested, RhsXprNested>::EnableAoclVML>> { \ typedef CwiseBinaryOp<scalar_##EIGENOP##_op<float, float>, LhsXprNested, RhsXprNested> SrcXprType; \ - static void run(DstXprType &dst, const SrcXprType &src, const assign_op<float, float> &) { \ + static void run(DstXprType &dst, const SrcXprType &src, const assign_op<float, float> &func) { \ + resize_if_allowed(dst, src, func); \ eigen_assert(dst.rows() == src.rows() && dst.cols() == src.cols()); \ Eigen::Index n = dst.size(); \ if (n <= 0) return; \ @@ -241,7 +249,8 @@ assign_op<double, double>, Dense2Dense, \ std::enable_if_t<aocl_assign_binary_traits<DstXprType, LhsXprNested, RhsXprNested>::EnableAoclVML>> { \ typedef CwiseBinaryOp<scalar_##EIGENOP##_op<double, double>, LhsXprNested, RhsXprNested> SrcXprType; \ - static void run(DstXprType &dst, const SrcXprType &src, const assign_op<double, double> &) { \ + static void run(DstXprType &dst, const SrcXprType &src, const assign_op<double, double> &func) { \ + resize_if_allowed(dst, src, func); \ eigen_assert(dst.rows() == src.rows() && dst.cols() == src.cols()); \ Eigen::Index n = dst.size(); \ eigen_assert(n <= INT_MAX && "AOCL does not support arrays larger than INT_MAX"); \
diff --git a/Eigen/src/Core/BandMatrix.h b/Eigen/src/Core/BandMatrix.h index 77bf0f2..dc0fe09 100644 --- a/Eigen/src/Core/BandMatrix.h +++ b/Eigen/src/Core/BandMatrix.h
@@ -102,18 +102,27 @@ using BuildType = Block<CoefficientsType, 1, DiagonalSize>; using Type = std::conditional_t<Conjugate, CwiseUnaryOp<internal::scalar_conjugate_op<Scalar>, BuildType>, BuildType>; + using ConstBuildType = Block<const CoefficientsType, 1, DiagonalSize>; + using ConstType = std::conditional_t<Conjugate, CwiseUnaryOp<internal::scalar_conjugate_op<Scalar>, ConstBuildType>, + ConstBuildType>; }; - /** \returns a vector expression of the \a N -th sub or super diagonal */ + /** \returns a vector expression of the \a N -th sub or super diagonal; for a self-adjoint matrix that stores + * only the opposite diagonal, the conjugate of that one */ template <int N> inline typename DiagonalIntReturnType<N>::Type diagonal() { - return typename DiagonalIntReturnType<N>::BuildType(coeffs(), supers() - N, (std::max)(0, N), 1, diagonalLength(N)); + constexpr int ActualN = DiagonalIntReturnType<N>::ActualIndex; + return typename DiagonalIntReturnType<N>::Type(typename DiagonalIntReturnType<N>::BuildType( + coeffs(), supers() - ActualN, (std::max)(0, ActualN), 1, diagonalLength(ActualN))); } - /** \returns a vector expression of the \a N -th sub or super diagonal */ + /** \returns a vector expression of the \a N -th sub or super diagonal; for a self-adjoint matrix that stores + * only the opposite diagonal, the conjugate of that one */ template <int N> - inline const typename DiagonalIntReturnType<N>::Type diagonal() const { - return typename DiagonalIntReturnType<N>::BuildType(coeffs(), supers() - N, (std::max)(0, N), 1, diagonalLength(N)); + inline const typename DiagonalIntReturnType<N>::ConstType diagonal() const { + constexpr int ActualN = DiagonalIntReturnType<N>::ActualIndex; + return typename DiagonalIntReturnType<N>::ConstType(typename DiagonalIntReturnType<N>::ConstBuildType( + coeffs(), supers() - ActualN, (std::max)(0, ActualN), 1, diagonalLength(ActualN))); } /** \returns a vector expression of the \a i -th sub or super diagonal */ @@ -302,11 +311,11 @@ explicit TridiagonalMatrix(Index size = Size) : Base(size, size, Options & SelfAdjoint ? 0 : 1, 1) {} inline typename Base::template DiagonalIntReturnType<1>::Type super() { return Base::template diagonal<1>(); } - inline const typename Base::template DiagonalIntReturnType<1>::Type super() const { + inline const typename Base::template DiagonalIntReturnType<1>::ConstType super() const { return Base::template diagonal<1>(); } inline typename Base::template DiagonalIntReturnType<-1>::Type sub() { return Base::template diagonal<-1>(); } - inline const typename Base::template DiagonalIntReturnType<-1>::Type sub() const { + inline const typename Base::template DiagonalIntReturnType<-1>::ConstType sub() const { return Base::template diagonal<-1>(); } };
diff --git a/Eigen/src/Core/ProductEvaluators.h b/Eigen/src/Core/ProductEvaluators.h index ec3775d..7da3add 100644 --- a/Eigen/src/Core/ProductEvaluators.h +++ b/Eigen/src/Core/ProductEvaluators.h
@@ -1907,6 +1907,27 @@ : Base(xpr, xpr.rhs().nestedExpression(), xpr.lhs().diagonal()) {} }; +/** \internal + * dst.noalias() += lhs * rhs and -= for the permutation and transposition products below, which only overwrite + * their destination: the product is evaluated into a temporary, as it is without noalias(). + */ +template <typename Lhs, typename Rhs, typename Derived> +struct product_accumulate_through_temporary { + template <typename Dest> + static EIGEN_DEVICE_FUNC void addTo(Dest& dst, const Lhs& lhs, const Rhs& rhs) { + typename Product<Lhs, Rhs>::PlainObject tmp(lhs.rows(), rhs.cols()); + Derived::evalTo(tmp, lhs, rhs); + dst += tmp; + } + + template <typename Dest> + static EIGEN_DEVICE_FUNC void subTo(Dest& dst, const Lhs& lhs, const Rhs& rhs) { + typename Product<Lhs, Rhs>::PlainObject tmp(lhs.rows(), rhs.cols()); + Derived::evalTo(tmp, lhs, rhs); + dst -= tmp; + } +}; + /*************************************************************************** * Products with permutation matrices ***************************************************************************/ @@ -2003,12 +2024,25 @@ struct permutation_matrix_product<ExpressionType, Side, Transposed, TriangularShape> : permutation_densified_matrix_product<ExpressionType, Side, Transposed> {}; +// The scalar factor of a self-adjoint view, folded in from s * (A.selfadjointView() * P), scales the whole view as +// in the other self-adjoint products, not only its stored triangle. template <typename ExpressionType, int Side, bool Transposed> -struct permutation_matrix_product<ExpressionType, Side, Transposed, SelfAdjointShape> - : permutation_densified_matrix_product<ExpressionType, Side, Transposed> {}; +struct permutation_matrix_product<ExpressionType, Side, Transposed, SelfAdjointShape> { + using DenseType = typename ExpressionType::DenseMatrixType; + + template <typename Dest, typename PermutationType> + static void run(Dest& dst, const PermutationType& perm, const ExpressionType& xpr) { + DenseType dense(xpr.rows(), xpr.cols()); + const typename ExpressionType::Scalar factor = densified_view_operand<ExpressionType>::run(dense, xpr); + EIGEN_IF_CONSTEXPR (blas_traits<typename ExpressionType::MatrixType>::HasScalarFactor) dense *= factor; + permutation_matrix_product<DenseType, Side, Transposed, DenseShape>::run(dst, perm, dense); + } +}; template <typename Lhs, typename Rhs, int ProductTag, typename MatrixShape> -struct generic_product_impl<Lhs, Rhs, PermutationShape, MatrixShape, ProductTag> { +struct generic_product_impl<Lhs, Rhs, PermutationShape, MatrixShape, ProductTag> + : product_accumulate_through_temporary<Lhs, Rhs, + generic_product_impl<Lhs, Rhs, PermutationShape, MatrixShape, ProductTag>> { template <typename Dest> static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs) { permutation_matrix_product<Rhs, OnTheLeft, false, MatrixShape>::run(dst, lhs, rhs); @@ -2016,7 +2050,9 @@ }; template <typename Lhs, typename Rhs, int ProductTag, typename MatrixShape> -struct generic_product_impl<Lhs, Rhs, MatrixShape, PermutationShape, ProductTag> { +struct generic_product_impl<Lhs, Rhs, MatrixShape, PermutationShape, ProductTag> + : product_accumulate_through_temporary<Lhs, Rhs, + generic_product_impl<Lhs, Rhs, MatrixShape, PermutationShape, ProductTag>> { template <typename Dest> static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs) { permutation_matrix_product<Lhs, OnTheRight, false, MatrixShape>::run(dst, rhs, lhs); @@ -2024,7 +2060,9 @@ }; template <typename Lhs, typename Rhs, int ProductTag, typename MatrixShape> -struct generic_product_impl<Inverse<Lhs>, Rhs, PermutationShape, MatrixShape, ProductTag> { +struct generic_product_impl<Inverse<Lhs>, Rhs, PermutationShape, MatrixShape, ProductTag> + : product_accumulate_through_temporary< + Inverse<Lhs>, Rhs, generic_product_impl<Inverse<Lhs>, Rhs, PermutationShape, MatrixShape, ProductTag>> { template <typename Dest> static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void evalTo(Dest& dst, const Inverse<Lhs>& lhs, const Rhs& rhs) { permutation_matrix_product<Rhs, OnTheLeft, true, MatrixShape>::run(dst, lhs.nestedExpression(), rhs); @@ -2032,7 +2070,9 @@ }; template <typename Lhs, typename Rhs, int ProductTag, typename MatrixShape> -struct generic_product_impl<Lhs, Inverse<Rhs>, MatrixShape, PermutationShape, ProductTag> { +struct generic_product_impl<Lhs, Inverse<Rhs>, MatrixShape, PermutationShape, ProductTag> + : product_accumulate_through_temporary< + Lhs, Inverse<Rhs>, generic_product_impl<Lhs, Inverse<Rhs>, MatrixShape, PermutationShape, ProductTag>> { template <typename Dest> static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void evalTo(Dest& dst, const Lhs& lhs, const Inverse<Rhs>& rhs) { permutation_matrix_product<Lhs, OnTheRight, true, MatrixShape>::run(dst, rhs.nestedExpression(), lhs); @@ -2079,7 +2119,9 @@ }; template <typename Lhs, typename Rhs, int ProductTag, typename MatrixShape> -struct generic_product_impl<Lhs, Rhs, TranspositionsShape, MatrixShape, ProductTag> { +struct generic_product_impl<Lhs, Rhs, TranspositionsShape, MatrixShape, ProductTag> + : product_accumulate_through_temporary< + Lhs, Rhs, generic_product_impl<Lhs, Rhs, TranspositionsShape, MatrixShape, ProductTag>> { template <typename Dest> static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs) { transposition_matrix_product<Rhs, OnTheLeft, false, MatrixShape>::run(dst, lhs, rhs); @@ -2087,7 +2129,9 @@ }; template <typename Lhs, typename Rhs, int ProductTag, typename MatrixShape> -struct generic_product_impl<Lhs, Rhs, MatrixShape, TranspositionsShape, ProductTag> { +struct generic_product_impl<Lhs, Rhs, MatrixShape, TranspositionsShape, ProductTag> + : product_accumulate_through_temporary< + Lhs, Rhs, generic_product_impl<Lhs, Rhs, MatrixShape, TranspositionsShape, ProductTag>> { template <typename Dest> static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs) { transposition_matrix_product<Lhs, OnTheRight, false, MatrixShape>::run(dst, rhs, lhs); @@ -2095,7 +2139,10 @@ }; template <typename Lhs, typename Rhs, int ProductTag, typename MatrixShape> -struct generic_product_impl<Transpose<Lhs>, Rhs, TranspositionsShape, MatrixShape, ProductTag> { +struct generic_product_impl<Transpose<Lhs>, Rhs, TranspositionsShape, MatrixShape, ProductTag> + : product_accumulate_through_temporary< + Transpose<Lhs>, Rhs, + generic_product_impl<Transpose<Lhs>, Rhs, TranspositionsShape, MatrixShape, ProductTag>> { template <typename Dest> static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void evalTo(Dest& dst, const Transpose<Lhs>& lhs, const Rhs& rhs) { transposition_matrix_product<Rhs, OnTheLeft, true, MatrixShape>::run(dst, lhs.nestedExpression(), rhs); @@ -2103,7 +2150,10 @@ }; template <typename Lhs, typename Rhs, int ProductTag, typename MatrixShape> -struct generic_product_impl<Lhs, Transpose<Rhs>, MatrixShape, TranspositionsShape, ProductTag> { +struct generic_product_impl<Lhs, Transpose<Rhs>, MatrixShape, TranspositionsShape, ProductTag> + : product_accumulate_through_temporary< + Lhs, Transpose<Rhs>, + generic_product_impl<Lhs, Transpose<Rhs>, MatrixShape, TranspositionsShape, ProductTag>> { template <typename Dest> static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void evalTo(Dest& dst, const Lhs& lhs, const Transpose<Rhs>& rhs) { transposition_matrix_product<Lhs, OnTheRight, true, MatrixShape>::run(dst, rhs.nestedExpression(), lhs);
diff --git a/Eigen/src/Core/SolveTriangular.h b/Eigen/src/Core/SolveTriangular.h index 04553e9..25e7be5 100644 --- a/Eigen/src/Core/SolveTriangular.h +++ b/Eigen/src/Core/SolveTriangular.h
@@ -164,6 +164,24 @@ } }; +// What a triangular solve reads its matrix from: the view itself when it must be evaluated, which reads only the +// triangle, or else its nested expression. +template <bool EvalLhs> +struct triangular_solve_lhs_source { + template <typename TriangularType> + EIGEN_DEVICE_FUNC static const TriangularType& run(const TriangularType& tri) { + return tri; + } +}; + +template <> +struct triangular_solve_lhs_source<false> { + template <typename TriangularType> + EIGEN_DEVICE_FUNC static const typename TriangularType::NestedExpression& run(const TriangularType& tri) { + return tri.nestedExpression(); + } +}; + } // end namespace internal /*************************************************************************** @@ -193,9 +211,20 @@ typename internal::plain_matrix_type<OtherDerived>::type>; using OtherCopy = std::conditional_t<copy, OtherPlainObject, OtherDerived&>; OtherCopy otherCopy(other); + using OtherCopyCleaned = std::remove_reference_t<OtherCopy>; - internal::triangular_solver_selector<MatrixType, std::remove_reference_t<OtherCopy>, Side, Mode>::run( - derived().nestedExpression(), otherCopy); + // The kernels read the matrix blas_traits extracts, without its scalar factor. A unit diagonal does not + // scale with that factor, so a scaled matrix is evaluated rather than its factor applied to the solution. + // Only the triangle is evaluated: the other coefficients, and an implicit unit diagonal, need not be readable. + // The unrolled solver reads the scaled coefficients itself and needs no evaluated copy. + constexpr bool EvalLhs = + internal::blas_traits<MatrixType>::HasScalarFactor && + int(internal::trsolve_traits<MatrixType, OtherCopyCleaned, Side>::Unrolling) != int(CompleteUnrolling); + using Lhs = std::conditional_t<EvalLhs, typename MatrixType::PlainObject, MatrixType>; + std::conditional_t<EvalLhs, const Lhs, const MatrixType&> lhs( + internal::triangular_solve_lhs_source<EvalLhs>::run(derived())); + + internal::triangular_solver_selector<Lhs, OtherCopyCleaned, Side, Mode>::run(lhs, otherCopy); if (copy) other = otherCopy; }
diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h index a7ad5f2..2aea7f0 100644 --- a/Eigen/src/Core/TriangularMatrix.h +++ b/Eigen/src/Core/TriangularMatrix.h
@@ -480,10 +480,11 @@ /** \returns the determinant of the triangular matrix * \sa MatrixBase::determinant() */ EIGEN_DEVICE_FUNC Scalar determinant() const { - EIGEN_IF_CONSTEXPR (Mode & UnitDiag) { + EIGEN_IF_CONSTEXPR (int(Mode) & int(UnitDiag)) { return 1; - } else EIGEN_IF_CONSTEXPR (Mode & ZeroDiag) { - return 0; + } else EIGEN_IF_CONSTEXPR (int(Mode) & int(ZeroDiag)) { + // A zero diagonal makes the product vanish unless it is empty. + return m_matrix.rows() == 0 ? 1 : 0; } else { return m_matrix.diagonal().prod(); }
diff --git a/Eigen/src/Core/products/GeneralMatrixMatrixTriangular_BLAS.h b/Eigen/src/Core/products/GeneralMatrixMatrixTriangular_BLAS.h index 9f24ace..e7f824c 100644 --- a/Eigen/src/Core/products/GeneralMatrixMatrixTriangular_BLAS.h +++ b/Eigen/src/Core/products/GeneralMatrixMatrixTriangular_BLAS.h
@@ -54,8 +54,9 @@ static EIGEN_STRONG_INLINE void run(Index size, Index depth, const Scalar* lhs, Index lhsStride, \ const Scalar* rhs, Index rhsStride, Scalar* res, Index resIncr, \ Index resStride, Scalar alpha, level3_blocking<Scalar, Scalar>& blocking) { \ - EIGEN_IF_CONSTEXPR ((UpLo & (Lower | Upper)) == UpLo) { \ - if (lhs == rhs) { \ + /* SYRK computes lhs * lhs^T: rhs must be the same storage read in the opposite order. */ \ + EIGEN_IF_CONSTEXPR ((UpLo & (Lower | Upper)) == UpLo && LhsStorageOrder != RhsStorageOrder) { \ + if (lhs == rhs && lhsStride == rhsStride) { \ general_matrix_matrix_rankupdate<Index, Scalar, LhsStorageOrder, ConjugateLhs, ColMajor, UpLo>::run( \ size, depth, lhs, lhsStride, rhs, rhsStride, res, resStride, alpha, blocking); \ return; \
diff --git a/Eigen/src/Core/products/GeneralMatrixVector_BLAS.h b/Eigen/src/Core/products/GeneralMatrixVector_BLAS.h index a36006d..23b0695 100644 --- a/Eigen/src/Core/products/GeneralMatrixVector_BLAS.h +++ b/Eigen/src/Core/products/GeneralMatrixVector_BLAS.h
@@ -114,10 +114,12 @@ x_ptr = x_tmp.data(); \ incx = 1; \ } else { \ - x_ptr = rhs; \ + /* BLAS reads a negative increment from the lowest address, Eigen from entry 0. */ \ + x_ptr = rhsIncr < 0 ? rhs + (cols - 1) * rhsIncr : rhs; \ } \ + auto* y_ptr = resIncr < 0 ? res + (rows - 1) * resIncr : res; \ BLASFUNC(&trans, &m, &n, (const BLASTYPE*)&numext::real_ref(alpha), (const BLASTYPE*)lhs, &lda, \ - (const BLASTYPE*)x_ptr, &incx, (const BLASTYPE*)&numext::real_ref(beta), (BLASTYPE*)res, &incy); \ + (const BLASTYPE*)x_ptr, &incx, (const BLASTYPE*)&numext::real_ref(beta), (BLASTYPE*)y_ptr, &incy); \ } \ };
diff --git a/Eigen/src/Core/products/SelfadjointProduct.h b/Eigen/src/Core/products/SelfadjointProduct.h index d85a40b..a58925e 100644 --- a/Eigen/src/Core/products/SelfadjointProduct.h +++ b/Eigen/src/Core/products/SelfadjointProduct.h
@@ -153,7 +153,9 @@ using ActualOtherType_ = internal::remove_all_t<ActualOtherType>; internal::add_const_on_value_type_t<ActualOtherType> actualOther = OtherBlasTraits::extract(other.derived()); - Scalar actualAlpha = alpha * OtherBlasTraits::extractScalarFactor(other.derived()); + // (s w) (s w)^* = s conj(s) w w^* + const Scalar otherFactor = OtherBlasTraits::extractScalarFactor(other.derived()); + Scalar actualAlpha = alpha * otherFactor * numext::conj(otherFactor); enum { StorageOrder = (internal::traits<MatrixType>::Flags & RowMajorBit) ? RowMajor : ColMajor, @@ -188,7 +190,9 @@ using ActualOtherType_ = internal::remove_all_t<ActualOtherType>; internal::add_const_on_value_type_t<ActualOtherType> actualOther = OtherBlasTraits::extract(other.derived()); - Scalar actualAlpha = alpha * OtherBlasTraits::extractScalarFactor(other.derived()); + // (s W) (s W)^* = s conj(s) W W^* + const Scalar otherFactor = OtherBlasTraits::extractScalarFactor(other.derived()); + Scalar actualAlpha = alpha * otherFactor * numext::conj(otherFactor); enum { IsRowMajor = (internal::traits<MatrixType>::Flags & RowMajorBit) ? 1 : 0,
diff --git a/Eigen/src/Core/products/SelfadjointRank2Update.h b/Eigen/src/Core/products/SelfadjointRank2Update.h index 6d9b217..322eb64 100644 --- a/Eigen/src/Core/products/SelfadjointRank2Update.h +++ b/Eigen/src/Core/products/SelfadjointRank2Update.h
@@ -195,6 +195,12 @@ template <typename DerivedU, typename DerivedV> EIGEN_DEVICE_FUNC SelfAdjointView<MatrixType, UpLo>& SelfAdjointView<MatrixType, UpLo>::rankUpdate( const MatrixBase<DerivedU>& u, const MatrixBase<DerivedV>& v, const Scalar& alpha) { + // The rank-2 kernel below addresses the matrix with a unit inner stride. + if (MatrixType::InnerStrideAtCompileTime != 1 && nestedExpression().innerStride() != 1) { + nestedExpression().template triangularView<UpLo>() += + typename MatrixType::PlainObject(alpha * u * v.adjoint() + numext::conj(alpha) * v * u.adjoint()); + return *this; + } using UBlasTraits = internal::blas_traits<DerivedU>; using ActualUType = typename UBlasTraits::DirectLinearAccessType; using ActualUType_ = internal::remove_all_t<ActualUType>;
diff --git a/Eigen/src/Core/products/TriangularMatrixMatrix.h b/Eigen/src/Core/products/TriangularMatrixMatrix.h index 6adda5c..a8a4e49 100644 --- a/Eigen/src/Core/products/TriangularMatrixMatrix.h +++ b/Eigen/src/Core/products/TriangularMatrixMatrix.h
@@ -357,17 +357,18 @@ &dst.coeffRef(0, 0), dst.innerStride(), dst.outerStride(), // result info actualAlpha, blocking); - // Apply correction if the diagonal is unit and a scalar factor was nested: + // A factor s nested in a unit triangle reached the kernel through actualAlpha, so the unit diagonal + // contributed alpha * s * other instead of alpha * other: subtract alpha * (s - 1) * other. EIGEN_IF_CONSTEXPR ((Mode & UnitDiag) == UnitDiag) { EIGEN_IF_CONSTEXPR (LhsIsTriangular) { if (!numext::is_exactly_one(lhs_alpha)) { Index diagSize = (std::min)(lhs.rows(), lhs.cols()); - dst.topRows(diagSize) -= ((lhs_alpha - LhsScalar(1)) * a_rhs).topRows(diagSize); + dst.topRows(diagSize) -= ((alpha * (lhs_alpha - LhsScalar(1))) * a_rhs).topRows(diagSize); } } else { if (!numext::is_exactly_one(rhs_alpha)) { Index diagSize = (std::min)(rhs.rows(), rhs.cols()); - dst.leftCols(diagSize) -= (rhs_alpha - RhsScalar(1)) * a_lhs.leftCols(diagSize); + dst.leftCols(diagSize) -= (alpha * (rhs_alpha - RhsScalar(1))) * a_lhs.leftCols(diagSize); } } }
diff --git a/Eigen/src/Core/products/TriangularMatrixVector.h b/Eigen/src/Core/products/TriangularMatrixVector.h index bd685ea..70d571c 100644 --- a/Eigen/src/Core/products/TriangularMatrixVector.h +++ b/Eigen/src/Core/products/TriangularMatrixVector.h
@@ -255,13 +255,16 @@ namespace internal { +// A factor s nested in a unit triangle reached the kernel through the product's alpha, so the unit diagonal +// contributed alpha * s * rhs instead of alpha * rhs: subtract alpha * (s - 1) * rhs. template <int Mode, typename Lhs, typename Rhs, typename Dest, typename LhsScalar> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void trmv_correct_unit_diagonal(const Lhs& lhs, const Rhs& rhs, Dest& dest, + const typename Dest::Scalar& alpha, const LhsScalar& lhs_alpha) { EIGEN_IF_CONSTEXPR ((Mode & UnitDiag) == UnitDiag) { if (!numext::is_exactly_one(lhs_alpha)) { Index diagSize = (std::min)(lhs.rows(), lhs.cols()); - dest.head(diagSize) -= (lhs_alpha - LhsScalar(1)) * rhs.head(diagSize); + dest.head(diagSize) -= (alpha * (lhs_alpha - LhsScalar(1))) * rhs.head(diagSize); } } } @@ -310,7 +313,7 @@ destPolicy.copy_back(dest, actualDestPtr); - trmv_correct_unit_diagonal<Mode>(lhs, rhs, dest, lhs_alpha); + trmv_correct_unit_diagonal<Mode>(lhs, rhs, dest, alpha, lhs_alpha); } }; @@ -357,7 +360,7 @@ dest.innerStride(), actualAlpha); - trmv_correct_unit_diagonal<Mode>(lhs, rhs, dest, lhs_alpha); + trmv_correct_unit_diagonal<Mode>(lhs, rhs, dest, alpha, lhs_alpha); } };
diff --git a/Eigen/src/Core/products/TriangularMatrixVector_BLAS.h b/Eigen/src/Core/products/TriangularMatrixVector_BLAS.h index 9ba3dea..f6e3a8f 100644 --- a/Eigen/src/Core/products/TriangularMatrixVector_BLAS.h +++ b/Eigen/src/Core/products/TriangularMatrixVector_BLAS.h
@@ -133,8 +133,10 @@ (&uplo, &trans, &diag, &n, (const BLASTYPE*)lhs_, &lda, (BLASTYPE*)x, &incx); \ \ /* Add op(a_tr)rhs into res*/ \ + /* BLAS addresses a vector with a negative increment from its lowest address, Eigen from entry 0. */ \ EIGEN_CAT(EIGEN_CAT(BLASPREFIX, axpy), BLASPOSTFIX) \ - (&n, (const BLASTYPE*)&numext::real_ref(alpha), (const BLASTYPE*)x, &incx, (BLASTYPE*)res_, &incy); \ + (&n, (const BLASTYPE*)&numext::real_ref(alpha), (const BLASTYPE*)x, &incx, \ + (BLASTYPE*)(resIncr < 0 ? res_ + (size - 1) * resIncr : res_), &incy); \ /* Non-square case - doesn't fit to BLAS ?TRMV. Fall to default triangular product*/ \ if (size < (std::max)(rows, cols)) { \ EIGEN_IF_CONSTEXPR (ConjRhs) \ @@ -154,6 +156,7 @@ m = convert_index<BlasIndex>(size); \ n = convert_index<BlasIndex>(cols - size); \ } \ + if (resIncr < 0) y += (m - 1) * resIncr; \ EIGEN_CAT(EIGEN_CAT(BLASPREFIX, gemv), BLASPOSTFIX) \ (&trans, &m, &n, (const BLASTYPE*)&numext::real_ref(alpha), (const BLASTYPE*)a, &lda, (const BLASTYPE*)x, \ &incx, (const BLASTYPE*)&numext::real_ref(beta), (BLASTYPE*)y, &incy); \ @@ -231,8 +234,10 @@ (&uplo, &trans, &diag, &n, (const BLASTYPE*)lhs_, &lda, (BLASTYPE*)x, &incx); \ \ /* Add op(a_tr)rhs into res*/ \ + /* BLAS addresses a vector with a negative increment from its lowest address, Eigen from entry 0. */ \ EIGEN_CAT(EIGEN_CAT(BLASPREFIX, axpy), BLASPOSTFIX) \ - (&n, (const BLASTYPE*)&numext::real_ref(alpha), (const BLASTYPE*)x, &incx, (BLASTYPE*)res_, &incy); \ + (&n, (const BLASTYPE*)&numext::real_ref(alpha), (const BLASTYPE*)x, &incx, \ + (BLASTYPE*)(resIncr < 0 ? res_ + (size - 1) * resIncr : res_), &incy); \ /* Non-square case - doesn't fit to BLAS ?TRMV. Fall to default triangular product*/ \ if (size < (std::max)(rows, cols)) { \ EIGEN_IF_CONSTEXPR (ConjRhs) \ @@ -252,6 +257,7 @@ m = convert_index<BlasIndex>(size); \ n = convert_index<BlasIndex>(cols - size); \ } \ + if (resIncr < 0) y += (m - 1) * resIncr; \ EIGEN_CAT(EIGEN_CAT(BLASPREFIX, gemv), BLASPOSTFIX) \ (&trans, &n, &m, (const BLASTYPE*)&numext::real_ref(alpha), (const BLASTYPE*)a, &lda, (const BLASTYPE*)x, \ &incx, (const BLASTYPE*)&numext::real_ref(beta), (BLASTYPE*)y, &incy); \
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 06a7a86..874d5b8 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt
@@ -539,6 +539,7 @@ ei_add_test(num_dimensions) ei_add_test(stl_iterators) ei_add_test(blasutil) +ei_add_test(aocl_vml) ei_add_test(random_matrix) ei_add_test(initializer_list_construction) ei_add_test(diagonal_matrix_variadic_ctor)
diff --git a/test/aocl_vml.cpp b/test/aocl_vml.cpp new file mode 100644 index 0000000..8a12899 --- /dev/null +++ b/test/aocl_vml.cpp
@@ -0,0 +1,92 @@ +// SPDX-FileCopyrightText: The Eigen Authors +// SPDX-License-Identifier: MPL-2.0 + +// Tests the dispatch of coefficient-wise assignments to the AOCL vector math library (Assign_AOCL.h). Without +// AOCL, its entry points are stubbed with their documented semantics, dst[i] = f(src[i]), so that what is tested +// is which assignments Eigen dispatches and how. +#if defined(__has_include) +#if !__has_include("amdlibm.h") +#define EIGEN_TEST_AOCL_STUBS +#endif +#endif + +#ifdef EIGEN_TEST_AOCL_STUBS +#include <cmath> +#define EIGEN_TEST_AOCL_UNARY_STUB(NAME, FUNC) \ + static void NAME(int n, double* x, double* y) { \ + for (int i = 0; i < n; ++i) y[i] = FUNC(x[i]); \ + } +#define EIGEN_TEST_AOCL_BINARY_STUB(NAME, FUNC) \ + static void NAME(int n, double* x, double* z, double* y) { \ + for (int i = 0; i < n; ++i) y[i] = FUNC(x[i], z[i]); \ + } +static double aocl_stub_add(double a, double b) { return a + b; } +EIGEN_TEST_AOCL_UNARY_STUB(amd_vrda_exp2, std::exp2) +EIGEN_TEST_AOCL_UNARY_STUB(amd_vrda_exp, std::exp) +EIGEN_TEST_AOCL_UNARY_STUB(amd_vrda_sin, std::sin) +EIGEN_TEST_AOCL_UNARY_STUB(amd_vrda_cos, std::cos) +EIGEN_TEST_AOCL_UNARY_STUB(amd_vrda_sqrt, std::sqrt) +EIGEN_TEST_AOCL_UNARY_STUB(amd_vrda_cbrt, std::cbrt) +EIGEN_TEST_AOCL_UNARY_STUB(amd_vrda_fabs, std::fabs) +EIGEN_TEST_AOCL_UNARY_STUB(amd_vrda_log, std::log) +EIGEN_TEST_AOCL_UNARY_STUB(amd_vrda_log10, std::log10) +EIGEN_TEST_AOCL_UNARY_STUB(amd_vrda_log2, std::log2) +EIGEN_TEST_AOCL_BINARY_STUB(amd_vrda_add, aocl_stub_add) +EIGEN_TEST_AOCL_BINARY_STUB(amd_vrda_pow, std::pow) +EIGEN_TEST_AOCL_BINARY_STUB(amd_vrda_fmax, std::fmax) +EIGEN_TEST_AOCL_BINARY_STUB(amd_vrda_fmin, std::fmin) +#endif + +#define EIGEN_USE_AOCL_ALL +#include "main.h" + +// Assignment to a plain object resizes it, as without AOCL. +template <int> +void aocl_vml_resize() { + const Index n = 2 * EIGEN_AOCL_VML_THRESHOLD + 1; + const VectorXd x = VectorXd::Random(n), y = VectorXd::Random(n); + VectorXd r; + r = x.array().exp(); + VERIFY_IS_EQUAL(r.size(), n); + VERIFY_IS_APPROX(r, x.unaryExpr([](double v) { return std::exp(v); })); + r.resize(0); + r = x + y; + VERIFY_IS_EQUAL(r.size(), n); + VERIFY_IS_APPROX(r, x.binaryExpr(y, [](double a, double b) { return a + b; })); +} + +// Operands that are not one contiguous array: blocks of a larger matrix, strided maps, and block destinations, +// whose coefficients outside the block must stay untouched. +template <int> +void aocl_vml_non_contiguous() { + const Index n = 2 * EIGEN_AOCL_VML_THRESHOLD + 3, m = n - 7; + const MatrixXd a = MatrixXd::Random(n, n), b = MatrixXd::Random(n, n); + const MatrixXd aBlock = a.block(1, 2, m, m), bBlock = b.block(3, 0, m, m); + const auto exp = [](double v) { return std::exp(v); }; + const auto sum = [](double x, double y) { return x + y; }; + + MatrixXd r = a.block(1, 2, m, m).array().exp(); + VERIFY_IS_APPROX(r, aBlock.unaryExpr(exp)); + r = a.block(1, 2, m, m) + b.block(3, 0, m, m); + VERIFY_IS_APPROX(r, aBlock.binaryExpr(bBlock, sum)); + + MatrixXd d = MatrixXd::Zero(n, n); + d.block(2, 1, m, m).array() = aBlock.array().exp(); + VERIFY_IS_APPROX(MatrixXd(d.block(2, 1, m, m)), aBlock.unaryExpr(exp)); + d.block(2, 1, m, m).setZero(); + VERIFY(d.isZero()); + d.block(2, 1, m, m) = aBlock + bBlock; + VERIFY_IS_APPROX(MatrixXd(d.block(2, 1, m, m)), aBlock.binaryExpr(bBlock, sum)); + d.block(2, 1, m, m).setZero(); + VERIFY(d.isZero()); + + const Map<const VectorXd, 0, InnerStride<2> > strided(a.data(), n); + const VectorXd stridedCopy = strided; + VectorXd v = strided.array().exp(); + VERIFY_IS_APPROX(v, stridedCopy.unaryExpr(exp)); +} + +EIGEN_DECLARE_TEST(aocl_vml) { + CALL_SUBTEST_1(aocl_vml_resize<0>()); + CALL_SUBTEST_2(aocl_vml_non_contiguous<0>()); +}
diff --git a/test/bandmatrix.cpp b/test/bandmatrix.cpp index 79d87bf..2e2f820 100644 --- a/test/bandmatrix.cpp +++ b/test/bandmatrix.cpp
@@ -56,7 +56,39 @@ using Eigen::internal::BandMatrix; +// diagonal<N>() of a const band matrix, and the sub and super diagonals of a const tridiagonal one. +template <int> +void bandmatrix_const_diagonals() { + BandMatrix<double, Dynamic, Dynamic, 2, 1> band(6, 7, 2, 1); + band.coeffs().setRandom(); + const BandMatrix<double, Dynamic, Dynamic, 2, 1>& constBand = band; + VERIFY_IS_EQUAL(RowVectorXd(constBand.diagonal<2>()), RowVectorXd(band.diagonal(2))); + VERIFY_IS_EQUAL(RowVectorXd(constBand.diagonal<-1>()), RowVectorXd(band.diagonal(-1))); + + Eigen::internal::TridiagonalMatrix<double, Dynamic, 0> tri(5); + tri.coeffs().setRandom(); + const Eigen::internal::TridiagonalMatrix<double, Dynamic, 0>& constTri = tri; + VERIFY_IS_EQUAL(RowVectorXd(constTri.super()), RowVectorXd(tri.diagonal(1))); + VERIFY_IS_EQUAL(RowVectorXd(constTri.sub()), RowVectorXd(tri.diagonal(-1))); +} + +// A self-adjoint tridiagonal matrix stores the subdiagonal only; its superdiagonal is the conjugate. +template <typename Scalar> +void bandmatrix_selfadjoint_diagonals() { + using RowVector = Matrix<Scalar, 1, Dynamic>; + Eigen::internal::TridiagonalMatrix<Scalar, Dynamic, SelfAdjoint> tri(5); + tri.coeffs().setRandom(); + const Eigen::internal::TridiagonalMatrix<Scalar, Dynamic, SelfAdjoint>& constTri = tri; + const RowVector sub = tri.diagonal(-1); + VERIFY_IS_EQUAL(RowVector(tri.sub()), sub); + VERIFY_IS_EQUAL(RowVector(tri.super()), RowVector(sub.conjugate())); + VERIFY_IS_EQUAL(RowVector(constTri.super()), RowVector(sub.conjugate())); +} + EIGEN_DECLARE_TEST(bandmatrix) { + CALL_SUBTEST(bandmatrix_const_diagonals<0>()); + CALL_SUBTEST(bandmatrix_selfadjoint_diagonals<double>()); + CALL_SUBTEST(bandmatrix_selfadjoint_diagonals<std::complex<double>>()); for (int i = 0; i < 10 * g_repeat; i++) { Index rows = internal::random<Index>(1, 10); Index cols = internal::random<Index>(1, 10);
diff --git a/test/determinant.cpp b/test/determinant.cpp index 741a7f4..4bafa01 100644 --- a/test/determinant.cpp +++ b/test/determinant.cpp
@@ -155,6 +155,14 @@ // check empty matrix VERIFY_IS_APPROX(m2.block(0, 0, 0, 0).determinant(), Scalar(1)); + + // triangular views + VERIFY_IS_APPROX(m2.template triangularView<Upper>().determinant(), m2.diagonal().prod()); + VERIFY_IS_EQUAL(m2.template triangularView<UnitLower>().determinant(), Scalar(1)); + VERIFY_IS_EQUAL(m2.template triangularView<StrictlyUpper>().determinant(), Scalar(0)); + VERIFY_IS_EQUAL(m2.block(0, 0, 0, 0).template triangularView<StrictlyLower>().determinant(), Scalar(1)); + VERIFY_IS_EQUAL(m2.block(0, 0, 0, 0).template triangularView<UnitUpper>().determinant(), Scalar(1)); + VERIFY_IS_EQUAL(m2.block(0, 0, 0, 0).template triangularView<Lower>().determinant(), Scalar(1)); } EIGEN_DECLARE_TEST(determinant) {
diff --git a/test/permutationmatrices.cpp b/test/permutationmatrices.cpp index 743daeb..d31ec48 100644 --- a/test/permutationmatrices.cpp +++ b/test/permutationmatrices.cpp
@@ -167,6 +167,26 @@ VERIFY_IS_APPROX(m_original * rt.adjoint(), m_original * rp); m_permuted = m_original; VERIFY_IS_APPROX(m_permuted = lt.adjoint() * m_permuted, lp.inverse() * m_original); + + // Accumulating permutation and transposition products without aliasing. + const MatrixType m_base = MatrixType::Random(rows, cols); + lm = lp; + rm = rp; + m_permuted = m_base; + m_permuted.noalias() += lp * m_original; + VERIFY_IS_APPROX(m_permuted, m_base + lm * m_original); + m_permuted = m_base; + m_permuted.noalias() -= m_original * rp.inverse(); + VERIFY_IS_APPROX(m_permuted, m_base - m_original * rm.transpose()); + m_permuted = m_base; + m_permuted.noalias() += Scalar(2) * (m_original * rp); + VERIFY_IS_APPROX(m_permuted, m_base + Scalar(2) * (m_original * rm)); + m_permuted = m_base; + m_permuted.noalias() -= lt * m_original; + VERIFY_IS_APPROX(m_permuted, m_base - lm * m_original); + m_permuted = m_base; + m_permuted.noalias() += m_original * rt.transpose(); + VERIFY_IS_APPROX(m_permuted, m_base + m_original * rm); } template <typename T>
diff --git a/test/product_extra.cpp b/test/product_extra.cpp index ac1f22e..00ddc52 100644 --- a/test/product_extra.cpp +++ b/test/product_extra.cpp
@@ -98,6 +98,12 @@ VERIFY_IS_APPROX(m1 * mapvc2, m1 * vc2.reverse().eval()); VERIFY_IS_APPROX(map1.adjoint() * v1.transpose(), m1.adjoint().reverse() * v1.transpose()); VERIFY_IS_APPROX(m1.adjoint() * mapv1.transpose(), m1.adjoint() * v1.reverse().transpose()); + + // a destination with a negative inner stride + ColVectorType reversedRes(cols); + Map<ColVectorType, Unaligned, InnerStride<-1> > mapvcres(&reversedRes(cols - 1), cols, InnerStride<-1>(-1)); + mapvcres.noalias() = m1.adjoint() * v1.transpose(); + VERIFY_IS_APPROX(reversedRes, (m1.adjoint() * v1.transpose()).eval().reverse()); } // regression test
diff --git a/test/product_selfadjoint.cpp b/test/product_selfadjoint.cpp index e017e67..a7fd13e 100644 --- a/test/product_selfadjoint.cpp +++ b/test/product_selfadjoint.cpp
@@ -156,6 +156,22 @@ VERIFY_IS_APPROX(m2, m3.template triangularView<Lower>().toDenseMatrix()); } + // rank2 update of a matrix with a non-unit inner stride + { + using VectorX = Matrix<Scalar, Dynamic, 1>; + const Index innerSize = MatrixType::IsRowMajor ? cols : rows; + VectorX buffer = VectorX::Zero(2 * rows * cols); + Map<MatrixType, 0, Stride<Dynamic, 2> > map1(buffer.data(), rows, cols, Stride<Dynamic, 2>(2 * innerSize, 2)); + map1 = m1.template triangularView<Upper>(); + map1.template selfadjointView<Upper>().rankUpdate(v1, v2, s1); + m3 = map1; + VERIFY_IS_APPROX(m3, (m1 + s1 * v1 * v2.adjoint() + numext::conj(s1) * v2 * v1.adjoint()) + .template triangularView<Upper>() + .toDenseMatrix()); + const Map<VectorX, 0, InnerStride<2> > gaps(buffer.data() + 1, rows * cols); + VERIFY(gaps.isZero()); + } + // matrix-vector m2 = m1.template triangularView<Lower>(); VERIFY_IS_APPROX(m1 * m4, m2.template selfadjointView<Lower>() * m4);
diff --git a/test/product_syrk.cpp b/test/product_syrk.cpp index 1035e86..0ce1bf8 100644 --- a/test/product_syrk.cpp +++ b/test/product_syrk.cpp
@@ -141,6 +141,34 @@ .template triangularView<Upper>() .toDenseMatrix())); + // The product of a matrix with itself, rather than with its adjoint, shares its storage without being a rank update. + m2.setZero(); + m2.template triangularView<Lower>() = m1 * m1; + VERIFY_IS_APPROX(m2, (m1 * m1).eval().template triangularView<Lower>().toDenseMatrix()); + rm2.setZero(); + rm2.template triangularView<Upper>() += s1 * m1 * m1; + VERIFY_IS_APPROX(rm2, (s1 * m1 * m1).eval().template triangularView<Upper>().toDenseMatrix()); + + // A scalar factor of the update enters twice: (s u) (s u)^* = s conj(s) u u^*. + { + const Scalar s2 = internal::random<Scalar>(); + const Rhs2 scaled2 = s2 * rhs2; + using VectorX = Matrix<Scalar, Dynamic, 1>; + const VectorX scaledCol = -(s2 * m1.col(c)), conjCol = (s2 * m1.col(c)).conjugate(); + m2.setZero(); + VERIFY_IS_APPROX((m2.template selfadjointView<Lower>().rankUpdate(s2 * rhs2, s1)._expression()), + ((s1 * scaled2 * scaled2.adjoint()).eval().template triangularView<Lower>().toDenseMatrix())); + rm2.setZero(); + VERIFY_IS_APPROX((rm2.template selfadjointView<Upper>().rankUpdate(rhs2 * s2, s1)._expression()), + ((s1 * scaled2 * scaled2.adjoint()).eval().template triangularView<Upper>().toDenseMatrix())); + m2.setZero(); + VERIFY_IS_APPROX((m2.template selfadjointView<Upper>().rankUpdate(-(s2 * m1.col(c)), s1)._expression()), + ((s1 * scaledCol * scaledCol.adjoint()).eval().template triangularView<Upper>().toDenseMatrix())); + rm2.setZero(); + VERIFY_IS_APPROX((rm2.template selfadjointView<Lower>().rankUpdate((s2 * m1.col(c)).conjugate(), s1)._expression()), + ((s1 * conjCol * conjCol.adjoint()).eval().template triangularView<Lower>().toDenseMatrix())); + } + // destination with a non-default inner-stride // see bug 1741 {
diff --git a/test/product_trmm.cpp b/test/product_trmm.cpp index 6d53745..0079356 100644 --- a/test/product_trmm.cpp +++ b/test/product_trmm.cpp
@@ -88,6 +88,17 @@ VERIFY_IS_APPROX(ge_xs = (s1 * mat).transpose().template triangularView<Mode>() * ge_left.adjoint(), s1triTr * ge_left.adjoint()); + // Subtracting a product whose triangle has a folded scalar factor, which a unit diagonal does not share. + ge_xs.setRandom(); + ge_xs_save = ge_xs; + VERIFY_IS_APPROX( + (ge_xs_save - s1triTr * (s2 * ge_left.adjoint())).eval(), + ge_xs.noalias() -= (s1 * mat.transpose()).template triangularView<Mode>() * (s2 * ge_left.adjoint())); + ge_sx.setRandom(ge_left.rows(), mat.cols()); + ge_sx_save = ge_sx; + VERIFY_IS_APPROX((ge_sx_save - (s2 * ge_left) * s1tri).eval(), + ge_sx.noalias() -= (s2 * ge_left) * (s1 * mat).template triangularView<Mode>()); + // TODO check with sub-matrix expressions ? // Triangular times triangular or self-adjoint view: the right factor is densified and the product
diff --git a/test/product_trmv.cpp b/test/product_trmv.cpp index f912cfb..067df4b 100644 --- a/test/product_trmv.cpp +++ b/test/product_trmv.cpp
@@ -79,6 +79,44 @@ // TODO check with sub-matrices } +// A destination with a negative inner stride, which BLAS addresses from its lowest entry and Eigen from entry 0. +// Square and both non-square shapes reach the TRMV, AXPY and GEMV calls of the BLAS backend. +template <typename Scalar, int StorageOrder, int Mode> +void trmv_negative_stride_destination(Index rows, Index cols) { + using Mat = Matrix<Scalar, Dynamic, Dynamic, StorageOrder>; + using Vec = Matrix<Scalar, Dynamic, 1>; + const Mat m = Mat::Random(rows, cols); + const Vec v = Vec::Random(cols); + const Mat dense = m.template triangularView<Mode>(); + const Vec ref = dense * v; + + // The reversed destination sits between two guard blocks of the same length. + Vec buffer = Vec::Random(3 * rows); + const Vec initial = buffer; + Map<Vec, Unaligned, InnerStride<-1>> dst(buffer.data() + 2 * rows - 1, rows, InnerStride<-1>(-1)); + dst.noalias() = m.template triangularView<Mode>() * v; + VERIFY_IS_APPROX(Vec(dst), ref); + VERIFY_IS_EQUAL(buffer.head(rows), initial.head(rows)); + VERIFY_IS_EQUAL(buffer.tail(rows), initial.tail(rows)); +} + +template <typename Scalar, int StorageOrder> +void trmv_negative_stride_destination_shapes() { + const Index shapes[][2] = {{7, 7}, {9, 5}, {5, 9}}; + for (const auto& shape : shapes) { + trmv_negative_stride_destination<Scalar, StorageOrder, Lower>(shape[0], shape[1]); + trmv_negative_stride_destination<Scalar, StorageOrder, Upper>(shape[0], shape[1]); + trmv_negative_stride_destination<Scalar, StorageOrder, UnitLower>(shape[0], shape[1]); + trmv_negative_stride_destination<Scalar, StorageOrder, UnitUpper>(shape[0], shape[1]); + } +} + +template <typename Scalar> +void trmv_negative_stride_destination_all() { + trmv_negative_stride_destination_shapes<Scalar, ColMajor>(); + trmv_negative_stride_destination_shapes<Scalar, RowMajor>(); +} + EIGEN_DECLARE_TEST(product_trmv) { int s = 0; for (int i = 0; i < g_repeat; i++) { @@ -95,4 +133,8 @@ CALL_SUBTEST_6(trmv(Matrix<float, Dynamic, Dynamic, RowMajor>(s, s))); TEST_SET_BUT_UNUSED_VARIABLE(s); } + CALL_SUBTEST_7(trmv_negative_stride_destination_all<float>()); + CALL_SUBTEST_7(trmv_negative_stride_destination_all<double>()); + CALL_SUBTEST_7(trmv_negative_stride_destination_all<std::complex<float>>()); + CALL_SUBTEST_7(trmv_negative_stride_destination_all<std::complex<double>>()); }
diff --git a/test/product_trsolve.cpp b/test/product_trsolve.cpp index 5b130bd..145fef7 100644 --- a/test/product_trsolve.cpp +++ b/test/product_trsolve.cpp
@@ -14,6 +14,44 @@ #include "main.h" +// Coefficients of a Mode triangle that flag every read outside the stored triangle, including the diagonal of a +// unit triangle. +template <typename Scalar, int Mode> +struct triangle_only_functor { + bool* touched; + Scalar operator()(Index i, Index j) const { + const bool outside = ((Mode & Upper) ? i > j : i < j) || (i == j && (Mode & UnitDiag)); + if (outside) *touched = true; + return i == j ? Scalar(2) : Scalar(1) / Scalar(Index(2) + i + j); + } +}; + +template <int Mode, typename Scalar> +void trsolve_scaled_triangle_only(Index size) { + using Mat = Matrix<Scalar, Dynamic, Dynamic>; + using Vec = Matrix<Scalar, Dynamic, 1>; + bool touched = false; + const auto a = Mat::NullaryExpr(size, size, triangle_only_functor<Scalar, Mode>{&touched}); + const Mat dense = a.template triangularView<Mode>(); + VERIFY(!touched); + const Scalar s(3); + const Vec b = Vec::Random(size); + Vec x = b; + (s * a).template triangularView<Mode>().solveInPlace(x); + VERIFY(!touched); + Vec ref = b; + (s * dense).eval().template triangularView<Mode>().solveInPlace(ref); + VERIFY_IS_APPROX(x, ref); +} + +template <typename Scalar> +void trsolve_scaled_reads_triangle_only(Index size) { + trsolve_scaled_triangle_only<Lower, Scalar>(size); + trsolve_scaled_triangle_only<Upper, Scalar>(size); + trsolve_scaled_triangle_only<UnitLower, Scalar>(size); + trsolve_scaled_triangle_only<UnitUpper, Scalar>(size); +} + #define VERIFY_TRSM(TRI, XB) \ { \ (XB).setRandom(); \ @@ -81,6 +119,17 @@ VERIFY_TRSM_ONTHERIGHT(rmLhs.template triangularView<Lower>(), cmRhs); VERIFY_TRSM_ONTHERIGHT(rmLhs.conjugate().template triangularView<UnitUpper>(), rmRhs); + // Triangles with a scalar factor, which a unit diagonal does not share. + const Scalar s = Scalar(RealScalar(1.5)) + RealScalar(0.25) * internal::random<Scalar>(); + VERIFY_TRSM((s * cmLhs).template triangularView<Lower>(), cmRhs); + VERIFY_TRSM((-rmLhs).template triangularView<Upper>(), rmRhs); + VERIFY_TRSM(s * cmLhs.template triangularView<Upper>(), rmRhs); + VERIFY_TRSM((cmLhs * s).adjoint().template triangularView<UnitLower>(), cmRhs); + VERIFY_TRSM((s * rmLhs).template triangularView<UnitUpper>(), cmRhs); + VERIFY_TRSM_ONTHERIGHT((s * cmLhs).template triangularView<UnitLower>(), cmRhs); + VERIFY_TRSM_ONTHERIGHT((s * rmLhs).template triangularView<Upper>(), rmRhs); + trsolve_scaled_reads_triangle_only<Scalar>(size); + int c = internal::random<int>(0, cols - 1); VERIFY_TRSM(rmLhs.template triangularView<Lower>(), rmRhs.col(c)); VERIFY_TRSM(cmLhs.template triangularView<Lower>(), rmRhs.col(c)); @@ -341,6 +390,40 @@ VERIFY_TRSM_NO_MALLOC_ONTHERIGHT(lhs.adjoint().template triangularView<Lower>(), xRight, refRight); } +// A right-hand side of fixed size at most 8 takes the unrolled solver, which reads the coefficients of a scaled +// triangle itself: it needs no evaluated copy, and so no allocation, even when the matrix is dynamic. +template <typename Scalar, int Mode> +void trsolve_scaled_unrolled_no_malloc() { + using Mat = Matrix<Scalar, Dynamic, Dynamic>; + using Vec = Matrix<Scalar, 4, 1>; + using RowVec = Matrix<Scalar, 1, 4>; + using RealScalar = typename NumTraits<Scalar>::Real; + STATIC_CHECK((int(internal::trsolve_traits<Mat, Vec, OnTheLeft>::Unrolling) == CompleteUnrolling)); + STATIC_CHECK((int(internal::trsolve_traits<Mat, RowVec, OnTheRight>::Unrolling) == CompleteUnrolling)); + + Mat lhs = Mat::Random(4, 4) * RealScalar(0.1); + lhs.diagonal().array() += RealScalar(1); + const Scalar s(3); + Vec x, ref = Vec::Random(); + VERIFY_TRSM_NO_MALLOC((s * lhs).template triangularView<Mode>(), x, ref); + RowVec xRight, refRight = RowVec::Random(); + VERIFY_TRSM_NO_MALLOC_ONTHERIGHT((s * lhs).template triangularView<Mode>(), xRight, refRight); + + bool touched = false; + const auto a = Mat::NullaryExpr(4, 4, triangle_only_functor<Scalar, Mode>{&touched}); + VERIFY_TRSM_NO_MALLOC((s * a).template triangularView<Mode>(), x, ref); + VERIFY_TRSM_NO_MALLOC_ONTHERIGHT((s * a).template triangularView<Mode>(), xRight, refRight); + VERIFY(!touched); +} + +template <typename Scalar> +void trsolve_scaled_unrolled_no_malloc_all() { + trsolve_scaled_unrolled_no_malloc<Scalar, Lower>(); + trsolve_scaled_unrolled_no_malloc<Scalar, Upper>(); + trsolve_scaled_unrolled_no_malloc<Scalar, UnitLower>(); + trsolve_scaled_unrolled_no_malloc<Scalar, UnitUpper>(); +} + template <int> void trsolve_no_malloc_all() { // The packed solver walks the triangular matrix in panels of max(mr, nr) columns; these sizes @@ -358,6 +441,9 @@ trsolve_no_malloc<std::complex<double>, ColMajor>(size, cols); } } + trsolve_scaled_unrolled_no_malloc_all<float>(); + trsolve_scaled_unrolled_no_malloc_all<double>(); + trsolve_scaled_unrolled_no_malloc_all<std::complex<double>>(); } // The blocked solvers take their k-block depth, and a solve on the left its column panels, from the
diff --git a/test/triangular.cpp b/test/triangular.cpp index 132cbfe..b237c67 100644 --- a/test/triangular.cpp +++ b/test/triangular.cpp
@@ -167,6 +167,19 @@ triangular_permutation_products_check(a.template triangularView<StrictlyUpper>(), perm); triangular_permutation_products_check(a.template selfadjointView<Lower>(), perm); triangular_permutation_products_check(a.template selfadjointView<Upper>(), perm); + + // A complex alpha folded into a self-adjoint view scales the whole view, as in its other products, rather than + // being conjugated with the mirrored triangle. + using Scalar = typename MatrixType::Scalar; + using RealScalar = typename NumTraits<Scalar>::Real; + const Scalar alpha = internal::random<Scalar>() + Scalar(RealScalar(2)); + const MatrixType lower = a.template selfadjointView<Lower>(), upper = a.template selfadjointView<Upper>(); + const MatrixType permDense = perm.toDenseMatrix().template cast<Scalar>(); + MatrixType result = alpha * (a.template selfadjointView<Lower>() * perm); + VERIFY_IS_APPROX(result, MatrixType(alpha * MatrixType(lower * permDense))); + result = a; + result += alpha * (a.template selfadjointView<Upper>() * perm.inverse()); + VERIFY_IS_APPROX(result, MatrixType(a + alpha * MatrixType(upper * permDense.transpose()))); } template <typename MatrixType>