diff options
-rw-r--r-- | Eigen/Array | 1 | ||||
-rw-r--r-- | Eigen/Core | 3 | ||||
-rw-r--r-- | Eigen/Geometry | 1 | ||||
-rw-r--r-- | Eigen/src/Array/PartialRedux.h | 38 | ||||
-rw-r--r-- | Eigen/src/Array/Replicate.h | 159 | ||||
-rw-r--r-- | Eigen/src/Core/Functors.h | 2 | ||||
-rw-r--r-- | Eigen/src/Core/MatrixBase.h | 23 | ||||
-rw-r--r-- | Eigen/src/Core/ReturnByValue.h | 26 | ||||
-rw-r--r-- | Eigen/src/Core/products/GeneralMatrixMatrix.h | 1 | ||||
-rw-r--r-- | Eigen/src/Core/products/GeneralMatrixVector.h | 43 | ||||
-rw-r--r-- | Eigen/src/Core/products/SelfadjointMatrixVector.h | 31 | ||||
-rw-r--r-- | Eigen/src/Core/util/Constants.h | 8 | ||||
-rw-r--r-- | Eigen/src/Core/util/ForwardDeclarations.h | 6 | ||||
-rw-r--r-- | Eigen/src/Core/util/XprHelper.h | 13 | ||||
-rw-r--r-- | Eigen/src/Geometry/Homogeneous.h | 168 | ||||
-rw-r--r-- | test/CMakeLists.txt | 4 | ||||
-rw-r--r-- | test/array_replicate.cpp | 86 | ||||
-rw-r--r-- | test/array_reverse.cpp (renamed from test/reverse.cpp) | 2 | ||||
-rw-r--r-- | test/geo_homogeneous.cpp | 79 | ||||
-rw-r--r-- | test/packetmath.cpp | 4 |
20 files changed, 623 insertions, 75 deletions
diff --git a/Eigen/Array b/Eigen/Array index eb61e7152..57bd39e9e 100644 --- a/Eigen/Array +++ b/Eigen/Array @@ -35,6 +35,7 @@ namespace Eigen { #include "src/Array/PartialRedux.h" #include "src/Array/Random.h" #include "src/Array/Norms.h" +#include "src/Array/Replicate.h" #include "src/Array/Reverse.h" } // namespace Eigen diff --git a/Eigen/Core b/Eigen/Core index e75b90931..89531e8da 100644 --- a/Eigen/Core +++ b/Eigen/Core @@ -129,7 +129,6 @@ namespace Eigen { #include "src/Core/CwiseUnaryOp.h" #include "src/Core/CwiseNullaryOp.h" #include "src/Core/Dot.h" -#include "src/Core/Product.h" #include "src/Core/DiagonalProduct.h" #include "src/Core/SolveTriangular.h" #include "src/Core/MapBase.h" @@ -146,9 +145,11 @@ namespace Eigen { #include "src/Core/Swap.h" #include "src/Core/CommaInitializer.h" #include "src/Core/Part.h" +#include "src/Core/Product.h" #include "src/Core/products/GeneralMatrixMatrix.h" #include "src/Core/products/GeneralMatrixVector.h" #include "src/Core/products/SelfadjointMatrixVector.h" +#include "src/Core/products/SelfadjointRank2Update.h" } // namespace Eigen diff --git a/Eigen/Geometry b/Eigen/Geometry index 617b25eb6..61aa70a91 100644 --- a/Eigen/Geometry +++ b/Eigen/Geometry @@ -32,6 +32,7 @@ namespace Eigen { */ #include "src/Geometry/OrthoMethods.h" +#include "src/Geometry/Homogeneous.h" #include "src/Geometry/RotationBase.h" #include "src/Geometry/Rotation2D.h" #include "src/Geometry/Quaternion.h" diff --git a/Eigen/src/Array/PartialRedux.h b/Eigen/src/Array/PartialRedux.h index d1ed33c38..37fe2b70f 100644 --- a/Eigen/src/Array/PartialRedux.h +++ b/Eigen/src/Array/PartialRedux.h @@ -175,8 +175,6 @@ template<typename ExpressionType, int Direction> class PartialRedux > Type; }; - typedef typename ExpressionType::PlainMatrixType CrossReturnType; - inline PartialRedux(const ExpressionType& matrix) : m_matrix(matrix) {} /** \internal */ @@ -282,9 +280,45 @@ template<typename ExpressionType, int Direction> class PartialRedux { return Reverse<ExpressionType, Direction>( _expression() ); } + +/////////// Geometry module /////////// + + const Homogeneous<ExpressionType,Direction> homogeneous() const; + + const Replicate<ExpressionType,Direction==Vertical?Dynamic:1,Direction==Horizontal?Dynamic:1> + replicate(int factor) const; + + template<int Factor> + const Replicate<ExpressionType,Direction==Vertical?Factor:1,Direction==Horizontal?Factor:1> + replicate() const; + typedef typename ExpressionType::PlainMatrixType CrossReturnType; template<typename OtherDerived> const CrossReturnType cross(const MatrixBase<OtherDerived>& other) const; + + enum { + HNormalized_Size = Direction==Vertical ? ei_traits<ExpressionType>::RowsAtCompileTime + : ei_traits<ExpressionType>::ColsAtCompileTime, + HNormalized_SizeMinusOne = HNormalized_Size==Dynamic ? Dynamic : HNormalized_Size-1 + }; + typedef Block<ExpressionType, + Direction==Vertical ? int(HNormalized_SizeMinusOne) + : int(ei_traits<ExpressionType>::RowsAtCompileTime), + Direction==Horizontal ? int(HNormalized_SizeMinusOne) + : int(ei_traits<ExpressionType>::ColsAtCompileTime)> + HNormalized_Block; + typedef Block<ExpressionType, + Direction==Vertical ? 1 : int(ei_traits<ExpressionType>::RowsAtCompileTime), + Direction==Horizontal ? 1 : int(ei_traits<ExpressionType>::ColsAtCompileTime)> + HNormalized_Factors; + typedef CwiseBinaryOp<ei_scalar_quotient_op<typename ei_traits<ExpressionType>::Scalar>, + NestByValue<HNormalized_Block>, + NestByValue<Replicate<NestByValue<HNormalized_Factors>, + Direction==Vertical ? HNormalized_SizeMinusOne : 1, + Direction==Horizontal ? HNormalized_SizeMinusOne : 1> > > + HNormalizedReturnType; + + const HNormalizedReturnType hnormalized() const; protected: ExpressionTypeNested m_matrix; diff --git a/Eigen/src/Array/Replicate.h b/Eigen/src/Array/Replicate.h new file mode 100644 index 000000000..09d64d067 --- /dev/null +++ b/Eigen/src/Array/Replicate.h @@ -0,0 +1,159 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#ifndef EIGEN_REPLICATE_H +#define EIGEN_REPLICATE_H + +/** \nonstableyet + * \class Replicate + * + * \brief Expression of the multiple replication of a matrix or vector + * + * \param MatrixType the type of the object we are replicating + * + * This class represents an expression of the multiple replication of a matrix or vector. + * It is the return type of MatrixBase::replicate() and most of the time + * this is the only way it is used. + * + * \sa MatrixBase::replicate() + */ +template<typename MatrixType,int RowFactor,int ColFactor> +struct ei_traits<Replicate<MatrixType,RowFactor,ColFactor> > +{ + typedef typename MatrixType::Scalar Scalar; + typedef typename ei_nested<MatrixType>::type MatrixTypeNested; + typedef typename ei_unref<MatrixTypeNested>::type _MatrixTypeNested; + enum { + RowsPlusOne = (MatrixType::RowsAtCompileTime != Dynamic) ? + int(MatrixType::RowsAtCompileTime) + 1 : Dynamic, + ColsPlusOne = (MatrixType::ColsAtCompileTime != Dynamic) ? + int(MatrixType::ColsAtCompileTime) + 1 : Dynamic, + RowsAtCompileTime = RowFactor==Dynamic || MatrixType::RowsAtCompileTime==Dynamic + ? Dynamic + : RowFactor * MatrixType::RowsAtCompileTime, + ColsAtCompileTime = ColFactor==Dynamic || MatrixType::ColsAtCompileTime==Dynamic + ? Dynamic + : ColFactor * MatrixType::ColsAtCompileTime, + MaxRowsAtCompileTime = RowsAtCompileTime, + MaxColsAtCompileTime = ColsAtCompileTime, + Flags = _MatrixTypeNested::Flags & HereditaryBits, + CoeffReadCost = _MatrixTypeNested::CoeffReadCost + }; +}; + +template<typename MatrixType,int RowFactor,int ColFactor> class Replicate + : public MatrixBase<Replicate<MatrixType,RowFactor,ColFactor> > +{ + public: + + EIGEN_GENERIC_PUBLIC_INTERFACE(Replicate) + + inline Replicate(const MatrixType& matrix) + : m_matrix(matrix) + { + ei_assert(RowFactor!=Dynamic && ColFactor!=Dynamic); + } + + inline Replicate(const MatrixType& matrix, int rowFactor, int colFactor) + : m_matrix(matrix), m_rowFactor(rowFactor), m_colFactor(colFactor) + {} + + inline int rows() const { return m_matrix.rows() * m_rowFactor.value(); } + inline int cols() const { return m_matrix.cols() * m_colFactor.value(); } + + inline Scalar coeff(int row, int col) const + { + return m_matrix.coeff(row%m_matrix.rows(), col%m_matrix.cols()); + } + + protected: + const typename MatrixType::Nested m_matrix; + const ei_int_if_dynamic<RowFactor> m_rowFactor; + const ei_int_if_dynamic<ColFactor> m_colFactor; +}; + +/** \nonstableyet + * \return an expression of the replication of \c *this + * + * Example: \include MatrixBase_replicate.cpp + * Output: \verbinclude MatrixBase_replicate.out + * + * \sa PartialRedux::replicate(), MatrixBase::replicate(int,int), class Replicate + */ +template<typename Derived> +template<int RowFactor, int ColFactor> +inline const Replicate<Derived,RowFactor,ColFactor> +MatrixBase<Derived>::replicate() const +{ + return derived(); +} + +/** \nonstableyet + * \return an expression of the replication of \c *this + * + * Example: \include MatrixBase_replicate_int_int.cpp + * Output: \verbinclude MatrixBase_replicate_int_int.out + * + * \sa PartialRedux::replicate(), MatrixBase::replicate<int,int>(), class Replicate + */ +template<typename Derived> +inline const Replicate<Derived,Dynamic,Dynamic> +MatrixBase<Derived>::replicate(int rowFactor,int colFactor) const +{ + return Replicate<Derived,Dynamic,Dynamic>(derived(),rowFactor,colFactor); +} + +/** \nonstableyet + * \return an expression of the replication of each column (or row) of \c *this + * + * Example: \include DirectionWise_replicate_int.cpp + * Output: \verbinclude DirectionWise_replicate_int.out + * + * \sa PartialRedux::replicate(), MatrixBase::replicate(), class Replicate + */ +template<typename ExpressionType, int Direction> +const Replicate<ExpressionType,Direction==Vertical?Dynamic:1,Direction==Horizontal?Dynamic:1> +PartialRedux<ExpressionType,Direction>::replicate(int factor) const +{ + return Replicate<ExpressionType,Direction==Vertical?Dynamic:1,Direction==Horizontal?Dynamic:1> + (_expression(),Direction==Vertical?factor:1,Direction==Horizontal?factor:1); +} + +/** \nonstableyet + * \return an expression of the replication of each column (or row) of \c *this + * + * Example: \include DirectionWise_replicate.cpp + * Output: \verbinclude DirectionWise_replicate.out + * + * \sa PartialRedux::replicate(int), MatrixBase::replicate(), class Replicate + */ +template<typename ExpressionType, int Direction> +template<int Factor> +const Replicate<ExpressionType,Direction==Vertical?Factor:1,Direction==Horizontal?Factor:1> +PartialRedux<ExpressionType,Direction>::replicate() const +{ + return _expression(); +} + +#endif // EIGEN_REPLICATE_H diff --git a/Eigen/src/Core/Functors.h b/Eigen/src/Core/Functors.h index 57a055662..41e94e28d 100644 --- a/Eigen/src/Core/Functors.h +++ b/Eigen/src/Core/Functors.h @@ -299,7 +299,7 @@ struct ei_functor_traits<ei_scalar_imag_op<Scalar> > * indeed it seems better to declare m_other as a PacketScalar and do the ei_pset1() once * in the constructor. However, in practice: * - GCC does not like m_other as a PacketScalar and generate a load every time it needs it - * - one the other hand GCC is able to moves the ei_pset1() away the loop :) + * - on the other hand GCC is able to moves the ei_pset1() away the loop :) * - simpler code ;) * (ICC and gcc 4.4 seems to perform well in both cases, the issue is visible with y = a*x + b*y) */ diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 3a7f1a3e0..9b54ac28c 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -360,10 +360,6 @@ template<typename Derived> class MatrixBase void transposeInPlace(); const AdjointReturnType adjoint() const; - Eigen::Reverse<Derived, BothDirections> reverse(); - const Eigen::Reverse<Derived, BothDirections> reverse() const; - void reverseInPlace(); - RowXpr row(int i); const RowXpr row(int i) const; @@ -589,6 +585,14 @@ template<typename Derived> class MatrixBase select(typename ElseDerived::Scalar thenScalar, const MatrixBase<ElseDerived>& elseMatrix) const; template<int p> RealScalar lpNorm() const; + + template<int RowFactor, int ColFactor> + const Replicate<Derived,RowFactor,ColFactor> replicate() const; + const Replicate<Derived,Dynamic,Dynamic> replicate(int rowFacor,int colFactor) const; + + Eigen::Reverse<Derived, BothDirections> reverse(); + const Eigen::Reverse<Derived, BothDirections> reverse() const; + void reverseInPlace(); /////////// LU module /////////// @@ -620,6 +624,17 @@ template<typename Derived> class MatrixBase PlainMatrixType unitOrthogonal(void) const; Matrix<Scalar,3,1> eulerAngles(int a0, int a1, int a2) const; const ScalarMultipleReturnType operator*(const UniformScaling<Scalar>& s) const; + enum { + SizeMinusOne = SizeAtCompileTime==Dynamic ? Dynamic : SizeAtCompileTime-1 + }; + typedef Block<Derived, + ei_traits<Derived>::ColsAtCompileTime==1 ? SizeMinusOne : 1, + ei_traits<Derived>::ColsAtCompileTime==1 ? 1 : SizeMinusOne> StartMinusOne; + typedef CwiseUnaryOp<ei_scalar_quotient1_op<typename ei_traits<Derived>::Scalar>, + NestByValue<StartMinusOne> > HNormalizedReturnType; + + const HNormalizedReturnType hnormalized() const; + const Homogeneous<Derived,MatrixBase<Derived>::ColsAtCompileTime==1?Vertical:Horizontal> homogeneous() const; /////////// Sparse module /////////// diff --git a/Eigen/src/Core/ReturnByValue.h b/Eigen/src/Core/ReturnByValue.h index e6d2b0d01..3613782ac 100644 --- a/Eigen/src/Core/ReturnByValue.h +++ b/Eigen/src/Core/ReturnByValue.h @@ -28,30 +28,40 @@ /** \class ReturnByValue * */ -template<typename Functor,typename EvalType> -struct ei_traits<ReturnByValue<Functor,EvalType> > : public ei_traits<EvalType> +template<typename Functor, typename _Scalar,int _Rows,int _Cols,int _Options,int _MaxRows,int _MaxCols> +struct ei_traits<ReturnByValue<Functor,Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols> > > + : public ei_traits<Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols> > { enum { - Flags = ei_traits<EvalType>::Flags | EvalBeforeNestingBit + Flags = ei_traits<Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols> >::Flags | EvalBeforeNestingBit }; }; -template<typename Functor,typename EvalType,int n> -struct ei_nested<ReturnByValue<Functor,EvalType>, n, EvalType> +template<typename Functor,typename EvalTypeDerived,int n> +struct ei_nested<ReturnByValue<Functor,MatrixBase<EvalTypeDerived> >, n, EvalTypeDerived> { - typedef EvalType type; + typedef EvalTypeDerived type; }; template<typename Functor, typename EvalType> class ReturnByValue - : public MatrixBase<ReturnByValue<Functor,EvalType> > { public: - EIGEN_GENERIC_PUBLIC_INTERFACE(ReturnByValue) template<typename Dest> inline void evalTo(Dest& dst) const { static_cast<const Functor*>(this)->evalTo(dst); } }; +template<typename Functor, typename _Scalar,int _Rows,int _Cols,int _Options,int _MaxRows,int _MaxCols> + class ReturnByValue<Functor,Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols> > + : public MatrixBase<ReturnByValue<Functor,Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols> > > +{ + public: + EIGEN_GENERIC_PUBLIC_INTERFACE(ReturnByValue) + template<typename Dest> + inline void evalTo(Dest& dst) const + { static_cast<const Functor* const>(this)->evalTo(dst); } +}; + template<typename Derived> template<typename OtherDerived,typename OtherEvalType> Derived& MatrixBase<Derived>::operator=(const ReturnByValue<OtherDerived,OtherEvalType>& other) diff --git a/Eigen/src/Core/products/GeneralMatrixMatrix.h b/Eigen/src/Core/products/GeneralMatrixMatrix.h index 7ec33d7f4..3543de9b8 100644 --- a/Eigen/src/Core/products/GeneralMatrixMatrix.h +++ b/Eigen/src/Core/products/GeneralMatrixMatrix.h @@ -262,6 +262,7 @@ static void ei_cache_friendly_product( blB += 4*nr*PacketSize; blA += 4*mr; } + // process remaining peeled loop for(int k=peeled_kc; k<actual_kc; k++) { PacketType B0, B1, B2, B3, A0, A1; diff --git a/Eigen/src/Core/products/GeneralMatrixVector.h b/Eigen/src/Core/products/GeneralMatrixVector.h index 089db4682..2df6a652d 100644 --- a/Eigen/src/Core/products/GeneralMatrixVector.h +++ b/Eigen/src/Core/products/GeneralMatrixVector.h @@ -38,19 +38,18 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_colmajor_times_vector( const Scalar* lhs, int lhsStride, const RhsType& rhs, Scalar* res) -{ +{asm("#ei_cache_friendly_product_colmajor_times_vector"); #ifdef _EIGEN_ACCUMULATE_PACKETS #error _EIGEN_ACCUMULATE_PACKETS has already been defined #endif - - #define _EIGEN_ACCUMULATE_PACKETS(A0,A13,A2,OFFSET) \ - ei_pstore(&res[j OFFSET], \ - ei_padd(ei_pload(&res[j OFFSET]), \ + #define _EIGEN_ACCUMULATE_PACKETS(A0,A13,A2) \ + ei_pstore(&res[j], \ + ei_padd(ei_pload(&res[j]), \ ei_padd( \ - ei_padd(ei_pmul(ptmp0,EIGEN_CAT(ei_pload , A0)(&lhs0[j OFFSET])), \ - ei_pmul(ptmp1,EIGEN_CAT(ei_pload , A13)(&lhs1[j OFFSET]))), \ - ei_padd(ei_pmul(ptmp2,EIGEN_CAT(ei_pload , A2)(&lhs2[j OFFSET])), \ - ei_pmul(ptmp3,EIGEN_CAT(ei_pload , A13)(&lhs3[j OFFSET]))) ))) + ei_padd(ei_pmul(ptmp0,EIGEN_CAT(ei_ploa , A0)(&lhs0[j])), \ + ei_pmul(ptmp1,EIGEN_CAT(ei_ploa , A13)(&lhs1[j]))), \ + ei_padd(ei_pmul(ptmp2,EIGEN_CAT(ei_ploa , A2)(&lhs2[j])), \ + ei_pmul(ptmp3,EIGEN_CAT(ei_ploa , A13)(&lhs3[j]))) ))) typedef typename ei_packet_traits<Scalar>::type Packet; const int PacketSize = sizeof(Packet)/sizeof(Scalar); @@ -125,11 +124,11 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_colmajor_times_vector( { case AllAligned: for (int j = alignedStart; j<alignedSize; j+=PacketSize) - _EIGEN_ACCUMULATE_PACKETS(EIGEN_EMPTY,EIGEN_EMPTY,EIGEN_EMPTY,EIGEN_EMPTY); + _EIGEN_ACCUMULATE_PACKETS(d,d,d); break; case EvenAligned: for (int j = alignedStart; j<alignedSize; j+=PacketSize) - _EIGEN_ACCUMULATE_PACKETS(EIGEN_EMPTY,u,EIGEN_EMPTY,EIGEN_EMPTY); + _EIGEN_ACCUMULATE_PACKETS(d,du,d); break; case FirstAligned: if(peels>1) @@ -165,11 +164,11 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_colmajor_times_vector( } } for (int j = peeledSize; j<alignedSize; j+=PacketSize) - _EIGEN_ACCUMULATE_PACKETS(EIGEN_EMPTY,u,u,EIGEN_EMPTY); + _EIGEN_ACCUMULATE_PACKETS(d,du,du); break; default: for (int j = alignedStart; j<alignedSize; j+=PacketSize) - _EIGEN_ACCUMULATE_PACKETS(u,u,u,EIGEN_EMPTY); + _EIGEN_ACCUMULATE_PACKETS(du,du,du); break; } } @@ -233,12 +232,12 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_rowmajor_times_vector( #error _EIGEN_ACCUMULATE_PACKETS has already been defined #endif - #define _EIGEN_ACCUMULATE_PACKETS(A0,A13,A2,OFFSET) {\ + #define _EIGEN_ACCUMULATE_PACKETS(A0,A13,A2) {\ Packet b = ei_pload(&rhs[j]); \ - ptmp0 = ei_pmadd(b, EIGEN_CAT(ei_pload,A0) (&lhs0[j]), ptmp0); \ - ptmp1 = ei_pmadd(b, EIGEN_CAT(ei_pload,A13)(&lhs1[j]), ptmp1); \ - ptmp2 = ei_pmadd(b, EIGEN_CAT(ei_pload,A2) (&lhs2[j]), ptmp2); \ - ptmp3 = ei_pmadd(b, EIGEN_CAT(ei_pload,A13)(&lhs3[j]), ptmp3); } + ptmp0 = ei_pmadd(b, EIGEN_CAT(ei_ploa,A0) (&lhs0[j]), ptmp0); \ + ptmp1 = ei_pmadd(b, EIGEN_CAT(ei_ploa,A13)(&lhs1[j]), ptmp1); \ + ptmp2 = ei_pmadd(b, EIGEN_CAT(ei_ploa,A2) (&lhs2[j]), ptmp2); \ + ptmp3 = ei_pmadd(b, EIGEN_CAT(ei_ploa,A13)(&lhs3[j]), ptmp3); } typedef typename ei_packet_traits<Scalar>::type Packet; const int PacketSize = sizeof(Packet)/sizeof(Scalar); @@ -319,11 +318,11 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_rowmajor_times_vector( { case AllAligned: for (int j = alignedStart; j<alignedSize; j+=PacketSize) - _EIGEN_ACCUMULATE_PACKETS(EIGEN_EMPTY,EIGEN_EMPTY,EIGEN_EMPTY,EIGEN_EMPTY); + _EIGEN_ACCUMULATE_PACKETS(d,d,d); break; case EvenAligned: for (int j = alignedStart; j<alignedSize; j+=PacketSize) - _EIGEN_ACCUMULATE_PACKETS(EIGEN_EMPTY,u,EIGEN_EMPTY,EIGEN_EMPTY); + _EIGEN_ACCUMULATE_PACKETS(d,du,d); break; case FirstAligned: if (peels>1) @@ -362,11 +361,11 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_rowmajor_times_vector( } } for (int j = peeledSize; j<alignedSize; j+=PacketSize) - _EIGEN_ACCUMULATE_PACKETS(EIGEN_EMPTY,u,u,EIGEN_EMPTY); + _EIGEN_ACCUMULATE_PACKETS(d,du,du); break; default: for (int j = alignedStart; j<alignedSize; j+=PacketSize) - _EIGEN_ACCUMULATE_PACKETS(u,u,u,EIGEN_EMPTY); + _EIGEN_ACCUMULATE_PACKETS(du,du,du); break; } tmp0 += ei_predux(ptmp0); diff --git a/Eigen/src/Core/products/SelfadjointMatrixVector.h b/Eigen/src/Core/products/SelfadjointMatrixVector.h index d539ef414..3c3ab7222 100644 --- a/Eigen/src/Core/products/SelfadjointMatrixVector.h +++ b/Eigen/src/Core/products/SelfadjointMatrixVector.h @@ -140,37 +140,6 @@ static EIGEN_DONT_INLINE void ei_product_selfadjoint_vector( } res[j] += t2; } - /* - // colmajor - upper - for (int j=0;j<size;j++) - { - register const Scalar* __restrict__ A0 = lhs + j*lhsStride; - - Scalar t1 = rhs[j]; - Scalar t2 = 0; - for (int i=0; i<j; i+=PacketSize) { - res[i] += t1 * A0[i]; - t2 += A0[i] * rhs[i]; - } - res[j] += t1 * A0[j]; - res[j] += t2; - } - - // rowmajor - lower - for (int j=0;j<size;j++) - { - register const Scalar* __restrict__ A0 = lhs + j*lhsStride; - - Scalar t1 = rhs[j]; - Scalar t2 = 0; - for (int i=0; i<j; i+=PacketSize) { - res[i] += t1 * A0[i]; - t2 += A0[i] * rhs[i]; - } - res[j] += t1 * A0[j]; - res[j] += t2; - } - */ } diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h index 603e595e8..5d0aafda9 100644 --- a/Eigen/src/Core/util/Constants.h +++ b/Eigen/src/Core/util/Constants.h @@ -239,4 +239,12 @@ enum { HasDirectAccess = DirectAccessBit }; +enum TransformTraits { + Affine = 0x1, + Isometry = 0x2, + AffineSquare = Affine|0x10, + AffineCompact = Affine|0x20, + Projective = 0x30 +}; + #endif // EIGEN_CONSTANTS_H diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index 3dc6762f8..a8953ead3 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -40,7 +40,6 @@ template<typename MatrixType, int BlockRows=Dynamic, int BlockCols=Dynamic, int int _DirectAccessStatus = ei_traits<MatrixType>::Flags&DirectAccessBit ? DirectAccessBit : ei_traits<MatrixType>::Flags&SparseBit> class Block; template<typename MatrixType> class Transpose; -template<typename MatrixType, int Direction = BothDirections> class Reverse; template<typename MatrixType> class Conjugate; template<typename NullaryOp, typename MatrixType> class CwiseNullaryOp; template<typename UnaryOp, typename MatrixType> class CwiseUnaryOp; @@ -104,6 +103,8 @@ void ei_cache_friendly_product( template<typename ConditionMatrixType, typename ThenMatrixType, typename ElseMatrixType> class Select; template<typename MatrixType, typename BinaryOp, int Direction> class PartialReduxExpr; template<typename ExpressionType, int Direction> class PartialRedux; +template<typename MatrixType,int RowFactor,int ColFactor> class Replicate; +template<typename MatrixType, int Direction = BothDirections> class Reverse; template<typename MatrixType> class LU; template<typename MatrixType> class QR; @@ -117,11 +118,12 @@ template<typename Lhs, typename Rhs> class Cross; template<typename Scalar> class Quaternion; template<typename Scalar> class Rotation2D; template<typename Scalar> class AngleAxis; -template<typename Scalar,int Dim> class Transform; +template<typename Scalar,int Dim/*,int Mode=AffineSquare*/> class Transform; template <typename _Scalar, int _AmbientDim> class ParametrizedLine; template <typename _Scalar, int _AmbientDim> class Hyperplane; template<typename Scalar,int Dim> class Translation; template<typename Scalar> class UniformScaling; +template<typename MatrixType,int Direction> class Homogeneous; // Sparse module: template<typename Lhs, typename Rhs, int ProductMode> class SparseProduct; diff --git a/Eigen/src/Core/util/XprHelper.h b/Eigen/src/Core/util/XprHelper.h index 12d6f9a3a..82e32f1f5 100644 --- a/Eigen/src/Core/util/XprHelper.h +++ b/Eigen/src/Core/util/XprHelper.h @@ -211,6 +211,19 @@ template<typename ExpressionType, int RowsOrSize=Dynamic, int Cols=Dynamic> stru typedef Block<ExpressionType, RowsOrSize, Cols> Type; }; +template<typename ExpressionType> struct HNormalizedReturnType { + + enum { + SizeAtCompileTime = ExpressionType::SizeAtCompileTime, + SizeMinusOne = SizeAtCompileTime==Dynamic ? Dynamic : SizeAtCompileTime-1 + }; + typedef Block<ExpressionType, + ei_traits<ExpressionType>::ColsAtCompileTime==1 ? SizeMinusOne : 1, + ei_traits<ExpressionType>::ColsAtCompileTime==1 ? 1 : SizeMinusOne> StartMinusOne; + typedef CwiseUnaryOp<ei_scalar_quotient1_op<typename ei_traits<ExpressionType>::Scalar>, + NestByValue<StartMinusOne> > Type; +}; + template<typename CurrentType, typename NewType> struct ei_cast_return_type { typedef typename ei_meta_if<ei_is_same_type<CurrentType,NewType>::ret,const CurrentType&,NewType>::ret type; diff --git a/Eigen/src/Geometry/Homogeneous.h b/Eigen/src/Geometry/Homogeneous.h new file mode 100644 index 000000000..7e1114474 --- /dev/null +++ b/Eigen/src/Geometry/Homogeneous.h @@ -0,0 +1,168 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#ifndef EIGEN_HOMOGENEOUS_H +#define EIGEN_HOMOGENEOUS_H + +/** \geometry_module \ingroup Geometry_Module + * \nonstableyet + * \class Homogeneous + * + * \brief Expression of one (or a set of) homogeneous vector(s) + * + * \param MatrixType the type of the object in which we are making homogeneous + * + * This class represents an expression of one (or a set of) homogeneous vector(s). + * It is the return type of MatrixBase::homogeneous() and most of the time + * this is the only way it is used. + * + * \sa MatrixBase::homogeneous() + */ +template<typename MatrixType,int Direction> +struct ei_traits<Homogeneous<MatrixType,Direction> > +{ + typedef typename MatrixType::Scalar Scalar; + typedef typename ei_nested<MatrixType>::type MatrixTypeNested; + typedef typename ei_unref<MatrixTypeNested>::type _MatrixTypeNested; + enum { + RowsPlusOne = (MatrixType::RowsAtCompileTime != Dynamic) ? + int(MatrixType::RowsAtCompileTime) + 1 : Dynamic, + ColsPlusOne = (MatrixType::ColsAtCompileTime != Dynamic) ? + int(MatrixType::ColsAtCompileTime) + 1 : Dynamic, + RowsAtCompileTime = Direction==Vertical ? RowsPlusOne : MatrixType::RowsAtCompileTime, + ColsAtCompileTime = Direction==Horizontal ? ColsPlusOne : MatrixType::ColsAtCompileTime, + MaxRowsAtCompileTime = RowsAtCompileTime, + MaxColsAtCompileTime = ColsAtCompileTime, + Flags = _MatrixTypeNested::Flags & HereditaryBits, + CoeffReadCost = _MatrixTypeNested::CoeffReadCost + }; +}; + +template<typename MatrixType,int Direction> class Homogeneous + : public MatrixBase<Homogeneous<MatrixType,Direction> > +{ + public: + + EIGEN_GENERIC_PUBLIC_INTERFACE(Homogeneous) + + inline Homogeneous(const MatrixType& matrix) + : m_matrix(matrix) + {} + + inline int rows() const { return m_matrix.rows() + (Direction==Vertical ? 1 : 0); } + inline int cols() const { return m_matrix.cols() + (Direction==Horizontal ? 1 : 0); } + + inline Scalar coeff(int row, int col) const + { + if( (Direction==Vertical && row==m_matrix.rows()) + || (Direction==Horizontal && col==m_matrix.cols())) + return 1; + return m_matrix.coeff(row, col); + } + + protected: + const typename MatrixType::Nested m_matrix; +}; + +/** \geometry_module + * \nonstableyet + * \return an expression of the equivalent homogeneous vector + * + * \vectoronly + * + * Example: \include MatrixBase_homogeneous.cpp + * Output: \verbinclude MatrixBase_homogeneous.out + * + * \sa class Homogeneous + */ +template<typename Derived> +inline const Homogeneous<Derived,MatrixBase<Derived>::ColsAtCompileTime==1?Vertical:Horizontal> +MatrixBase<Derived>::homogeneous() const +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived); + return derived(); +} + +/** \geometry_module + * \nonstableyet + * \returns a matrix expression of homogeneous column (or row) vectors + * + * Example: \include PartialRedux_homogeneous.cpp + * Output: \verbinclude PartialRedux_homogeneous.out + * + * \sa MatrixBase::homogeneous() */ +template<typename ExpressionType, int Direction> +inline const Homogeneous<ExpressionType,Direction> +PartialRedux<ExpressionType,Direction>::homogeneous() const +{ + return _expression(); +} + +/** \geometry_module + * \nonstableyet + * \returns an expression of the homogeneous normalized vector of \c *this + * + * Example: \include MatrixBase_hnormalized.cpp + * Output: \verbinclude MatrixBase_hnormalized.out + * + * \sa PartialRedux::hnormalized() */ +template<typename Derived> +inline const typename MatrixBase<Derived>::HNormalizedReturnType +MatrixBase<Derived>::hnormalized() const +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived); + return StartMinusOne(derived(),0,0, + ColsAtCompileTime==1?size()-1:1, + ColsAtCompileTime==1?1:size()-1).nestByValue() / coeff(size()-1); +} + +/** \geometry_module + * \nonstableyet + * \returns an expression of the homogeneous normalized vector of \c *this + * + * Example: \include DirectionWise_hnormalized.cpp + * Output: \verbinclude DirectionWise_hnormalized.out + * + * \sa MatrixBase::hnormalized() */ +template<typename ExpressionType, int Direction> +inline const typename PartialRedux<ExpressionType,Direction>::HNormalizedReturnType +PartialRedux<ExpressionType,Direction>::hnormalized() const +{ + return HNormalized_Block(_expression(),0,0, + Direction==Vertical ? _expression().rows()-1 : _expression().rows(), + Direction==Horizontal ? _expression().cols()-1 : _expression().cols()).nestByValue() + .cwise()/ + Replicate<NestByValue<HNormalized_Factors>, + Direction==Vertical ? HNormalized_SizeMinusOne : 1, + Direction==Horizontal ? HNormalized_SizeMinusOne : 1> + (HNormalized_Factors(_expression(), + Direction==Vertical ? _expression().rows()-1:0, + Direction==Horizontal ? _expression().cols()-1:0, + Direction==Vertical ? 1 : _expression().rows(), + Direction==Horizontal ? 1 : _expression().cols()).nestByValue(), + Direction==Vertical ? _expression().rows()-1 : 1, + Direction==Horizontal ? _expression().cols()-1 : 1).nestByValue(); +} + +#endif // EIGEN_HOMOGENEOUS_H diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 09c0d9a93..cc959d69a 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -105,6 +105,8 @@ ei_add_test(commainitializer) ei_add_test(smallvectors) ei_add_test(map) ei_add_test(array) +ei_add_test(array_replicate) +ei_add_test(array_reverse) ei_add_test(triangular) ei_add_test(cholesky " " "${GSL_LIBRARIES}") ei_add_test(lu ${EI_OFLAG}) @@ -114,6 +116,7 @@ ei_add_test(qr) ei_add_test(eigensolver " " "${GSL_LIBRARIES}") ei_add_test(svd) ei_add_test(geo_orthomethods) +ei_add_test(geo_homogeneous) ei_add_test(geo_quaternion) ei_add_test(geo_transformations) ei_add_test(geo_eulerangles) @@ -135,7 +138,6 @@ ei_add_test(sparse_vector) ei_add_test(sparse_basic) ei_add_test(sparse_product) ei_add_test(sparse_solvers " " "${SPARSE_LIBS}") -ei_add_test(reverse) ei_add_property(EIGEN_TESTING_SUMMARY "CXX: ${CMAKE_CXX_COMPILER}\n") diff --git a/test/array_replicate.cpp b/test/array_replicate.cpp new file mode 100644 index 000000000..e8d0f33cb --- /dev/null +++ b/test/array_replicate.cpp @@ -0,0 +1,86 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#include "main.h" +#include <Eigen/Array> + +template<typename MatrixType> void replicate(const MatrixType& m) +{ + /* this test covers the following files: + Replicate.cpp + */ + + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; + typedef Matrix<Scalar, Dynamic, Dynamic> MatrixX; + typedef Matrix<Scalar, Dynamic, 1> VectorX; + + int rows = m.rows(); + int cols = m.cols(); + + MatrixType m1 = MatrixType::Random(rows, cols), + m2 = MatrixType::Random(rows, cols); + + VectorType v1 = VectorType::Random(rows); + + MatrixX x1, x2; + VectorX vx1; + + int f1 = ei_random<int>(1,10), + f2 = ei_random<int>(1,10); + + x1.resize(rows*f1,cols*f2); + for(int j=0; j<f2; j++) + for(int i=0; i<f1; i++) + x1.block(i*rows,j*cols,rows,cols) = m1; + VERIFY_IS_APPROX(x1, m1.replicate(f1,f2)); + + x2.resize(2*rows,3*cols); + x2 << m2, m2, m2, + m2, m2, m2; + VERIFY_IS_APPROX(x2, (m2.template replicate<2,3>())); + + x2.resize(rows,f1); + for (int j=0; j<f1; ++j) + x2.col(j) = v1; + VERIFY_IS_APPROX(x2, v1.rowwise().replicate(f1)); + + vx1.resize(rows*f2); + for (int j=0; j<f2; ++j) + vx1.segment(j*rows,rows) = v1; + VERIFY_IS_APPROX(vx1, v1.colwise().replicate(f2)); +} + +void test_array_replicate() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST( replicate(Matrix<float, 1, 1>()) ); + CALL_SUBTEST( replicate(Vector2f()) ); + CALL_SUBTEST( replicate(Vector3d()) ); + CALL_SUBTEST( replicate(Vector4f()) ); + CALL_SUBTEST( replicate(VectorXf(16)) ); + CALL_SUBTEST( replicate(VectorXcd(10)) ); + } +} diff --git a/test/reverse.cpp b/test/array_reverse.cpp index ce1a48ed8..c69ff73e8 100644 --- a/test/reverse.cpp +++ b/test/array_reverse.cpp @@ -160,7 +160,7 @@ template<typename MatrixType> void reverse(const MatrixType& m) */ } -void test_reverse() +void test_array_reverse() { for(int i = 0; i < g_repeat; i++) { CALL_SUBTEST( reverse(Matrix<float, 1, 1>()) ); diff --git a/test/geo_homogeneous.cpp b/test/geo_homogeneous.cpp new file mode 100644 index 000000000..aea6f746e --- /dev/null +++ b/test/geo_homogeneous.cpp @@ -0,0 +1,79 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#include "main.h" +#include <Eigen/Geometry> + +template<typename Scalar,int Size> void homogeneous(void) +{ + /* this test covers the following files: + Homogeneous.h + */ + + typedef Matrix<Scalar,Size,Size> MatrixType; + typedef Matrix<Scalar,Size,1> VectorType; + + typedef Matrix<Scalar,Size+1,Size> HMatrixType; + typedef Matrix<Scalar,Size+1,1> HVectorType; + + Scalar largeEps = test_precision<Scalar>(); + if (ei_is_same_type<Scalar,float>::ret) + largeEps = 1e-3f; + + Scalar eps = ei_random<Scalar>() * 1e-2; + + VectorType v0 = VectorType::Random(), + v1 = VectorType::Random(), + ones = VectorType::Ones(); + + HVectorType hv0 = HVectorType::Random(), + hv1 = HVectorType::Random(); + + MatrixType m0 = MatrixType::Random(), + m1 = MatrixType::Random(); + + HMatrixType hm0 = HMatrixType::Random(), + hm1 = HMatrixType::Random(); + + hv0 << v0, 1; + VERIFY_IS_APPROX(v0.homogeneous(), hv0); + VERIFY_IS_APPROX(v0, hv0.hnormalized()); + + hm0 << m0, ones.transpose(); + VERIFY_IS_APPROX(m0.colwise().homogeneous(), hm0); + VERIFY_IS_APPROX(m0, hm0.colwise().hnormalized()); + hm0.row(Size-1).setRandom(); + for(int j=0; j<Size; ++j) + m0.col(j) = hm0.col(j).start(Size) / hm0(Size,j); + VERIFY_IS_APPROX(m0, hm0.colwise().hnormalized()); +} + +void test_geo_homogeneous() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST(( homogeneous<float,1>() )); + CALL_SUBTEST(( homogeneous<double,3>() )); + CALL_SUBTEST(( homogeneous<double,8>() )); + } +} diff --git a/test/packetmath.cpp b/test/packetmath.cpp index f56532f91..c7694b660 100644 --- a/test/packetmath.cpp +++ b/test/packetmath.cpp @@ -137,8 +137,8 @@ template<typename Scalar> void packetmath() ref[0] = data1[0]; for (int i=0; i<PacketSize; ++i) - ref[0] = std::min(ref[0],data1[i]); - VERIFY(ei_isApprox(ref[0], ei_predux_min(ei_pload(data1))) && "ei_predux_max"); + ref[0] = std::max(ref[0],data1[i]); + VERIFY(ei_isApprox(ref[0], ei_predux_max(ei_pload(data1))) && "ei_predux_max"); for (int j=0; j<PacketSize; ++j) { |