diff options
-rw-r--r-- | Eigen/src/Core/AssignEvaluator.h | 36 | ||||
-rw-r--r-- | Eigen/src/Core/DiagonalMatrix.h | 5 | ||||
-rw-r--r-- | Eigen/src/Core/ProductEvaluators.h | 12 | ||||
-rw-r--r-- | Eigen/src/Core/Solve.h | 16 | ||||
-rw-r--r-- | Eigen/src/Core/Transpose.h | 5 | ||||
-rw-r--r-- | Eigen/src/Core/TriangularMatrix.h | 18 | ||||
-rw-r--r-- | Eigen/src/Geometry/Homogeneous.h | 10 | ||||
-rw-r--r-- | Eigen/src/IterativeLinearSolvers/SolveWithGuess.h | 6 | ||||
-rw-r--r-- | Eigen/src/LU/InverseImpl.h | 6 | ||||
-rw-r--r-- | Eigen/src/SparseCore/SparseAssign.h | 17 | ||||
-rw-r--r-- | Eigen/src/SparseCore/SparseProduct.h | 5 | ||||
-rw-r--r-- | test/product_extra.cpp | 45 |
12 files changed, 152 insertions, 29 deletions
diff --git a/Eigen/src/Core/AssignEvaluator.h b/Eigen/src/Core/AssignEvaluator.h index abad8c790..ffe1dd0ca 100644 --- a/Eigen/src/Core/AssignEvaluator.h +++ b/Eigen/src/Core/AssignEvaluator.h @@ -697,15 +697,21 @@ protected: ***************************************************************************/ template<typename DstXprType, typename SrcXprType, typename Functor> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void call_dense_assignment_loop(const DstXprType& dst, const SrcXprType& src, const Functor &func) +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void call_dense_assignment_loop(DstXprType& dst, const SrcXprType& src, const Functor &func) { - eigen_assert(dst.rows() == src.rows() && dst.cols() == src.cols()); - typedef evaluator<DstXprType> DstEvaluatorType; typedef evaluator<SrcXprType> SrcEvaluatorType; - DstEvaluatorType dstEvaluator(dst); SrcEvaluatorType srcEvaluator(src); + + // NOTE To properly handle A = (A*A.transpose())/s with A rectangular, + // we need to resize the destination after the source evaluator has been created. + Index dstRows = src.rows(); + Index dstCols = src.cols(); + if((dst.rows()!=dstRows) || (dst.cols()!=dstCols)) + dst.resize(dstRows, dstCols); + + DstEvaluatorType dstEvaluator(dst); typedef generic_dense_assignment_kernel<DstEvaluatorType,SrcEvaluatorType,Functor> Kernel; Kernel kernel(dstEvaluator, srcEvaluator, func, dst.const_cast_derived()); @@ -714,7 +720,7 @@ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void call_dense_assignment_loop(const DstX } template<typename DstXprType, typename SrcXprType> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void call_dense_assignment_loop(const DstXprType& dst, const SrcXprType& src) +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void call_dense_assignment_loop(DstXprType& dst, const SrcXprType& src) { call_dense_assignment_loop(dst, src, internal::assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar>()); } @@ -796,11 +802,6 @@ void call_assignment_no_alias(Dst& dst, const Src& src, const Func& func) ) && int(Dst::SizeAtCompileTime) != 1 }; - Index dstRows = NeedToTranspose ? src.cols() : src.rows(); - Index dstCols = NeedToTranspose ? src.rows() : src.cols(); - if((dst.rows()!=dstRows) || (dst.cols()!=dstCols)) - dst.resize(dstRows, dstCols); - typedef typename internal::conditional<NeedToTranspose, Transpose<Dst>, Dst>::type ActualDstTypeCleaned; typedef typename internal::conditional<NeedToTranspose, Transpose<Dst>, Dst&>::type ActualDstType; ActualDstType actualDst(dst); @@ -823,15 +824,11 @@ template<typename Dst, typename Src, typename Func> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void call_assignment_no_alias_no_transpose(Dst& dst, const Src& src, const Func& func) { - Index dstRows = src.rows(); - Index dstCols = src.cols(); - if((dst.rows()!=dstRows) || (dst.cols()!=dstCols)) - dst.resize(dstRows, dstCols); - // TODO check whether this is the right place to perform these checks: EIGEN_STATIC_ASSERT_LVALUE(Dst) EIGEN_STATIC_ASSERT_SAME_MATRIX_SIZE(Dst,Src) - + EIGEN_CHECK_BINARY_COMPATIBILIY(Func,typename Dst::Scalar,typename Src::Scalar); + Assignment<Dst,Src,Func>::run(dst, src, func); } template<typename Dst, typename Src> @@ -853,8 +850,6 @@ struct Assignment<DstXprType, SrcXprType, Functor, Dense2Dense, Weak> EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE void run(DstXprType &dst, const SrcXprType &src, const Functor &func) { - eigen_assert(dst.rows() == src.rows() && dst.cols() == src.cols()); - #ifndef EIGEN_NO_DEBUG internal::check_for_aliasing(dst, src); #endif @@ -873,6 +868,11 @@ struct Assignment<DstXprType, SrcXprType, Functor, EigenBase2EigenBase, Weak> EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar> &/*func*/) { + Index dstRows = src.rows(); + Index dstCols = src.cols(); + if((dst.rows()!=dstRows) || (dst.cols()!=dstCols)) + dst.resize(dstRows, dstCols); + eigen_assert(dst.rows() == src.rows() && dst.cols() == src.cols()); src.evalTo(dst); } diff --git a/Eigen/src/Core/DiagonalMatrix.h b/Eigen/src/Core/DiagonalMatrix.h index c682c6d7f..ecfdce8ef 100644 --- a/Eigen/src/Core/DiagonalMatrix.h +++ b/Eigen/src/Core/DiagonalMatrix.h @@ -320,6 +320,11 @@ struct Assignment<DstXprType, SrcXprType, Functor, Diagonal2Dense> { static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar> &/*func*/) { + Index dstRows = src.rows(); + Index dstCols = src.cols(); + if((dst.rows()!=dstRows) || (dst.cols()!=dstCols)) + dst.resize(dstRows, dstCols); + dst.setZero(); dst.diagonal() = src.diagonal(); } diff --git a/Eigen/src/Core/ProductEvaluators.h b/Eigen/src/Core/ProductEvaluators.h index dfd76990d..c9e2e1a07 100644 --- a/Eigen/src/Core/ProductEvaluators.h +++ b/Eigen/src/Core/ProductEvaluators.h @@ -140,6 +140,10 @@ struct Assignment<DstXprType, Product<Lhs,Rhs,Options>, internal::assign_op<Scal static EIGEN_STRONG_INLINE void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,Scalar> &) { + Index dstRows = src.rows(); + Index dstCols = src.cols(); + if((dst.rows()!=dstRows) || (dst.cols()!=dstCols)) + dst.resize(dstRows, dstCols); // FIXME shall we handle nested_eval here? generic_product_impl<Lhs, Rhs>::evalTo(dst, src.lhs(), src.rhs()); } @@ -154,6 +158,10 @@ struct Assignment<DstXprType, Product<Lhs,Rhs,Options>, internal::add_assign_op< static EIGEN_STRONG_INLINE void run(DstXprType &dst, const SrcXprType &src, const internal::add_assign_op<Scalar,Scalar> &) { + Index dstRows = src.rows(); + Index dstCols = src.cols(); + if((dst.rows()!=dstRows) || (dst.cols()!=dstCols)) + dst.resize(dstRows, dstCols); // FIXME shall we handle nested_eval here? generic_product_impl<Lhs, Rhs>::addTo(dst, src.lhs(), src.rhs()); } @@ -168,6 +176,10 @@ struct Assignment<DstXprType, Product<Lhs,Rhs,Options>, internal::sub_assign_op< static EIGEN_STRONG_INLINE void run(DstXprType &dst, const SrcXprType &src, const internal::sub_assign_op<Scalar,Scalar> &) { + Index dstRows = src.rows(); + Index dstCols = src.cols(); + if((dst.rows()!=dstRows) || (dst.cols()!=dstCols)) + dst.resize(dstRows, dstCols); // FIXME shall we handle nested_eval here? generic_product_impl<Lhs, Rhs>::subTo(dst, src.lhs(), src.rhs()); } diff --git a/Eigen/src/Core/Solve.h b/Eigen/src/Core/Solve.h index 8fc69c4b8..960a58597 100644 --- a/Eigen/src/Core/Solve.h +++ b/Eigen/src/Core/Solve.h @@ -139,7 +139,11 @@ struct Assignment<DstXprType, Solve<DecType,RhsType>, internal::assign_op<Scalar typedef Solve<DecType,RhsType> SrcXprType; static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,Scalar> &) { - // FIXME shall we resize dst here? + Index dstRows = src.rows(); + Index dstCols = src.cols(); + if((dst.rows()!=dstRows) || (dst.cols()!=dstCols)) + dst.resize(dstRows, dstCols); + src.dec()._solve_impl(src.rhs(), dst); } }; @@ -151,6 +155,11 @@ struct Assignment<DstXprType, Solve<Transpose<const DecType>,RhsType>, internal: typedef Solve<Transpose<const DecType>,RhsType> SrcXprType; static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,Scalar> &) { + Index dstRows = src.rows(); + Index dstCols = src.cols(); + if((dst.rows()!=dstRows) || (dst.cols()!=dstCols)) + dst.resize(dstRows, dstCols); + src.dec().nestedExpression().template _solve_impl_transposed<false>(src.rhs(), dst); } }; @@ -163,6 +172,11 @@ struct Assignment<DstXprType, Solve<CwiseUnaryOp<internal::scalar_conjugate_op<t typedef Solve<CwiseUnaryOp<internal::scalar_conjugate_op<typename DecType::Scalar>, const Transpose<const DecType> >,RhsType> SrcXprType; static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,Scalar> &) { + Index dstRows = src.rows(); + Index dstCols = src.cols(); + if((dst.rows()!=dstRows) || (dst.cols()!=dstCols)) + dst.resize(dstRows, dstCols); + src.dec().nestedExpression().nestedExpression().template _solve_impl_transposed<true>(src.rhs(), dst); } }; diff --git a/Eigen/src/Core/Transpose.h b/Eigen/src/Core/Transpose.h index bc232526a..79b767bcc 100644 --- a/Eigen/src/Core/Transpose.h +++ b/Eigen/src/Core/Transpose.h @@ -78,6 +78,11 @@ template<typename MatrixType> class Transpose typename internal::remove_reference<MatrixTypeNested>::type& nestedExpression() { return m_matrix; } + /** \internal */ + void resize(Index nrows, Index ncols) { + m_matrix.resize(ncols,nrows); + } + protected: typename internal::ref_selector<MatrixType>::non_const_type m_matrix; }; diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h index 71f5d4f29..641c20417 100644 --- a/Eigen/src/Core/TriangularMatrix.h +++ b/Eigen/src/Core/TriangularMatrix.h @@ -775,15 +775,18 @@ public: template<int Mode, bool SetOpposite, typename DstXprType, typename SrcXprType, typename Functor> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE -void call_triangular_assignment_loop(const DstXprType& dst, const SrcXprType& src, const Functor &func) +void call_triangular_assignment_loop(DstXprType& dst, const SrcXprType& src, const Functor &func) { - eigen_assert(dst.rows() == src.rows() && dst.cols() == src.cols()); - typedef evaluator<DstXprType> DstEvaluatorType; typedef evaluator<SrcXprType> SrcEvaluatorType; - DstEvaluatorType dstEvaluator(dst); SrcEvaluatorType srcEvaluator(src); + + Index dstRows = src.rows(); + Index dstCols = src.cols(); + if((dst.rows()!=dstRows) || (dst.cols()!=dstCols)) + dst.resize(dstRows, dstCols); + DstEvaluatorType dstEvaluator(dst); typedef triangular_dense_assignment_kernel< Mode&(Lower|Upper),Mode&(UnitDiag|ZeroDiag|SelfAdjoint),SetOpposite, DstEvaluatorType,SrcEvaluatorType,Functor> Kernel; @@ -800,7 +803,7 @@ void call_triangular_assignment_loop(const DstXprType& dst, const SrcXprType& sr template<int Mode, bool SetOpposite, typename DstXprType, typename SrcXprType> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE -void call_triangular_assignment_loop(const DstXprType& dst, const SrcXprType& src) +void call_triangular_assignment_loop(DstXprType& dst, const SrcXprType& src) { call_triangular_assignment_loop<Mode,SetOpposite>(dst, src, internal::assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar>()); } @@ -936,6 +939,11 @@ struct Assignment<DstXprType, Product<Lhs,Rhs,DefaultProduct>, internal::assign_ typedef Product<Lhs,Rhs,DefaultProduct> SrcXprType; static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,typename SrcXprType::Scalar> &) { + Index dstRows = src.rows(); + Index dstCols = src.cols(); + if((dst.rows()!=dstRows) || (dst.cols()!=dstCols)) + dst.resize(dstRows, dstCols); + dst.setZero(); dst._assignProduct(src, 1); } diff --git a/Eigen/src/Geometry/Homogeneous.h b/Eigen/src/Geometry/Homogeneous.h index 804e5da73..b4d7946e3 100644 --- a/Eigen/src/Geometry/Homogeneous.h +++ b/Eigen/src/Geometry/Homogeneous.h @@ -334,6 +334,11 @@ struct Assignment<DstXprType, Homogeneous<ArgType,Vertical>, internal::assign_op typedef Homogeneous<ArgType,Vertical> SrcXprType; EIGEN_DEVICE_FUNC static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,typename ArgType::Scalar> &) { + Index dstRows = src.rows(); + Index dstCols = src.cols(); + if((dst.rows()!=dstRows) || (dst.cols()!=dstCols)) + dst.resize(dstRows, dstCols); + dst.template topRows<ArgType::RowsAtCompileTime>(src.nestedExpression().rows()) = src.nestedExpression(); dst.row(dst.rows()-1).setOnes(); } @@ -346,6 +351,11 @@ struct Assignment<DstXprType, Homogeneous<ArgType,Horizontal>, internal::assign_ typedef Homogeneous<ArgType,Horizontal> SrcXprType; EIGEN_DEVICE_FUNC static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,typename ArgType::Scalar> &) { + Index dstRows = src.rows(); + Index dstCols = src.cols(); + if((dst.rows()!=dstRows) || (dst.cols()!=dstCols)) + dst.resize(dstRows, dstCols); + dst.template leftCols<ArgType::ColsAtCompileTime>(src.nestedExpression().cols()) = src.nestedExpression(); dst.col(dst.cols()-1).setOnes(); } diff --git a/Eigen/src/IterativeLinearSolvers/SolveWithGuess.h b/Eigen/src/IterativeLinearSolvers/SolveWithGuess.h index 0498db396..0ace45177 100644 --- a/Eigen/src/IterativeLinearSolvers/SolveWithGuess.h +++ b/Eigen/src/IterativeLinearSolvers/SolveWithGuess.h @@ -98,7 +98,11 @@ struct Assignment<DstXprType, SolveWithGuess<DecType,RhsType,GuessType>, interna typedef SolveWithGuess<DecType,RhsType,GuessType> SrcXprType; static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,Scalar> &) { - // FIXME shall we resize dst here? + Index dstRows = src.rows(); + Index dstCols = src.cols(); + if((dst.rows()!=dstRows) || (dst.cols()!=dstCols)) + dst.resize(dstRows, dstCols); + dst = src.guess(); src.dec()._solve_with_guess_impl(src.rhs(), dst/*, src.guess()*/); } diff --git a/Eigen/src/LU/InverseImpl.h b/Eigen/src/LU/InverseImpl.h index 3134632e1..018f99b58 100644 --- a/Eigen/src/LU/InverseImpl.h +++ b/Eigen/src/LU/InverseImpl.h @@ -292,7 +292,11 @@ struct Assignment<DstXprType, Inverse<XprType>, internal::assign_op<typename Dst typedef Inverse<XprType> SrcXprType; static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename XprType::Scalar> &) { - // FIXME shall we resize dst here? + Index dstRows = src.rows(); + Index dstCols = src.cols(); + if((dst.rows()!=dstRows) || (dst.cols()!=dstCols)) + dst.resize(dstRows, dstCols); + const int Size = EIGEN_PLAIN_ENUM_MIN(XprType::ColsAtCompileTime,DstXprType::ColsAtCompileTime); EIGEN_ONLY_USED_FOR_DEBUG(Size); eigen_assert(( (Size<=1) || (Size>4) || (extract_data(src.nestedExpression())!=extract_data(dst))) diff --git a/Eigen/src/SparseCore/SparseAssign.h b/Eigen/src/SparseCore/SparseAssign.h index fa5386599..83776645b 100644 --- a/Eigen/src/SparseCore/SparseAssign.h +++ b/Eigen/src/SparseCore/SparseAssign.h @@ -139,13 +139,16 @@ struct Assignment<DstXprType, SrcXprType, Functor, Sparse2Dense> { static void run(DstXprType &dst, const SrcXprType &src, const Functor &func) { - eigen_assert(dst.rows() == src.rows() && dst.cols() == src.cols()); - if(internal::is_same<Functor,internal::assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar> >::value) dst.setZero(); internal::evaluator<SrcXprType> srcEval(src); + Index dstRows = src.rows(); + Index dstCols = src.cols(); + if((dst.rows()!=dstRows) || (dst.cols()!=dstCols)) + dst.resize(dstRows, dstCols); internal::evaluator<DstXprType> dstEval(dst); + const Index outerEvaluationSize = (internal::evaluator<SrcXprType>::Flags&RowMajorBit) ? src.rows() : src.cols(); for (Index j=0; j<outerEvaluationSize; ++j) for (typename internal::evaluator<SrcXprType>::InnerIterator i(srcEval,j); i; ++i) @@ -161,6 +164,11 @@ struct Assignment<DstXprType, Solve<DecType,RhsType>, internal::assign_op<Scalar typedef Solve<DecType,RhsType> SrcXprType; static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,Scalar> &) { + Index dstRows = src.rows(); + Index dstCols = src.cols(); + if((dst.rows()!=dstRows) || (dst.cols()!=dstCols)) + dst.resize(dstRows, dstCols); + src.dec()._solve_impl(src.rhs(), dst); } }; @@ -179,6 +187,11 @@ struct Assignment<DstXprType, SrcXprType, Functor, Diagonal2Sparse> template<int Options> static void run(SparseMatrix<Scalar,Options,StorageIndex> &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar> &/*func*/) { + Index dstRows = src.rows(); + Index dstCols = src.cols(); + if((dst.rows()!=dstRows) || (dst.cols()!=dstCols)) + dst.resize(dstRows, dstCols); + Index size = src.diagonal().size(); dst.makeCompressed(); dst.resizeNonZeros(size); diff --git a/Eigen/src/SparseCore/SparseProduct.h b/Eigen/src/SparseCore/SparseProduct.h index 7a5ad0635..4cbf68781 100644 --- a/Eigen/src/SparseCore/SparseProduct.h +++ b/Eigen/src/SparseCore/SparseProduct.h @@ -104,6 +104,11 @@ struct Assignment<DstXprType, Product<Lhs,Rhs,AliasFreeProduct>, internal::assig typedef Product<Lhs,Rhs,AliasFreeProduct> SrcXprType; static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar> &) { + Index dstRows = src.rows(); + Index dstCols = src.cols(); + if((dst.rows()!=dstRows) || (dst.cols()!=dstCols)) + dst.resize(dstRows, dstCols); + generic_product_impl<Lhs, Rhs>::evalTo(dst,src.lhs(),src.rhs()); } }; diff --git a/test/product_extra.cpp b/test/product_extra.cpp index e4990ac8c..03d5c3657 100644 --- a/test/product_extra.cpp +++ b/test/product_extra.cpp @@ -256,7 +256,49 @@ Index compute_block_size() return ret; } - +template<typename> +void aliasing_with_resize() +{ + Index m = internal::random<Index>(10,50); + Index n = internal::random<Index>(10,50); + MatrixXd A, B, C(m,n), D(m,m); + VectorXd a, b, c(n); + C.setRandom(); + D.setRandom(); + c.setRandom(); + double s = internal::random<double>(1,10); + + A = C; + B = A * A.transpose(); + A = A * A.transpose(); + VERIFY_IS_APPROX(A,B); + + A = C; + B = (A * A.transpose())/s; + A = (A * A.transpose())/s; + VERIFY_IS_APPROX(A,B); + + A = C; + B = (A * A.transpose()) + D; + A = (A * A.transpose()) + D; + VERIFY_IS_APPROX(A,B); + + A = C; + B = D + (A * A.transpose()); + A = D + (A * A.transpose()); + VERIFY_IS_APPROX(A,B); + + A = C; + B = s * (A * A.transpose()); + A = s * (A * A.transpose()); + VERIFY_IS_APPROX(A,B); + + A = C; + a = c; + b = (A * a)/s; + a = (A * a)/s; + VERIFY_IS_APPROX(a,b); +} template<int> void bug_1308() @@ -318,5 +360,6 @@ void test_product_extra() CALL_SUBTEST_7( compute_block_size<float>() ); CALL_SUBTEST_7( compute_block_size<double>() ); CALL_SUBTEST_7( compute_block_size<std::complex<double> >() ); + CALL_SUBTEST_8( aliasing_with_resize<void>() ); } |