diff options
Diffstat (limited to 'Eigen/src')
53 files changed, 787 insertions, 414 deletions
diff --git a/Eigen/src/Array/Array.h b/Eigen/src/Array/Array.h index 533d638a4..91a091152 100644 --- a/Eigen/src/Array/Array.h +++ b/Eigen/src/Array/Array.h @@ -41,7 +41,7 @@ class Array EIGEN_DENSE_PUBLIC_INTERFACE(Array) enum { Options = _Options }; - typedef typename Base::PlainMatrixType PlainMatrixType; + typedef typename Base::PlainObject PlainObject; protected: using Base::m_storage; @@ -61,7 +61,7 @@ class Array * the usage of 'using'. This should be done only for operator=. */ template<typename OtherDerived> - EIGEN_STRONG_INLINE Array& operator=(const AnyMatrixBase<OtherDerived> &other) + EIGEN_STRONG_INLINE Array& operator=(const EigenBase<OtherDerived> &other) { return Base::operator=(other); } @@ -196,9 +196,9 @@ class Array other.evalTo(*this); } - /** \sa MatrixBase::operator=(const AnyMatrixBase<OtherDerived>&) */ + /** \sa MatrixBase::operator=(const EigenBase<OtherDerived>&) */ template<typename OtherDerived> - EIGEN_STRONG_INLINE Array(const AnyMatrixBase<OtherDerived> &other) + EIGEN_STRONG_INLINE Array(const EigenBase<OtherDerived> &other) : Base(other.derived().rows() * other.derived().cols(), other.derived().rows(), other.derived().cols()) { Base::_check_template_params(); diff --git a/Eigen/src/Array/ArrayBase.h b/Eigen/src/Array/ArrayBase.h index 21f6fefb1..97807e5fc 100644 --- a/Eigen/src/Array/ArrayBase.h +++ b/Eigen/src/Array/ArrayBase.h @@ -97,7 +97,7 @@ template<typename Derived> class ArrayBase /** \internal the plain matrix type corresponding to this expression. Note that is not necessarily * exactly the return type of eval(): in the case of plain matrices, the return type of eval() is a const * reference to a matrix, not a matrix! It is however guaranteed that the return type of eval() is either - * PlainMatrixType or const PlainMatrixType&. + * PlainObject or const PlainObject&. */ typedef Array<typename ei_traits<Derived>::Scalar, ei_traits<Derived>::RowsAtCompileTime, @@ -105,7 +105,7 @@ template<typename Derived> class ArrayBase AutoAlign | (ei_traits<Derived>::Flags&RowMajorBit ? RowMajor : ColMajor), ei_traits<Derived>::MaxRowsAtCompileTime, ei_traits<Derived>::MaxColsAtCompileTime - > PlainMatrixType; + > PlainObject; /** \internal Represents a matrix with all coefficients equal to one another*/ diff --git a/Eigen/src/Array/VectorwiseOp.h b/Eigen/src/Array/VectorwiseOp.h index 697a07d32..06d999b14 100644 --- a/Eigen/src/Array/VectorwiseOp.h +++ b/Eigen/src/Array/VectorwiseOp.h @@ -462,7 +462,7 @@ template<typename ExpressionType, int Direction> class VectorwiseOp const Homogeneous<ExpressionType,Direction> homogeneous() const; - typedef typename ExpressionType::PlainMatrixType CrossReturnType; + typedef typename ExpressionType::PlainObject CrossReturnType; template<typename OtherDerived> const CrossReturnType cross(const MatrixBase<OtherDerived>& other) const; diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index b794b0c43..8699fe7e0 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -62,14 +62,21 @@ template<typename _MatrixType> class LDLT typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType; typedef Matrix<int, 1, MatrixType::RowsAtCompileTime> IntRowVectorType; - /** - * \brief Default Constructor. - * - * The default constructor is useful in cases in which the user intends to - * perform decompositions via LDLT::compute(const MatrixType&). - */ + /** \brief Default Constructor. + * + * The default constructor is useful in cases in which the user intends to + * perform decompositions via LDLT::compute(const MatrixType&). + */ LDLT() : m_matrix(), m_p(), m_transpositions(), m_isInitialized(false) {} + /** \brief Default Constructor with memory preallocation + * + * Like the default constructor but with preallocation of the internal data + * according to the specified problem \a size. + * \sa LDLT() + */ + LDLT(int size) : m_matrix(size,size), m_p(size), m_transpositions(size), m_isInitialized(false) {} + LDLT(const MatrixType& matrix) : m_matrix(matrix.rows(), matrix.cols()), m_p(matrix.rows()), @@ -148,6 +155,8 @@ template<typename _MatrixType> class LDLT return m_matrix; } + MatrixType reconstructedMatrix() const; + inline int rows() const { return m_matrix.rows(); } inline int cols() const { return m_matrix.cols(); } @@ -175,6 +184,10 @@ LDLT<MatrixType>& LDLT<MatrixType>::compute(const MatrixType& a) m_matrix = a; + m_p.resize(size); + m_transpositions.resize(size); + m_isInitialized = false; + if (size <= 1) { m_p.setZero(); m_transpositions.setZero(); @@ -202,11 +215,8 @@ LDLT<MatrixType>& LDLT<MatrixType>::compute(const MatrixType& a) { // The biggest overall is the point of reference to which further diagonals // are compared; if any diagonal is negligible compared - // to the largest overall, the algorithm bails. This cutoff is suggested - // in "Analysis of the Cholesky Decomposition of a Semi-definite Matrix" by - // Nicholas J. Higham. Also see "Accuracy and Stability of Numerical - // Algorithms" page 217, also by Higham. - cutoff = ei_abs(NumTraits<Scalar>::epsilon() * RealScalar(size) * biggest_in_corner); + // to the largest overall, the algorithm bails. + cutoff = ei_abs(NumTraits<Scalar>::epsilon() * biggest_in_corner); m_sign = ei_real(m_matrix.diagonal().coeff(index_of_biggest_in_corner)) > 0 ? 1 : -1; } @@ -231,26 +241,21 @@ LDLT<MatrixType>& LDLT<MatrixType>::compute(const MatrixType& a) continue; } - RealScalar Djj = ei_real(m_matrix.coeff(j,j) - m_matrix.row(j).head(j) - .dot(m_matrix.col(j).head(j))); + RealScalar Djj = ei_real(m_matrix.coeff(j,j) - m_matrix.row(j).head(j).dot(m_matrix.col(j).head(j))); m_matrix.coeffRef(j,j) = Djj; - // Finish early if the matrix is not full rank. - if(ei_abs(Djj) < cutoff) - { - for(int i = j; i < size; i++) m_transpositions.coeffRef(i) = i; - break; - } - int endSize = size - j - 1; if (endSize > 0) { _temporary.tail(endSize).noalias() = m_matrix.block(j+1,0, endSize, j) * m_matrix.col(j).head(j).conjugate(); m_matrix.row(j).tail(endSize) = m_matrix.row(j).tail(endSize).conjugate() - - _temporary.tail(endSize).transpose(); + - _temporary.tail(endSize).transpose(); - m_matrix.col(j).tail(endSize) = m_matrix.row(j).tail(endSize) / Djj; + if(ei_abs(Djj) > cutoff) + { + m_matrix.col(j).tail(endSize) = m_matrix.row(j).tail(endSize) / Djj; + } } } @@ -315,14 +320,39 @@ bool LDLT<MatrixType>::solveInPlace(MatrixBase<Derived> &bAndX) const return true; } +/** \returns the matrix represented by the decomposition, + * i.e., it returns the product: P^T L D L^* P. + * This function is provided for debug purpose. */ +template<typename MatrixType> +MatrixType LDLT<MatrixType>::reconstructedMatrix() const +{ + ei_assert(m_isInitialized && "LDLT is not initialized."); + const int size = m_matrix.rows(); + MatrixType res(size,size); + res.setIdentity(); + + // PI + for(int i = 0; i < size; ++i) res.row(m_transpositions.coeff(i)).swap(res.row(i)); + // L^* P + res = matrixL().adjoint() * res; + // D(L^*P) + res = vectorD().asDiagonal() * res; + // L(DL^*P) + res = matrixL() * res; + // P^T (LDL^*P) + for (int i = size-1; i >= 0; --i) res.row(m_transpositions.coeff(i)).swap(res.row(i)); + + return res; +} + /** \cholesky_module * \returns the Cholesky decomposition with full pivoting without square root of \c *this */ template<typename Derived> -inline const LDLT<typename MatrixBase<Derived>::PlainMatrixType> +inline const LDLT<typename MatrixBase<Derived>::PlainObject> MatrixBase<Derived>::ldlt() const { - return derived(); + return LDLT<PlainObject>(derived()); } #endif // EIGEN_LDLT_H diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h index 8a149a316..2e8df7661 100644 --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h @@ -117,7 +117,7 @@ template<typename _MatrixType, int _UpLo> class LLT && "LLT::solve(): invalid number of rows of the right hand side matrix b"); return ei_solve_retval<LLT, Rhs>(*this, b.derived()); } - + template<typename Derived> bool solveInPlace(MatrixBase<Derived> &bAndX) const; @@ -133,6 +133,8 @@ template<typename _MatrixType, int _UpLo> class LLT return m_matrix; } + MatrixType reconstructedMatrix() const; + inline int rows() const { return m_matrix.rows(); } inline int cols() const { return m_matrix.cols(); } @@ -295,24 +297,34 @@ bool LLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const return true; } +/** \returns the matrix represented by the decomposition, + * i.e., it returns the product: L L^*. + * This function is provided for debug purpose. */ +template<typename MatrixType, int _UpLo> +MatrixType LLT<MatrixType,_UpLo>::reconstructedMatrix() const +{ + ei_assert(m_isInitialized && "LLT is not initialized."); + return matrixL() * matrixL().adjoint().toDenseMatrix(); +} + /** \cholesky_module * \returns the LLT decomposition of \c *this */ template<typename Derived> -inline const LLT<typename MatrixBase<Derived>::PlainMatrixType> +inline const LLT<typename MatrixBase<Derived>::PlainObject> MatrixBase<Derived>::llt() const { - return LLT<PlainMatrixType>(derived()); + return LLT<PlainObject>(derived()); } /** \cholesky_module * \returns the LLT decomposition of \c *this */ template<typename MatrixType, unsigned int UpLo> -inline const LLT<typename SelfAdjointView<MatrixType, UpLo>::PlainMatrixType, UpLo> +inline const LLT<typename SelfAdjointView<MatrixType, UpLo>::PlainObject, UpLo> SelfAdjointView<MatrixType, UpLo>::llt() const { - return LLT<PlainMatrixType,UpLo>(m_matrix); + return LLT<PlainObject,UpLo>(m_matrix); } #endif // EIGEN_LLT_H diff --git a/Eigen/src/Core/BandMatrix.h b/Eigen/src/Core/BandMatrix.h index 538e6dd76..432df0b34 100644 --- a/Eigen/src/Core/BandMatrix.h +++ b/Eigen/src/Core/BandMatrix.h @@ -57,7 +57,7 @@ struct ei_traits<BandMatrix<_Scalar,Rows,Cols,Supers,Subs,Options> > }; template<typename _Scalar, int Rows, int Cols, int Supers, int Subs, int Options> -class BandMatrix : public AnyMatrixBase<BandMatrix<_Scalar,Rows,Cols,Supers,Subs,Options> > +class BandMatrix : public EigenBase<BandMatrix<_Scalar,Rows,Cols,Supers,Subs,Options> > { public: diff --git a/Eigen/src/Core/DenseBase.h b/Eigen/src/Core/DenseBase.h index 2078f023b..5682d7278 100644 --- a/Eigen/src/Core/DenseBase.h +++ b/Eigen/src/Core/DenseBase.h @@ -40,7 +40,7 @@ template<typename Derived> class DenseBase : public ei_special_scalar_op_base<Derived,typename ei_traits<Derived>::Scalar, typename NumTraits<typename ei_traits<Derived>::Scalar>::Real> #else - : public AnyMatrixBase<Derived> + : public EigenBase<Derived> #endif // not EIGEN_PARSED_BY_DOXYGEN { public: @@ -53,8 +53,8 @@ template<typename Derived> class DenseBase typedef typename ei_traits<Derived>::Scalar Scalar; typedef typename ei_packet_traits<Scalar>::type PacketScalar; - using AnyMatrixBase<Derived>::derived; - using AnyMatrixBase<Derived>::const_cast_derived; + using EigenBase<Derived>::derived; + using EigenBase<Derived>::const_cast_derived; #endif // not EIGEN_PARSED_BY_DOXYGEN enum { @@ -292,13 +292,13 @@ template<typename Derived> class DenseBase Derived& operator=(const DenseBase& other); template<typename OtherDerived> - Derived& operator=(const AnyMatrixBase<OtherDerived> &other); + Derived& operator=(const EigenBase<OtherDerived> &other); template<typename OtherDerived> - Derived& operator+=(const AnyMatrixBase<OtherDerived> &other); + Derived& operator+=(const EigenBase<OtherDerived> &other); template<typename OtherDerived> - Derived& operator-=(const AnyMatrixBase<OtherDerived> &other); + Derived& operator-=(const EigenBase<OtherDerived> &other); template<typename OtherDerived> Derived& operator=(const ReturnByValue<OtherDerived>& func); diff --git a/Eigen/src/Core/DenseStorageBase.h b/Eigen/src/Core/DenseStorageBase.h index e93e439e6..89e6e7112 100644 --- a/Eigen/src/Core/DenseStorageBase.h +++ b/Eigen/src/Core/DenseStorageBase.h @@ -44,7 +44,7 @@ class DenseStorageBase : public _Base<Derived> public: enum { Options = _Options }; typedef _Base<Derived> Base; - typedef typename Base::PlainMatrixType PlainMatrixType; + typedef typename Base::PlainObject PlainObject; typedef typename Base::Scalar Scalar; typedef typename Base::PacketScalar PacketScalar; using Base::RowsAtCompileTime; @@ -338,19 +338,19 @@ class DenseStorageBase : public _Base<Derived> // EIGEN_INITIALIZE_BY_ZERO_IF_THAT_OPTION_IS_ENABLED } - /** \copydoc MatrixBase::operator=(const AnyMatrixBase<OtherDerived>&) + /** \copydoc MatrixBase::operator=(const EigenBase<OtherDerived>&) */ template<typename OtherDerived> - EIGEN_STRONG_INLINE Derived& operator=(const AnyMatrixBase<OtherDerived> &other) + EIGEN_STRONG_INLINE Derived& operator=(const EigenBase<OtherDerived> &other) { resize(other.derived().rows(), other.derived().cols()); Base::operator=(other.derived()); return this->derived(); } - /** \sa MatrixBase::operator=(const AnyMatrixBase<OtherDerived>&) */ + /** \sa MatrixBase::operator=(const EigenBase<OtherDerived>&) */ template<typename OtherDerived> - EIGEN_STRONG_INLINE DenseStorageBase(const AnyMatrixBase<OtherDerived> &other) + EIGEN_STRONG_INLINE DenseStorageBase(const EigenBase<OtherDerived> &other) : m_storage(other.derived().rows() * other.derived().cols(), other.derived().rows(), other.derived().cols()) { _check_template_params(); @@ -527,7 +527,7 @@ struct ei_conservative_resize_like_impl { if (_this.rows() == rows && _this.cols() == cols) return; EIGEN_STATIC_ASSERT_DYNAMIC_SIZE(Derived) - typename Derived::PlainMatrixType tmp(rows,cols); + typename Derived::PlainObject tmp(rows,cols); const int common_rows = std::min(rows, _this.rows()); const int common_cols = std::min(cols, _this.cols()); tmp.block(0,0,common_rows,common_cols) = _this.block(0,0,common_rows,common_cols); @@ -546,7 +546,7 @@ struct ei_conservative_resize_like_impl EIGEN_STATIC_ASSERT_DYNAMIC_SIZE(Derived) EIGEN_STATIC_ASSERT_DYNAMIC_SIZE(OtherDerived) - typename Derived::PlainMatrixType tmp(other); + typename Derived::PlainObject tmp(other); const int common_rows = std::min(tmp.rows(), _this.rows()); const int common_cols = std::min(tmp.cols(), _this.cols()); tmp.block(0,0,common_rows,common_cols) = _this.block(0,0,common_rows,common_cols); @@ -560,7 +560,7 @@ struct ei_conservative_resize_like_impl<Derived,OtherDerived,true> static void run(DenseBase<Derived>& _this, int size) { if (_this.size() == size) return; - typename Derived::PlainMatrixType tmp(size); + typename Derived::PlainObject tmp(size); const int common_size = std::min<int>(_this.size(),size); tmp.segment(0,common_size) = _this.segment(0,common_size); _this.derived().swap(tmp); @@ -571,7 +571,7 @@ struct ei_conservative_resize_like_impl<Derived,OtherDerived,true> if (_this.rows() == other.rows() && _this.cols() == other.cols()) return; // segment(...) will check whether Derived/OtherDerived are vectors! - typename Derived::PlainMatrixType tmp(other); + typename Derived::PlainObject tmp(other); const int common_size = std::min<int>(_this.size(),tmp.size()); tmp.segment(0,common_size) = _this.segment(0,common_size); _this.derived().swap(tmp); diff --git a/Eigen/src/Core/DiagonalMatrix.h b/Eigen/src/Core/DiagonalMatrix.h index 08c046611..774b0d7ae 100644 --- a/Eigen/src/Core/DiagonalMatrix.h +++ b/Eigen/src/Core/DiagonalMatrix.h @@ -28,7 +28,7 @@ #ifndef EIGEN_PARSED_BY_DOXYGEN template<typename Derived> -class DiagonalBase : public AnyMatrixBase<Derived> +class DiagonalBase : public EigenBase<Derived> { public: typedef typename ei_traits<Derived>::DiagonalVectorType DiagonalVectorType; diff --git a/Eigen/src/Core/Dot.h b/Eigen/src/Core/Dot.h index f0c520b1f..201bd23ca 100644 --- a/Eigen/src/Core/Dot.h +++ b/Eigen/src/Core/Dot.h @@ -299,7 +299,7 @@ inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real MatrixBase< * \sa norm(), normalize() */ template<typename Derived> -inline const typename MatrixBase<Derived>::PlainMatrixType +inline const typename MatrixBase<Derived>::PlainObject MatrixBase<Derived>::normalized() const { typedef typename ei_nested<Derived>::type Nested; diff --git a/Eigen/src/Core/AnyMatrixBase.h b/Eigen/src/Core/EigenBase.h index a5d1cfe9f..cf1ce4376 100644 --- a/Eigen/src/Core/AnyMatrixBase.h +++ b/Eigen/src/Core/EigenBase.h @@ -23,21 +23,21 @@ // License and a copy of the GNU General Public License along with // Eigen. If not, see <http://www.gnu.org/licenses/>. -#ifndef EIGEN_ANYMATRIXBASE_H -#define EIGEN_ANYMATRIXBASE_H +#ifndef EIGEN_EIGENBASE_H +#define EIGEN_EIGENBASE_H /** Common base class for all classes T such that MatrixBase has an operator=(T) and a constructor MatrixBase(T). * - * In other words, an AnyMatrixBase object is an object that can be copied into a MatrixBase. + * In other words, an EigenBase object is an object that can be copied into a MatrixBase. * * Besides MatrixBase-derived classes, this also includes special matrix classes such as diagonal matrices, etc. * * Notice that this class is trivial, it is only used to disambiguate overloaded functions. */ -template<typename Derived> struct AnyMatrixBase +template<typename Derived> struct EigenBase { -// typedef typename ei_plain_matrix_type<Derived>::type PlainMatrixType; +// typedef typename ei_plain_matrix_type<Derived>::type PlainObject; /** \returns a reference to the derived object */ Derived& derived() { return *static_cast<Derived*>(this); } @@ -45,7 +45,7 @@ template<typename Derived> struct AnyMatrixBase const Derived& derived() const { return *static_cast<const Derived*>(this); } inline Derived& const_cast_derived() const - { return *static_cast<Derived*>(const_cast<AnyMatrixBase*>(this)); } + { return *static_cast<Derived*>(const_cast<EigenBase*>(this)); } /** \returns the number of rows. \sa cols(), RowsAtCompileTime */ inline int rows() const { return derived().rows(); } @@ -61,7 +61,7 @@ template<typename Derived> struct AnyMatrixBase { // This is the default implementation, // derived class can reimplement it in a more optimized way. - typename Dest::PlainMatrixType res(rows(),cols()); + typename Dest::PlainObject res(rows(),cols()); evalTo(res); dst += res; } @@ -71,7 +71,7 @@ template<typename Derived> struct AnyMatrixBase { // This is the default implementation, // derived class can reimplement it in a more optimized way. - typename Dest::PlainMatrixType res(rows(),cols()); + typename Dest::PlainObject res(rows(),cols()); evalTo(res); dst -= res; } @@ -108,7 +108,7 @@ template<typename Derived> struct AnyMatrixBase */ template<typename Derived> template<typename OtherDerived> -Derived& DenseBase<Derived>::operator=(const AnyMatrixBase<OtherDerived> &other) +Derived& DenseBase<Derived>::operator=(const EigenBase<OtherDerived> &other) { other.derived().evalTo(derived()); return derived(); @@ -116,7 +116,7 @@ Derived& DenseBase<Derived>::operator=(const AnyMatrixBase<OtherDerived> &other) template<typename Derived> template<typename OtherDerived> -Derived& DenseBase<Derived>::operator+=(const AnyMatrixBase<OtherDerived> &other) +Derived& DenseBase<Derived>::operator+=(const EigenBase<OtherDerived> &other) { other.derived().addTo(derived()); return derived(); @@ -124,7 +124,7 @@ Derived& DenseBase<Derived>::operator+=(const AnyMatrixBase<OtherDerived> &other template<typename Derived> template<typename OtherDerived> -Derived& DenseBase<Derived>::operator-=(const AnyMatrixBase<OtherDerived> &other) +Derived& DenseBase<Derived>::operator-=(const EigenBase<OtherDerived> &other) { other.derived().subTo(derived()); return derived(); @@ -137,7 +137,7 @@ Derived& DenseBase<Derived>::operator-=(const AnyMatrixBase<OtherDerived> &other template<typename Derived> template<typename OtherDerived> inline Derived& -MatrixBase<Derived>::operator*=(const AnyMatrixBase<OtherDerived> &other) +MatrixBase<Derived>::operator*=(const EigenBase<OtherDerived> &other) { other.derived().applyThisOnTheRight(derived()); return derived(); @@ -146,7 +146,7 @@ MatrixBase<Derived>::operator*=(const AnyMatrixBase<OtherDerived> &other) /** replaces \c *this by \c *this * \a other. It is equivalent to MatrixBase::operator*=() */ template<typename Derived> template<typename OtherDerived> -inline void MatrixBase<Derived>::applyOnTheRight(const AnyMatrixBase<OtherDerived> &other) +inline void MatrixBase<Derived>::applyOnTheRight(const EigenBase<OtherDerived> &other) { other.derived().applyThisOnTheRight(derived()); } @@ -154,9 +154,9 @@ inline void MatrixBase<Derived>::applyOnTheRight(const AnyMatrixBase<OtherDerive /** replaces \c *this by \c *this * \a other. */ template<typename Derived> template<typename OtherDerived> -inline void MatrixBase<Derived>::applyOnTheLeft(const AnyMatrixBase<OtherDerived> &other) +inline void MatrixBase<Derived>::applyOnTheLeft(const EigenBase<OtherDerived> &other) { other.derived().applyThisOnTheLeft(derived()); } -#endif // EIGEN_ANYMATRIXBASE_H +#endif // EIGEN_EIGENBASE_H diff --git a/Eigen/src/Core/Flagged.h b/Eigen/src/Core/Flagged.h index 9d14aceaa..9413b74fa 100644 --- a/Eigen/src/Core/Flagged.h +++ b/Eigen/src/Core/Flagged.h @@ -110,7 +110,7 @@ template<typename ExpressionType, unsigned int Added, unsigned int Removed> clas const ExpressionType& _expression() const { return m_matrix; } template<typename OtherDerived> - typename ExpressionType::PlainMatrixType solveTriangular(const MatrixBase<OtherDerived>& other) const; + typename ExpressionType::PlainObject solveTriangular(const MatrixBase<OtherDerived>& other) const; template<typename OtherDerived> void solveTriangularInPlace(const MatrixBase<OtherDerived>& other) const; diff --git a/Eigen/src/Core/Matrix.h b/Eigen/src/Core/Matrix.h index b494b2f00..eae2711f4 100644 --- a/Eigen/src/Core/Matrix.h +++ b/Eigen/src/Core/Matrix.h @@ -139,7 +139,7 @@ class Matrix EIGEN_DENSE_PUBLIC_INTERFACE(Matrix) - typedef typename Base::PlainMatrixType PlainMatrixType; + typedef typename Base::PlainObject PlainObject; enum { NeedsToAlign = (!(Options&DontAlign)) && SizeAtCompileTime!=Dynamic && ((sizeof(Scalar)*SizeAtCompileTime)%16)==0 }; @@ -181,10 +181,10 @@ class Matrix /** * \brief Copies the generic expression \a other into *this. - * \copydetails DenseBase::operator=(const AnyMatrixBase<OtherDerived> &other) + * \copydetails DenseBase::operator=(const EigenBase<OtherDerived> &other) */ template<typename OtherDerived> - EIGEN_STRONG_INLINE Matrix& operator=(const AnyMatrixBase<OtherDerived> &other) + EIGEN_STRONG_INLINE Matrix& operator=(const EigenBase<OtherDerived> &other) { return Base::operator=(other); } @@ -297,10 +297,10 @@ class Matrix } /** \brief Copy constructor for generic expressions. - * \sa MatrixBase::operator=(const AnyMatrixBase<OtherDerived>&) + * \sa MatrixBase::operator=(const EigenBase<OtherDerived>&) */ template<typename OtherDerived> - EIGEN_STRONG_INLINE Matrix(const AnyMatrixBase<OtherDerived> &other) + EIGEN_STRONG_INLINE Matrix(const EigenBase<OtherDerived> &other) : Base(other.derived().rows() * other.derived().cols(), other.derived().rows(), other.derived().cols()) { Base::_check_template_params(); diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 229195046..9c62163ba 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -121,7 +121,7 @@ template<typename Derived> class MatrixBase * * This is not necessarily exactly the return type of eval(). In the case of plain matrices, * the return type of eval() is a const reference to a matrix, not a matrix! It is however guaranteed - * that the return type of eval() is either PlainMatrixType or const PlainMatrixType&. + * that the return type of eval() is either PlainObject or const PlainObject&. */ typedef Matrix<typename ei_traits<Derived>::Scalar, ei_traits<Derived>::RowsAtCompileTime, @@ -129,8 +129,7 @@ template<typename Derived> class MatrixBase AutoAlign | (ei_traits<Derived>::Flags&RowMajorBit ? RowMajor : ColMajor), ei_traits<Derived>::MaxRowsAtCompileTime, ei_traits<Derived>::MaxColsAtCompileTime - > PlainMatrixType; - // typedef typename ei_plain_matrix_type<Derived>::type PlainMatrixType; + > PlainObject; #ifndef EIGEN_PARSED_BY_DOXYGEN /** \internal Represents a matrix with all coefficients equal to one another*/ @@ -193,13 +192,13 @@ template<typename Derived> class MatrixBase lazyProduct(const MatrixBase<OtherDerived> &other) const; template<typename OtherDerived> - Derived& operator*=(const AnyMatrixBase<OtherDerived>& other); + Derived& operator*=(const EigenBase<OtherDerived>& other); template<typename OtherDerived> - void applyOnTheLeft(const AnyMatrixBase<OtherDerived>& other); + void applyOnTheLeft(const EigenBase<OtherDerived>& other); template<typename OtherDerived> - void applyOnTheRight(const AnyMatrixBase<OtherDerived>& other); + void applyOnTheRight(const EigenBase<OtherDerived>& other); template<typename DiagonalDerived> const DiagonalProduct<Derived, DiagonalDerived, OnTheRight> @@ -212,7 +211,7 @@ template<typename Derived> class MatrixBase RealScalar stableNorm() const; RealScalar blueNorm() const; RealScalar hypotNorm() const; - const PlainMatrixType normalized() const; + const PlainObject normalized() const; void normalize(); const AdjointReturnType adjoint() const; @@ -301,9 +300,9 @@ template<typename Derived> class MatrixBase /////////// LU module /////////// - const FullPivLU<PlainMatrixType> fullPivLu() const; - const PartialPivLU<PlainMatrixType> partialPivLu() const; - const PartialPivLU<PlainMatrixType> lu() const; + const FullPivLU<PlainObject> fullPivLu() const; + const PartialPivLU<PlainObject> partialPivLu() const; + const PartialPivLU<PlainObject> lu() const; const ei_inverse_impl<Derived> inverse() const; template<typename ResultType> void computeInverseAndDetWithCheck( @@ -322,29 +321,29 @@ template<typename Derived> class MatrixBase /////////// Cholesky module /////////// - const LLT<PlainMatrixType> llt() const; - const LDLT<PlainMatrixType> ldlt() const; + const LLT<PlainObject> llt() const; + const LDLT<PlainObject> ldlt() const; /////////// QR module /////////// - const HouseholderQR<PlainMatrixType> householderQr() const; - const ColPivHouseholderQR<PlainMatrixType> colPivHouseholderQr() const; - const FullPivHouseholderQR<PlainMatrixType> fullPivHouseholderQr() const; + const HouseholderQR<PlainObject> householderQr() const; + const ColPivHouseholderQR<PlainObject> colPivHouseholderQr() const; + const FullPivHouseholderQR<PlainObject> fullPivHouseholderQr() const; EigenvaluesReturnType eigenvalues() const; RealScalar operatorNorm() const; /////////// SVD module /////////// - SVD<PlainMatrixType> svd() const; + SVD<PlainObject> svd() const; /////////// Geometry module /////////// template<typename OtherDerived> - PlainMatrixType cross(const MatrixBase<OtherDerived>& other) const; + PlainObject cross(const MatrixBase<OtherDerived>& other) const; template<typename OtherDerived> - PlainMatrixType cross3(const MatrixBase<OtherDerived>& other) const; - PlainMatrixType unitOrthogonal(void) const; + PlainObject cross3(const MatrixBase<OtherDerived>& other) const; + PlainObject unitOrthogonal(void) const; Matrix<Scalar,3,1> eulerAngles(int a0, int a1, int a2) const; const ScalarMultipleReturnType operator*(const UniformScaling<Scalar>& s) const; enum { diff --git a/Eigen/src/Core/PermutationMatrix.h b/Eigen/src/Core/PermutationMatrix.h index 284baf678..46884dc3f 100644 --- a/Eigen/src/Core/PermutationMatrix.h +++ b/Eigen/src/Core/PermutationMatrix.h @@ -47,7 +47,7 @@ * \sa class DiagonalMatrix */ template<int SizeAtCompileTime, int MaxSizeAtCompileTime = SizeAtCompileTime> class PermutationMatrix; -template<typename PermutationType, typename MatrixType, int Side> struct ei_permut_matrix_product_retval; +template<typename PermutationType, typename MatrixType, int Side, bool Transposed=false> struct ei_permut_matrix_product_retval; template<int SizeAtCompileTime, int MaxSizeAtCompileTime> struct ei_traits<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime> > @@ -55,7 +55,7 @@ struct ei_traits<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime> > {}; template<int SizeAtCompileTime, int MaxSizeAtCompileTime> -class PermutationMatrix : public AnyMatrixBase<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime> > +class PermutationMatrix : public EigenBase<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime> > { public: @@ -132,6 +132,9 @@ class PermutationMatrix : public AnyMatrixBase<PermutationMatrix<SizeAtCompileTi /** \returns the number of columns */ inline int cols() const { return m_indices.size(); } + /** \returns the size of a side of the respective square matrix, i.e., the number of indices */ + inline int size() const { return m_indices.size(); } + #ifndef EIGEN_PARSED_BY_DOXYGEN template<typename DenseDerived> void evalTo(MatrixBase<DenseDerived>& other) const @@ -144,7 +147,7 @@ class PermutationMatrix : public AnyMatrixBase<PermutationMatrix<SizeAtCompileTi /** \returns a Matrix object initialized from this permutation matrix. Notice that it * is inefficient to return this Matrix object by value. For efficiency, favor using - * the Matrix constructor taking AnyMatrixBase objects. + * the Matrix constructor taking EigenBase objects. */ DenseMatrixType toDenseMatrix() const { @@ -213,16 +216,29 @@ class PermutationMatrix : public AnyMatrixBase<PermutationMatrix<SizeAtCompileTi return *this; } - /**** inversion and multiplication helpers to hopefully get RVO ****/ + /** \returns the inverse permutation matrix. + * + * \note \note_try_to_help_rvo + */ + inline Transpose<PermutationMatrix> inverse() const + { return *this; } + /** \returns the tranpose permutation matrix. + * + * \note \note_try_to_help_rvo + */ + inline Transpose<PermutationMatrix> transpose() const + { return *this; } + + /**** multiplication helpers to hopefully get RVO ****/ #ifndef EIGEN_PARSED_BY_DOXYGEN - protected: - enum Inverse_t {Inverse}; - PermutationMatrix(Inverse_t, const PermutationMatrix& other) - : m_indices(other.m_indices.size()) + template<int OtherSize, int OtherMaxSize> + PermutationMatrix(const Transpose<PermutationMatrix<OtherSize,OtherMaxSize> >& other) + : m_indices(other.nestedPermutation().size()) { - for (int i=0; i<rows();++i) m_indices.coeffRef(other.m_indices.coeff(i)) = i; + for (int i=0; i<rows();++i) m_indices.coeffRef(other.nestedPermutation().indices().coeff(i)) = i; } + protected: enum Product_t {Product}; PermutationMatrix(Product_t, const PermutationMatrix& lhs, const PermutationMatrix& rhs) : m_indices(lhs.m_indices.size()) @@ -233,12 +249,7 @@ class PermutationMatrix : public AnyMatrixBase<PermutationMatrix<SizeAtCompileTi #endif public: - /** \returns the inverse permutation matrix. - * - * \note \note_try_to_help_rvo - */ - inline PermutationMatrix inverse() const - { return PermutationMatrix(Inverse, *this); } + /** \returns the product permutation matrix. * * \note \note_try_to_help_rvo @@ -247,6 +258,22 @@ class PermutationMatrix : public AnyMatrixBase<PermutationMatrix<SizeAtCompileTi inline PermutationMatrix operator*(const PermutationMatrix<OtherSize, OtherMaxSize>& other) const { return PermutationMatrix(Product, *this, other); } + /** \returns the product of a permutation with another inverse permutation. + * + * \note \note_try_to_help_rvo + */ + template<int OtherSize, int OtherMaxSize> + inline PermutationMatrix operator*(const Transpose<PermutationMatrix<OtherSize,OtherMaxSize> >& other) const + { return PermutationMatrix(Product, *this, other.eval()); } + + /** \returns the product of an inverse permutation with another permutation. + * + * \note \note_try_to_help_rvo + */ + template<int OtherSize, int OtherMaxSize> friend + inline PermutationMatrix operator*(const Transpose<PermutationMatrix<OtherSize,OtherMaxSize> >& other, const PermutationMatrix& perm) + { return PermutationMatrix(Product, other.eval(), perm); } + protected: IndicesType m_indices; @@ -277,15 +304,15 @@ operator*(const PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime> &perm (permutation, matrix.derived()); } -template<typename PermutationType, typename MatrixType, int Side> -struct ei_traits<ei_permut_matrix_product_retval<PermutationType, MatrixType, Side> > +template<typename PermutationType, typename MatrixType, int Side, bool Transposed> +struct ei_traits<ei_permut_matrix_product_retval<PermutationType, MatrixType, Side, Transposed> > { - typedef typename MatrixType::PlainMatrixType ReturnMatrixType; + typedef typename MatrixType::PlainObject ReturnType; }; -template<typename PermutationType, typename MatrixType, int Side> +template<typename PermutationType, typename MatrixType, int Side, bool Transposed> struct ei_permut_matrix_product_retval - : public ReturnByValue<ei_permut_matrix_product_retval<PermutationType, MatrixType, Side> > + : public ReturnByValue<ei_permut_matrix_product_retval<PermutationType, MatrixType, Side, Transposed> > { typedef typename ei_cleantype<typename MatrixType::Nested>::type MatrixTypeNestedCleaned; @@ -299,21 +326,46 @@ struct ei_permut_matrix_product_retval template<typename Dest> inline void evalTo(Dest& dst) const { const int n = Side==OnTheLeft ? rows() : cols(); - for(int i = 0; i < n; ++i) + + if(ei_is_same_type<MatrixTypeNestedCleaned,Dest>::ret && ei_extract_data(dst) == ei_extract_data(m_matrix)) { - Block< - Dest, - Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, - Side==OnTheRight ? 1 : Dest::ColsAtCompileTime - >(dst, Side==OnTheLeft ? m_permutation.indices().coeff(i) : i) - - = - - Block< - MatrixTypeNestedCleaned, - Side==OnTheLeft ? 1 : MatrixType::RowsAtCompileTime, - Side==OnTheRight ? 1 : MatrixType::ColsAtCompileTime - >(m_matrix, Side==OnTheRight ? m_permutation.indices().coeff(i) : i); + // apply the permutation inplace + Matrix<bool,PermutationType::RowsAtCompileTime,1,0,PermutationType::MaxRowsAtCompileTime> mask(m_permutation.size()); + mask.fill(false); + int r = 0; + while(r < m_permutation.size()) + { + // search for the next seed + while(r<m_permutation.size() && mask[r]) r++; + if(r>=m_permutation.size()) + break; + // we got one, let's follow it until we are back to the seed + int k0 = r++; + int kPrev = k0; + mask.coeffRef(k0) = true; + for(int k=m_permutation.indices().coeff(k0); k!=k0; k=m_permutation.indices().coeff(k)) + { + Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>(dst, k) + .swap(Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime> + (dst,((Side==OnTheLeft) ^ Transposed) ? k0 : kPrev)); + + mask.coeffRef(k) = true; + kPrev = k; + } + } + } + else + { + for(int i = 0; i < n; ++i) + { + Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime> + (dst, ((Side==OnTheLeft) ^ Transposed) ? m_permutation.indices().coeff(i) : i) + + = + + Block<MatrixTypeNestedCleaned,Side==OnTheLeft ? 1 : MatrixType::RowsAtCompileTime,Side==OnTheRight ? 1 : MatrixType::ColsAtCompileTime> + (m_matrix, ((Side==OnTheRight) ^ Transposed) ? m_permutation.indices().coeff(i) : i); + } } } @@ -322,4 +374,78 @@ struct ei_permut_matrix_product_retval const typename MatrixType::Nested m_matrix; }; +/* Template partial specialization for transposed/inverse permutations */ + +template<int SizeAtCompileTime, int MaxSizeAtCompileTime> +struct ei_traits<Transpose<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime> > > + : ei_traits<Matrix<int,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> > +{}; + +template<int SizeAtCompileTime, int MaxSizeAtCompileTime> +class Transpose<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime> > + : public EigenBase<Transpose<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime> > > +{ + typedef PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime> PermutationType; + typedef typename PermutationType::IndicesType IndicesType; + public: + + #ifndef EIGEN_PARSED_BY_DOXYGEN + typedef ei_traits<PermutationType> Traits; + typedef Matrix<int,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> + DenseMatrixType; + enum { + Flags = Traits::Flags, + CoeffReadCost = Traits::CoeffReadCost, + RowsAtCompileTime = Traits::RowsAtCompileTime, + ColsAtCompileTime = Traits::ColsAtCompileTime, + MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime, + MaxColsAtCompileTime = Traits::MaxColsAtCompileTime + }; + typedef typename Traits::Scalar Scalar; + #endif + + Transpose(const PermutationType& p) : m_permutation(p) {} + + inline int rows() const { return m_permutation.rows(); } + inline int cols() const { return m_permutation.cols(); } + + #ifndef EIGEN_PARSED_BY_DOXYGEN + template<typename DenseDerived> + void evalTo(MatrixBase<DenseDerived>& other) const + { + other.setZero(); + for (int i=0; i<rows();++i) + other.coeffRef(i, m_permutation.indices().coeff(i)) = typename DenseDerived::Scalar(1); + } + #endif + + /** \return the equivalent permutation matrix */ + PermutationType eval() const { return *this; } + + DenseMatrixType toDenseMatrix() const { return *this; } + + /** \returns the matrix with the inverse permutation applied to the columns. + */ + template<typename Derived> friend + inline const ei_permut_matrix_product_retval<PermutationType, Derived, OnTheRight, true> + operator*(const MatrixBase<Derived>& matrix, const Transpose& trPerm) + { + return ei_permut_matrix_product_retval<PermutationType, Derived, OnTheRight, true>(trPerm.m_permutation, matrix.derived()); + } + + /** \returns the matrix with the inverse permutation applied to the rows. + */ + template<typename Derived> + inline const ei_permut_matrix_product_retval<PermutationType, Derived, OnTheLeft, true> + operator*(const MatrixBase<Derived>& matrix) const + { + return ei_permut_matrix_product_retval<PermutationType, Derived, OnTheLeft, true>(m_permutation, matrix.derived()); + } + + const PermutationType& nestedPermutation() const { return m_permutation; } + + protected: + const PermutationType& m_permutation; +}; + #endif // EIGEN_PERMUTATIONMATRIX_H diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h index af05773ee..236e4f130 100644 --- a/Eigen/src/Core/Product.h +++ b/Eigen/src/Core/Product.h @@ -50,8 +50,8 @@ class GeneralProduct; template<int Rows, int Cols, int Depth> struct ei_product_type_selector; enum { - Large = Dynamic, - Small = Dynamic/2 + Large = 2, + Small = 3 }; template<typename Lhs, typename Rhs> struct ei_product_type @@ -95,10 +95,10 @@ template<> struct ei_product_type_selector<Small, Large, 1> template<> struct ei_product_type_selector<Large, Small, 1> { enum { ret = LazyCoeffBasedProductMode }; }; template<> struct ei_product_type_selector<1, Large,Small> { enum { ret = GemvProduct }; }; template<> struct ei_product_type_selector<1, Large,Large> { enum { ret = GemvProduct }; }; -template<> struct ei_product_type_selector<1, Small,Large> { enum { ret = GemvProduct }; }; +template<> struct ei_product_type_selector<1, Small,Large> { enum { ret = CoeffBasedProductMode }; }; template<> struct ei_product_type_selector<Large,1, Small> { enum { ret = GemvProduct }; }; template<> struct ei_product_type_selector<Large,1, Large> { enum { ret = GemvProduct }; }; -template<> struct ei_product_type_selector<Small,1, Large> { enum { ret = GemvProduct }; }; +template<> struct ei_product_type_selector<Small,1, Large> { enum { ret = CoeffBasedProductMode }; }; template<> struct ei_product_type_selector<Small,Small,Large> { enum { ret = GemmProduct }; }; template<> struct ei_product_type_selector<Large,Small,Large> { enum { ret = GemmProduct }; }; template<> struct ei_product_type_selector<Small,Large,Large> { enum { ret = GemmProduct }; }; diff --git a/Eigen/src/Core/ProductBase.h b/Eigen/src/Core/ProductBase.h index 481e7c760..789aecfb6 100644 --- a/Eigen/src/Core/ProductBase.h +++ b/Eigen/src/Core/ProductBase.h @@ -88,7 +88,7 @@ class ProductBase : public MatrixBase<Derived> public: - typedef typename Base::PlainMatrixType PlainMatrixType; + typedef typename Base::PlainObject PlainObject; ProductBase(const Lhs& lhs, const Rhs& rhs) : m_lhs(lhs), m_rhs(rhs) @@ -116,8 +116,8 @@ class ProductBase : public MatrixBase<Derived> const _LhsNested& lhs() const { return m_lhs; } const _RhsNested& rhs() const { return m_rhs; } - // Implicit convertion to the nested type (trigger the evaluation of the product) - operator const PlainMatrixType& () const + // Implicit conversion to the nested type (trigger the evaluation of the product) + operator const PlainObject& () const { m_result.resize(m_lhs.rows(), m_rhs.cols()); this->evalTo(m_result); @@ -139,7 +139,7 @@ class ProductBase : public MatrixBase<Derived> const LhsNested m_lhs; const RhsNested m_rhs; - mutable PlainMatrixType m_result; + mutable PlainObject m_result; private: @@ -152,10 +152,10 @@ class ProductBase : public MatrixBase<Derived> // here we need to overload the nested rule for products // such that the nested type is a const reference to a plain matrix -template<typename Lhs, typename Rhs, int Mode, int N, typename PlainMatrixType> -struct ei_nested<GeneralProduct<Lhs,Rhs,Mode>, N, PlainMatrixType> +template<typename Lhs, typename Rhs, int Mode, int N, typename PlainObject> +struct ei_nested<GeneralProduct<Lhs,Rhs,Mode>, N, PlainObject> { - typedef PlainMatrixType const& type; + typedef PlainObject const& type; }; template<typename NestedProduct> diff --git a/Eigen/src/Core/ReturnByValue.h b/Eigen/src/Core/ReturnByValue.h index 8d45fc31b..160b973bd 100644 --- a/Eigen/src/Core/ReturnByValue.h +++ b/Eigen/src/Core/ReturnByValue.h @@ -31,13 +31,13 @@ */ template<typename Derived> struct ei_traits<ReturnByValue<Derived> > - : public ei_traits<typename ei_traits<Derived>::ReturnMatrixType> + : public ei_traits<typename ei_traits<Derived>::ReturnType> { enum { // We're disabling the DirectAccess because e.g. the constructor of // the Block-with-DirectAccess expression requires to have a coeffRef method. // Also, we don't want to have to implement the stride stuff. - Flags = (ei_traits<typename ei_traits<Derived>::ReturnMatrixType>::Flags + Flags = (ei_traits<typename ei_traits<Derived>::ReturnType>::Flags | EvalBeforeNestingBit) & ~DirectAccessBit }; }; @@ -46,18 +46,18 @@ struct ei_traits<ReturnByValue<Derived> > * So the only way that nesting it in an expression can work, is by evaluating it into a plain matrix. * So ei_nested always gives the plain return matrix type. */ -template<typename Derived,int n,typename PlainMatrixType> -struct ei_nested<ReturnByValue<Derived>, n, PlainMatrixType> +template<typename Derived,int n,typename PlainObject> +struct ei_nested<ReturnByValue<Derived>, n, PlainObject> { - typedef typename ei_traits<Derived>::ReturnMatrixType type; + typedef typename ei_traits<Derived>::ReturnType type; }; template<typename Derived> class ReturnByValue - : public ei_traits<Derived>::ReturnMatrixType::template MakeBase<ReturnByValue<Derived> >::Type + : public ei_traits<Derived>::ReturnType::template MakeBase<ReturnByValue<Derived> >::Type { public: - typedef typename ei_traits<Derived>::ReturnMatrixType ReturnMatrixType; - typedef typename ReturnMatrixType::template MakeBase<ReturnByValue<Derived> >::Type Base; + typedef typename ei_traits<Derived>::ReturnType ReturnType; + typedef typename ReturnType::template MakeBase<ReturnByValue<Derived> >::Type Base; EIGEN_DENSE_PUBLIC_INTERFACE(ReturnByValue) template<typename Dest> diff --git a/Eigen/src/Core/SelfAdjointView.h b/Eigen/src/Core/SelfAdjointView.h index add5a3afb..0b57b9968 100644 --- a/Eigen/src/Core/SelfAdjointView.h +++ b/Eigen/src/Core/SelfAdjointView.h @@ -68,7 +68,7 @@ template<typename MatrixType, unsigned int UpLo> class SelfAdjointView enum { Mode = ei_traits<SelfAdjointView>::Mode }; - typedef typename MatrixType::PlainMatrixType PlainMatrixType; + typedef typename MatrixType::PlainObject PlainObject; inline SelfAdjointView(const MatrixType& matrix) : m_matrix(matrix) { ei_assert(ei_are_flags_consistent<Mode>::ret); } @@ -147,8 +147,8 @@ template<typename MatrixType, unsigned int UpLo> class SelfAdjointView /////////// Cholesky module /////////// - const LLT<PlainMatrixType, UpLo> llt() const; - const LDLT<PlainMatrixType> ldlt() const; + const LLT<PlainObject, UpLo> llt() const; + const LDLT<PlainObject> ldlt() const; protected: const typename MatrixType::Nested m_matrix; diff --git a/Eigen/src/Core/SelfCwiseBinaryOp.h b/Eigen/src/Core/SelfCwiseBinaryOp.h index 58aee182d..529a9994d 100644 --- a/Eigen/src/Core/SelfCwiseBinaryOp.h +++ b/Eigen/src/Core/SelfCwiseBinaryOp.h @@ -125,8 +125,8 @@ template<typename Derived> inline Derived& DenseBase<Derived>::operator*=(const Scalar& other) { SelfCwiseBinaryOp<ei_scalar_product_op<Scalar>, Derived> tmp(derived()); - typedef typename Derived::PlainMatrixType PlainMatrixType; - tmp = PlainMatrixType::Constant(rows(),cols(),other); + typedef typename Derived::PlainObject PlainObject; + tmp = PlainObject::Constant(rows(),cols(),other); return derived(); } @@ -134,8 +134,8 @@ template<typename Derived> inline Derived& DenseBase<Derived>::operator/=(const Scalar& other) { SelfCwiseBinaryOp<typename ei_meta_if<NumTraits<Scalar>::HasFloatingPoint,ei_scalar_product_op<Scalar>,ei_scalar_quotient_op<Scalar> >::ret, Derived> tmp(derived()); - typedef typename Derived::PlainMatrixType PlainMatrixType; - tmp = PlainMatrixType::Constant(rows(),cols(), NumTraits<Scalar>::HasFloatingPoint ? Scalar(1)/other : other); + typedef typename Derived::PlainObject PlainObject; + tmp = PlainObject::Constant(rows(),cols(), NumTraits<Scalar>::HasFloatingPoint ? Scalar(1)/other : other); return derived(); } diff --git a/Eigen/src/Core/Transpose.h b/Eigen/src/Core/Transpose.h index 47dae5776..1f064d1c2 100644 --- a/Eigen/src/Core/Transpose.h +++ b/Eigen/src/Core/Transpose.h @@ -296,25 +296,6 @@ struct ei_blas_traits<SelfCwiseBinaryOp<BinOp,NestedXpr> > static inline const XprType extract(const XprType& x) { return x; } }; - -template<typename T, int Access=ei_blas_traits<T>::ActualAccess> -struct ei_extract_data_selector { - static typename T::Scalar* run(const T& m) - { - return &ei_blas_traits<T>::extract(m).const_cast_derived().coeffRef(0,0); - } -}; - -template<typename T> -struct ei_extract_data_selector<T,NoDirectAccess> { - static typename T::Scalar* run(const T&) { return 0; } -}; - -template<typename T> typename T::Scalar* ei_extract_data(const T& m) -{ - return ei_extract_data_selector<T>::run(m); -} - template<typename Scalar, bool DestIsTranposed, typename OtherDerived> struct ei_check_transpose_aliasing_selector { diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h index c61a6d7cc..2230680d1 100644 --- a/Eigen/src/Core/TriangularMatrix.h +++ b/Eigen/src/Core/TriangularMatrix.h @@ -32,7 +32,7 @@ * * \brief Base class for triangular part in a matrix */ -template<typename Derived> class TriangularBase : public AnyMatrixBase<Derived> +template<typename Derived> class TriangularBase : public EigenBase<Derived> { public: @@ -149,7 +149,7 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView typedef TriangularBase<TriangularView> Base; typedef typename ei_traits<TriangularView>::Scalar Scalar; typedef _MatrixType MatrixType; - typedef typename MatrixType::PlainMatrixType DenseMatrixType; + typedef typename MatrixType::PlainObject DenseMatrixType; typedef typename MatrixType::Nested MatrixTypeNested; typedef typename ei_cleantype<MatrixTypeNested>::type _MatrixTypeNested; diff --git a/Eigen/src/Core/arch/SSE/PacketMath.h b/Eigen/src/Core/arch/SSE/PacketMath.h index a5a56f759..f78bf0dd3 100644 --- a/Eigen/src/Core/arch/SSE/PacketMath.h +++ b/Eigen/src/Core/arch/SSE/PacketMath.h @@ -122,7 +122,7 @@ template<> EIGEN_STRONG_INLINE Packet4f ei_pmul<Packet4f>(const Packet4f& a, con template<> EIGEN_STRONG_INLINE Packet2d ei_pmul<Packet2d>(const Packet2d& a, const Packet2d& b) { return _mm_mul_pd(a,b); } template<> EIGEN_STRONG_INLINE Packet4i ei_pmul<Packet4i>(const Packet4i& a, const Packet4i& b) { -#ifdef __SSE4_1__ +#ifdef EIGEN_VECTORIZE_SSE4_1 return _mm_mullo_epi32(a,b); #else // this version is slightly faster than 4 scalar products @@ -269,7 +269,7 @@ template<> EIGEN_STRONG_INLINE Packet2d ei_pabs(const Packet2d& a) } template<> EIGEN_STRONG_INLINE Packet4i ei_pabs(const Packet4i& a) { - #ifdef __SSSE3__ + #ifdef EIGEN_VECTORIZE_SSSE3 return _mm_abs_epi32(a); #else Packet4i aux = _mm_srai_epi32(a,31); @@ -278,7 +278,7 @@ template<> EIGEN_STRONG_INLINE Packet4i ei_pabs(const Packet4i& a) } -#ifdef __SSE3__ +#ifdef EIGEN_VECTORIZE_SSE3 // TODO implement SSE2 versions as well as integer versions template<> EIGEN_STRONG_INLINE Packet4f ei_preduxp<Packet4f>(const Packet4f* vecs) { @@ -439,7 +439,7 @@ template<> EIGEN_STRONG_INLINE int ei_predux_max<Packet4i>(const Packet4i& a) // } #endif -#ifdef __SSSE3__ +#ifdef EIGEN_VECTORIZE_SSSE3 // SSSE3 versions template<int Offset> struct ei_palign_impl<Offset,Packet4f> diff --git a/Eigen/src/Core/products/CoeffBasedProduct.h b/Eigen/src/Core/products/CoeffBasedProduct.h index f030d59b5..3343b1875 100644 --- a/Eigen/src/Core/products/CoeffBasedProduct.h +++ b/Eigen/src/Core/products/CoeffBasedProduct.h @@ -109,7 +109,7 @@ class CoeffBasedProduct typedef MatrixBase<CoeffBasedProduct> Base; EIGEN_DENSE_PUBLIC_INTERFACE(CoeffBasedProduct) - typedef typename Base::PlainMatrixType PlainMatrixType; + typedef typename Base::PlainObject PlainObject; private: @@ -181,8 +181,8 @@ class CoeffBasedProduct return res; } - // Implicit convertion to the nested type (trigger the evaluation of the product) - operator const PlainMatrixType& () const + // Implicit conversion to the nested type (trigger the evaluation of the product) + operator const PlainObject& () const { m_result.lazyAssign(*this); return m_result; @@ -205,15 +205,15 @@ class CoeffBasedProduct const LhsNested m_lhs; const RhsNested m_rhs; - mutable PlainMatrixType m_result; + mutable PlainObject m_result; }; // here we need to overload the nested rule for products // such that the nested type is a const reference to a plain matrix -template<typename Lhs, typename Rhs, int N, typename PlainMatrixType> -struct ei_nested<CoeffBasedProduct<Lhs,Rhs,EvalBeforeNestingBit|EvalBeforeAssigningBit>, N, PlainMatrixType> +template<typename Lhs, typename Rhs, int N, typename PlainObject> +struct ei_nested<CoeffBasedProduct<Lhs,Rhs,EvalBeforeNestingBit|EvalBeforeAssigningBit>, N, PlainObject> { - typedef PlainMatrixType const& type; + typedef PlainObject const& type; }; /*************************************************************************** diff --git a/Eigen/src/Core/products/GeneralBlockPanelKernel.h b/Eigen/src/Core/products/GeneralBlockPanelKernel.h index fe1987bdd..18e913b0e 100644 --- a/Eigen/src/Core/products/GeneralBlockPanelKernel.h +++ b/Eigen/src/Core/products/GeneralBlockPanelKernel.h @@ -27,6 +27,12 @@ #ifndef EIGEN_EXTERN_INSTANTIATIONS +#ifdef EIGEN_HAS_FUSE_CJMADD +#define CJMADD(A,B,C,T) C = cj.pmadd(A,B,C); +#else +#define CJMADD(A,B,C,T) T = A; T = cj.pmul(T,B); C = ei_padd(C,T); +#endif + // optimized GEneral packed Block * packed Panel product kernel template<typename Scalar, int mr, int nr, typename Conj> struct ei_gebp_kernel @@ -74,65 +80,111 @@ struct ei_gebp_kernel const Scalar* blB = &blockB[j2*strideB*PacketSize+offsetB*nr]; for(int k=0; k<peeled_kc; k+=4) { - PacketType B0, B1, B2, B3, A0, A1; - - A0 = ei_pload(&blA[0*PacketSize]); - A1 = ei_pload(&blA[1*PacketSize]); - B0 = ei_pload(&blB[0*PacketSize]); - B1 = ei_pload(&blB[1*PacketSize]); - C0 = cj.pmadd(A0, B0, C0); - if(nr==4) B2 = ei_pload(&blB[2*PacketSize]); - C4 = cj.pmadd(A1, B0, C4); - if(nr==4) B3 = ei_pload(&blB[3*PacketSize]); - B0 = ei_pload(&blB[(nr==4 ? 4 : 2)*PacketSize]); - C1 = cj.pmadd(A0, B1, C1); - C5 = cj.pmadd(A1, B1, C5); - B1 = ei_pload(&blB[(nr==4 ? 5 : 3)*PacketSize]); - if(nr==4) C2 = cj.pmadd(A0, B2, C2); - if(nr==4) C6 = cj.pmadd(A1, B2, C6); - if(nr==4) B2 = ei_pload(&blB[6*PacketSize]); - if(nr==4) C3 = cj.pmadd(A0, B3, C3); - A0 = ei_pload(&blA[2*PacketSize]); - if(nr==4) C7 = cj.pmadd(A1, B3, C7); - A1 = ei_pload(&blA[3*PacketSize]); - if(nr==4) B3 = ei_pload(&blB[7*PacketSize]); - C0 = cj.pmadd(A0, B0, C0); - C4 = cj.pmadd(A1, B0, C4); - B0 = ei_pload(&blB[(nr==4 ? 8 : 4)*PacketSize]); - C1 = cj.pmadd(A0, B1, C1); - C5 = cj.pmadd(A1, B1, C5); - B1 = ei_pload(&blB[(nr==4 ? 9 : 5)*PacketSize]); - if(nr==4) C2 = cj.pmadd(A0, B2, C2); - if(nr==4) C6 = cj.pmadd(A1, B2, C6); - if(nr==4) B2 = ei_pload(&blB[10*PacketSize]); - if(nr==4) C3 = cj.pmadd(A0, B3, C3); - A0 = ei_pload(&blA[4*PacketSize]); - if(nr==4) C7 = cj.pmadd(A1, B3, C7); - A1 = ei_pload(&blA[5*PacketSize]); - if(nr==4) B3 = ei_pload(&blB[11*PacketSize]); - - C0 = cj.pmadd(A0, B0, C0); - C4 = cj.pmadd(A1, B0, C4); - B0 = ei_pload(&blB[(nr==4 ? 12 : 6)*PacketSize]); - C1 = cj.pmadd(A0, B1, C1); - C5 = cj.pmadd(A1, B1, C5); - B1 = ei_pload(&blB[(nr==4 ? 13 : 7)*PacketSize]); - if(nr==4) C2 = cj.pmadd(A0, B2, C2); - if(nr==4) C6 = cj.pmadd(A1, B2, C6); - if(nr==4) B2 = ei_pload(&blB[14*PacketSize]); - if(nr==4) C3 = cj.pmadd(A0, B3, C3); - A0 = ei_pload(&blA[6*PacketSize]); - if(nr==4) C7 = cj.pmadd(A1, B3, C7); - A1 = ei_pload(&blA[7*PacketSize]); - if(nr==4) B3 = ei_pload(&blB[15*PacketSize]); - C0 = cj.pmadd(A0, B0, C0); - C4 = cj.pmadd(A1, B0, C4); - C1 = cj.pmadd(A0, B1, C1); - C5 = cj.pmadd(A1, B1, C5); - if(nr==4) C2 = cj.pmadd(A0, B2, C2); - if(nr==4) C6 = cj.pmadd(A1, B2, C6); - if(nr==4) C3 = cj.pmadd(A0, B3, C3); - if(nr==4) C7 = cj.pmadd(A1, B3, C7); + if(nr==2) + { + PacketType B0, T0, A0, A1; + + A0 = ei_pload(&blA[0*PacketSize]); + A1 = ei_pload(&blA[1*PacketSize]); + B0 = ei_pload(&blB[0*PacketSize]); + CJMADD(A0,B0,C0,T0); + CJMADD(A1,B0,C4,T0); + B0 = ei_pload(&blB[1*PacketSize]); + CJMADD(A0,B0,C1,T0); + CJMADD(A1,B0,C5,T0); + + A0 = ei_pload(&blA[2*PacketSize]); + A1 = ei_pload(&blA[3*PacketSize]); + B0 = ei_pload(&blB[2*PacketSize]); + CJMADD(A0,B0,C0,T0); + CJMADD(A1,B0,C4,T0); + B0 = ei_pload(&blB[3*PacketSize]); + CJMADD(A0,B0,C1,T0); + CJMADD(A1,B0,C5,T0); + + A0 = ei_pload(&blA[4*PacketSize]); + A1 = ei_pload(&blA[5*PacketSize]); + B0 = ei_pload(&blB[4*PacketSize]); + CJMADD(A0,B0,C0,T0); + CJMADD(A1,B0,C4,T0); + B0 = ei_pload(&blB[5*PacketSize]); + CJMADD(A0,B0,C1,T0); + CJMADD(A1,B0,C5,T0); + + A0 = ei_pload(&blA[6*PacketSize]); + A1 = ei_pload(&blA[7*PacketSize]); + B0 = ei_pload(&blB[6*PacketSize]); + CJMADD(A0,B0,C0,T0); + CJMADD(A1,B0,C4,T0); + B0 = ei_pload(&blB[7*PacketSize]); + CJMADD(A0,B0,C1,T0); + CJMADD(A1,B0,C5,T0); + } + else + { + + PacketType B0, B1, B2, B3, A0, A1; + PacketType T0, T1; + + A0 = ei_pload(&blA[0*PacketSize]); + A1 = ei_pload(&blA[1*PacketSize]); + B0 = ei_pload(&blB[0*PacketSize]); + B1 = ei_pload(&blB[1*PacketSize]); + + CJMADD(A0,B0,C0,T0); + if(nr==4) B2 = ei_pload(&blB[2*PacketSize]); + CJMADD(A1,B0,C4,T1); + if(nr==4) B3 = ei_pload(&blB[3*PacketSize]); + B0 = ei_pload(&blB[(nr==4 ? 4 : 2)*PacketSize]); + CJMADD(A0,B1,C1,T0); + CJMADD(A1,B1,C5,T1); + B1 = ei_pload(&blB[(nr==4 ? 5 : 3)*PacketSize]); + if(nr==4) { CJMADD(A0,B2,C2,T0); } + if(nr==4) { CJMADD(A1,B2,C6,T1); } + if(nr==4) B2 = ei_pload(&blB[6*PacketSize]); + if(nr==4) { CJMADD(A0,B3,C3,T0); } + A0 = ei_pload(&blA[2*PacketSize]); + if(nr==4) { CJMADD(A1,B3,C7,T1); } + A1 = ei_pload(&blA[3*PacketSize]); + if(nr==4) B3 = ei_pload(&blB[7*PacketSize]); + CJMADD(A0,B0,C0,T0); + CJMADD(A1,B0,C4,T1); + B0 = ei_pload(&blB[(nr==4 ? 8 : 4)*PacketSize]); + CJMADD(A0,B1,C1,T0); + CJMADD(A1,B1,C5,T1); + B1 = ei_pload(&blB[(nr==4 ? 9 : 5)*PacketSize]); + if(nr==4) { CJMADD(A0,B2,C2,T0); } + if(nr==4) { CJMADD(A1,B2,C6,T1); } + if(nr==4) B2 = ei_pload(&blB[10*PacketSize]); + if(nr==4) { CJMADD(A0,B3,C3,T0); } + A0 = ei_pload(&blA[4*PacketSize]); + if(nr==4) { CJMADD(A1,B3,C7,T1); } + A1 = ei_pload(&blA[5*PacketSize]); + if(nr==4) B3 = ei_pload(&blB[11*PacketSize]); + + CJMADD(A0,B0,C0,T0); + CJMADD(A1,B0,C4,T1); + B0 = ei_pload(&blB[(nr==4 ? 12 : 6)*PacketSize]); + CJMADD(A0,B1,C1,T0); + CJMADD(A1,B1,C5,T1); + B1 = ei_pload(&blB[(nr==4 ? 13 : 7)*PacketSize]); + if(nr==4) { CJMADD(A0,B2,C2,T0); } + if(nr==4) { CJMADD(A1,B2,C6,T1); } + if(nr==4) B2 = ei_pload(&blB[14*PacketSize]); + if(nr==4) { CJMADD(A0,B3,C3,T0); } + A0 = ei_pload(&blA[6*PacketSize]); + if(nr==4) { CJMADD(A1,B3,C7,T1); } + A1 = ei_pload(&blA[7*PacketSize]); + if(nr==4) B3 = ei_pload(&blB[15*PacketSize]); + CJMADD(A0,B0,C0,T0); + CJMADD(A1,B0,C4,T1); + CJMADD(A0,B1,C1,T0); + CJMADD(A1,B1,C5,T1); + if(nr==4) { CJMADD(A0,B2,C2,T0); } + if(nr==4) { CJMADD(A1,B2,C6,T1); } + if(nr==4) { CJMADD(A0,B3,C3,T0); } + if(nr==4) { CJMADD(A1,B3,C7,T1); } + } blB += 4*nr*PacketSize; blA += 4*mr; @@ -140,22 +192,40 @@ struct ei_gebp_kernel // process remaining peeled loop for(int k=peeled_kc; k<depth; k++) { - PacketType B0, B1, B2, B3, A0, A1; - - A0 = ei_pload(&blA[0*PacketSize]); - A1 = ei_pload(&blA[1*PacketSize]); - B0 = ei_pload(&blB[0*PacketSize]); - B1 = ei_pload(&blB[1*PacketSize]); - C0 = cj.pmadd(A0, B0, C0); - if(nr==4) B2 = ei_pload(&blB[2*PacketSize]); - C4 = cj.pmadd(A1, B0, C4); - if(nr==4) B3 = ei_pload(&blB[3*PacketSize]); - C1 = cj.pmadd(A0, B1, C1); - C5 = cj.pmadd(A1, B1, C5); - if(nr==4) C2 = cj.pmadd(A0, B2, C2); - if(nr==4) C6 = cj.pmadd(A1, B2, C6); - if(nr==4) C3 = cj.pmadd(A0, B3, C3); - if(nr==4) C7 = cj.pmadd(A1, B3, C7); + if(nr==2) + { + PacketType B0, T0, A0, A1; + + A0 = ei_pload(&blA[0*PacketSize]); + A1 = ei_pload(&blA[1*PacketSize]); + B0 = ei_pload(&blB[0*PacketSize]); + CJMADD(A0,B0,C0,T0); + CJMADD(A1,B0,C4,T0); + B0 = ei_pload(&blB[1*PacketSize]); + CJMADD(A0,B0,C1,T0); + CJMADD(A1,B0,C5,T0); + } + else + { + PacketType B0, B1, B2, B3, A0, A1, T0, T1; + + A0 = ei_pload(&blA[0*PacketSize]); + A1 = ei_pload(&blA[1*PacketSize]); + B0 = ei_pload(&blB[0*PacketSize]); + B1 = ei_pload(&blB[1*PacketSize]); + + CJMADD(A0,B0,C0,T0); + if(nr==4) B2 = ei_pload(&blB[2*PacketSize]); + CJMADD(A1,B0,C4,T1); + if(nr==4) B3 = ei_pload(&blB[3*PacketSize]); + B0 = ei_pload(&blB[(nr==4 ? 4 : 2)*PacketSize]); + CJMADD(A0,B1,C1,T0); + CJMADD(A1,B1,C5,T1); + if(nr==4) { CJMADD(A0,B2,C2,T0); } + if(nr==4) { CJMADD(A1,B2,C6,T1); } + if(nr==4) { CJMADD(A0,B3,C3,T0); } + if(nr==4) { CJMADD(A1,B3,C7,T1); } + } blB += nr*PacketSize; blA += mr; @@ -189,45 +259,79 @@ struct ei_gebp_kernel const Scalar* blB = &blockB[j2*strideB*PacketSize+offsetB*nr]; for(int k=0; k<peeled_kc; k+=4) { - PacketType B0, B1, B2, B3, A0; - - A0 = ei_pload(&blA[0*PacketSize]); - B0 = ei_pload(&blB[0*PacketSize]); - B1 = ei_pload(&blB[1*PacketSize]); - C0 = cj.pmadd(A0, B0, C0); - if(nr==4) B2 = ei_pload(&blB[2*PacketSize]); - if(nr==4) B3 = ei_pload(&blB[3*PacketSize]); - B0 = ei_pload(&blB[(nr==4 ? 4 : 2)*PacketSize]); - C1 = cj.pmadd(A0, B1, C1); - B1 = ei_pload(&blB[(nr==4 ? 5 : 3)*PacketSize]); - if(nr==4) C2 = cj.pmadd(A0, B2, C2); - if(nr==4) B2 = ei_pload(&blB[6*PacketSize]); - if(nr==4) C3 = cj.pmadd(A0, B3, C3); - A0 = ei_pload(&blA[1*PacketSize]); - if(nr==4) B3 = ei_pload(&blB[7*PacketSize]); - C0 = cj.pmadd(A0, B0, C0); - B0 = ei_pload(&blB[(nr==4 ? 8 : 4)*PacketSize]); - C1 = cj.pmadd(A0, B1, C1); - B1 = ei_pload(&blB[(nr==4 ? 9 : 5)*PacketSize]); - if(nr==4) C2 = cj.pmadd(A0, B2, C2); - if(nr==4) B2 = ei_pload(&blB[10*PacketSize]); - if(nr==4) C3 = cj.pmadd(A0, B3, C3); - A0 = ei_pload(&blA[2*PacketSize]); - if(nr==4) B3 = ei_pload(&blB[11*PacketSize]); - - C0 = cj.pmadd(A0, B0, C0); - B0 = ei_pload(&blB[(nr==4 ? 12 : 6)*PacketSize]); - C1 = cj.pmadd(A0, B1, C1); - B1 = ei_pload(&blB[(nr==4 ? 13 : 7)*PacketSize]); - if(nr==4) C2 = cj.pmadd(A0, B2, C2); - if(nr==4) B2 = ei_pload(&blB[14*PacketSize]); - if(nr==4) C3 = cj.pmadd(A0, B3, C3); - A0 = ei_pload(&blA[3*PacketSize]); - if(nr==4) B3 = ei_pload(&blB[15*PacketSize]); - C0 = cj.pmadd(A0, B0, C0); - C1 = cj.pmadd(A0, B1, C1); - if(nr==4) C2 = cj.pmadd(A0, B2, C2); - if(nr==4) C3 = cj.pmadd(A0, B3, C3); + if(nr==2) + { + PacketType B0, T0, A0; + + A0 = ei_pload(&blA[0*PacketSize]); + B0 = ei_pload(&blB[0*PacketSize]); + CJMADD(A0,B0,C0,T0); + B0 = ei_pload(&blB[1*PacketSize]); + CJMADD(A0,B0,C1,T0); + + A0 = ei_pload(&blA[1*PacketSize]); + B0 = ei_pload(&blB[2*PacketSize]); + CJMADD(A0,B0,C0,T0); + B0 = ei_pload(&blB[3*PacketSize]); + CJMADD(A0,B0,C1,T0); + + A0 = ei_pload(&blA[2*PacketSize]); + B0 = ei_pload(&blB[4*PacketSize]); + CJMADD(A0,B0,C0,T0); + B0 = ei_pload(&blB[5*PacketSize]); + CJMADD(A0,B0,C1,T0); + + A0 = ei_pload(&blA[3*PacketSize]); + B0 = ei_pload(&blB[6*PacketSize]); + CJMADD(A0,B0,C0,T0); + B0 = ei_pload(&blB[7*PacketSize]); + CJMADD(A0,B0,C1,T0); + } + else + { + + PacketType B0, B1, B2, B3, A0; + PacketType T0, T1; + + A0 = ei_pload(&blA[0*PacketSize]); + B0 = ei_pload(&blB[0*PacketSize]); + B1 = ei_pload(&blB[1*PacketSize]); + + CJMADD(A0,B0,C0,T0); + if(nr==4) B2 = ei_pload(&blB[2*PacketSize]); + if(nr==4) B3 = ei_pload(&blB[3*PacketSize]); + B0 = ei_pload(&blB[(nr==4 ? 4 : 2)*PacketSize]); + CJMADD(A0,B1,C1,T1); + B1 = ei_pload(&blB[(nr==4 ? 5 : 3)*PacketSize]); + if(nr==4) { CJMADD(A0,B2,C2,T0); } + if(nr==4) B2 = ei_pload(&blB[6*PacketSize]); + if(nr==4) { CJMADD(A0,B3,C3,T1); } + A0 = ei_pload(&blA[1*PacketSize]); + if(nr==4) B3 = ei_pload(&blB[7*PacketSize]); + CJMADD(A0,B0,C0,T0); + B0 = ei_pload(&blB[(nr==4 ? 8 : 4)*PacketSize]); + CJMADD(A0,B1,C1,T1); + B1 = ei_pload(&blB[(nr==4 ? 9 : 5)*PacketSize]); + if(nr==4) { CJMADD(A0,B2,C2,T0); } + if(nr==4) B2 = ei_pload(&blB[10*PacketSize]); + if(nr==4) { CJMADD(A0,B3,C3,T1); } + A0 = ei_pload(&blA[2*PacketSize]); + if(nr==4) B3 = ei_pload(&blB[11*PacketSize]); + + CJMADD(A0,B0,C0,T0); + B0 = ei_pload(&blB[(nr==4 ? 12 : 6)*PacketSize]); + CJMADD(A0,B1,C1,T1); + B1 = ei_pload(&blB[(nr==4 ? 13 : 7)*PacketSize]); + if(nr==4) { CJMADD(A0,B2,C2,T0); } + if(nr==4) B2 = ei_pload(&blB[14*PacketSize]); + if(nr==4) { CJMADD(A0,B3,C3,T1); } + A0 = ei_pload(&blA[3*PacketSize]); + if(nr==4) B3 = ei_pload(&blB[15*PacketSize]); + CJMADD(A0,B0,C0,T0); + CJMADD(A0,B1,C1,T1); + if(nr==4) { CJMADD(A0,B2,C2,T0); } + if(nr==4) { CJMADD(A0,B3,C3,T1); } + } blB += 4*nr*PacketSize; blA += 4*PacketSize; @@ -235,17 +339,32 @@ struct ei_gebp_kernel // process remaining peeled loop for(int k=peeled_kc; k<depth; k++) { - PacketType B0, B1, B2, B3, A0; - - A0 = ei_pload(&blA[0*PacketSize]); - B0 = ei_pload(&blB[0*PacketSize]); - B1 = ei_pload(&blB[1*PacketSize]); - C0 = cj.pmadd(A0, B0, C0); - if(nr==4) B2 = ei_pload(&blB[2*PacketSize]); - if(nr==4) B3 = ei_pload(&blB[3*PacketSize]); - C1 = cj.pmadd(A0, B1, C1); - if(nr==4) C2 = cj.pmadd(A0, B2, C2); - if(nr==4) C3 = cj.pmadd(A0, B3, C3); + if(nr==2) + { + PacketType B0, T0, A0; + + A0 = ei_pload(&blA[0*PacketSize]); + B0 = ei_pload(&blB[0*PacketSize]); + CJMADD(A0,B0,C0,T0); + B0 = ei_pload(&blB[1*PacketSize]); + CJMADD(A0,B0,C1,T0); + } + else + { + PacketType B0, B1, B2, B3, A0; + PacketType T0, T1; + + A0 = ei_pload(&blA[0*PacketSize]); + B0 = ei_pload(&blB[0*PacketSize]); + B1 = ei_pload(&blB[1*PacketSize]); + if(nr==4) B2 = ei_pload(&blB[2*PacketSize]); + if(nr==4) B3 = ei_pload(&blB[3*PacketSize]); + + CJMADD(A0,B0,C0,T0); + CJMADD(A0,B1,C1,T1); + if(nr==4) { CJMADD(A0,B2,C2,T0); } + if(nr==4) { CJMADD(A0,B3,C3,T1); } + } blB += nr*PacketSize; blA += PacketSize; @@ -268,17 +387,32 @@ struct ei_gebp_kernel const Scalar* blB = &blockB[j2*strideB*PacketSize+offsetB*nr]; for(int k=0; k<depth; k++) { - Scalar B0, B1, B2, B3, A0; - - A0 = blA[k]; - B0 = blB[0*PacketSize]; - B1 = blB[1*PacketSize]; - C0 = cj.pmadd(A0, B0, C0); - if(nr==4) B2 = blB[2*PacketSize]; - if(nr==4) B3 = blB[3*PacketSize]; - C1 = cj.pmadd(A0, B1, C1); - if(nr==4) C2 = cj.pmadd(A0, B2, C2); - if(nr==4) C3 = cj.pmadd(A0, B3, C3); + if(nr==2) + { + Scalar B0, T0, A0; + + A0 = blA[0*PacketSize]; + B0 = blB[0*PacketSize]; + CJMADD(A0,B0,C0,T0); + B0 = blB[1*PacketSize]; + CJMADD(A0,B0,C1,T0); + } + else + { + Scalar B0, B1, B2, B3, A0; + Scalar T0, T1; + + A0 = blA[k]; + B0 = blB[0*PacketSize]; + B1 = blB[1*PacketSize]; + if(nr==4) B2 = blB[2*PacketSize]; + if(nr==4) B3 = blB[3*PacketSize]; + + CJMADD(A0,B0,C0,T0); + CJMADD(A0,B1,C1,T1); + if(nr==4) { CJMADD(A0,B2,C2,T0); } + if(nr==4) { CJMADD(A0,B3,C3,T1); } + } blB += nr*PacketSize; } @@ -310,13 +444,13 @@ struct ei_gebp_kernel const Scalar* blB = &blockB[j2*strideB*PacketSize+offsetB]; for(int k=0; k<depth; k++) { - PacketType B0, A0, A1; + PacketType B0, A0, A1, T0, T1; A0 = ei_pload(&blA[0*PacketSize]); A1 = ei_pload(&blA[1*PacketSize]); B0 = ei_pload(&blB[0*PacketSize]); - C0 = cj.pmadd(A0, B0, C0); - C4 = cj.pmadd(A1, B0, C4); + CJMADD(A0,B0,C0,T0); + CJMADD(A1,B0,C4,T1); blB += PacketSize; blA += mr; @@ -334,7 +468,7 @@ struct ei_gebp_kernel #endif PacketType C0 = ei_ploadu(&res[(j2+0)*resStride + i]); - + const Scalar* blB = &blockB[j2*strideB*PacketSize+offsetB]; for(int k=0; k<depth; k++) { @@ -363,6 +497,8 @@ struct ei_gebp_kernel } }; +#undef CJMADD + // pack a block of the lhs // The travesal is as follow (mr==4): // 0 4 8 12 ... @@ -474,7 +610,7 @@ struct ei_gemm_pack_rhs<Scalar, nr, ColMajor, PanelMode> // skip what we have after if(PanelMode) count += PacketSize * nr * (stride-offset-depth); } - + // copy the remaining columns one at a time (nr==1) for(int j2=packet_cols; j2<cols; ++j2) { diff --git a/Eigen/src/Core/util/BlasUtil.h b/Eigen/src/Core/util/BlasUtil.h index 3777464dc..4d216d77a 100644 --- a/Eigen/src/Core/util/BlasUtil.h +++ b/Eigen/src/Core/util/BlasUtil.h @@ -166,7 +166,7 @@ template<typename XprType> struct ei_blas_traits }; typedef typename ei_meta_if<int(ActualAccess)==HasDirectAccess, ExtractType, - typename _ExtractType::PlainMatrixType + typename _ExtractType::PlainObject >::ret DirectLinearAccessType; static inline ExtractType extract(const XprType& x) { return x; } static inline Scalar extractScalarFactor(const XprType&) { return Scalar(1); } @@ -227,7 +227,7 @@ struct ei_blas_traits<Transpose<NestedXpr> > typedef Transpose<typename Base::_ExtractType> _ExtractType; typedef typename ei_meta_if<int(Base::ActualAccess)==HasDirectAccess, ExtractType, - typename ExtractType::PlainMatrixType + typename ExtractType::PlainObject >::ret DirectLinearAccessType; enum { IsTransposed = Base::IsTransposed ? 0 : 1 @@ -236,4 +236,22 @@ struct ei_blas_traits<Transpose<NestedXpr> > static inline Scalar extractScalarFactor(const XprType& x) { return Base::extractScalarFactor(x.nestedExpression()); } }; +template<typename T, int Access=ei_blas_traits<T>::ActualAccess> +struct ei_extract_data_selector { + static const typename T::Scalar* run(const T& m) + { + return &ei_blas_traits<T>::extract(m).const_cast_derived().coeffRef(0,0); // FIXME this should be .data() + } +}; + +template<typename T> +struct ei_extract_data_selector<T,NoDirectAccess> { + static typename T::Scalar* run(const T&) { return 0; } +}; + +template<typename T> const typename T::Scalar* ei_extract_data(const T& m) +{ + return ei_extract_data_selector<T>::run(m); +} + #endif // EIGEN_BLASUTIL_H diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index c2d45dc30..6096272fa 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -29,7 +29,7 @@ template<typename T> struct ei_traits; template<typename T> struct NumTraits; -template<typename Derived> struct AnyMatrixBase; +template<typename Derived> struct EigenBase; template<typename _Scalar, int _Rows, int _Cols, int _Options = EIGEN_DEFAULT_MATRIX_STORAGE_ORDER_OPTION | AutoAlign, diff --git a/Eigen/src/Core/util/Macros.h b/Eigen/src/Core/util/Macros.h index dc1aa150b..37ccef047 100644 --- a/Eigen/src/Core/util/Macros.h +++ b/Eigen/src/Core/util/Macros.h @@ -211,7 +211,7 @@ using Eigen::ei_cos; */ #if !EIGEN_ALIGN #define EIGEN_ALIGN_TO_BOUNDARY(n) -#elif (defined __GNUC__) +#elif (defined __GNUC__) || (defined __PGI) #define EIGEN_ALIGN_TO_BOUNDARY(n) __attribute__((aligned(n))) #elif (defined _MSC_VER) #define EIGEN_ALIGN_TO_BOUNDARY(n) __declspec(align(n)) diff --git a/Eigen/src/Core/util/XprHelper.h b/Eigen/src/Core/util/XprHelper.h index 8ddf4450a..eceb5ab2a 100644 --- a/Eigen/src/Core/util/XprHelper.h +++ b/Eigen/src/Core/util/XprHelper.h @@ -147,7 +147,7 @@ template<typename T, typename StorageType = typename ei_traits<T>::StorageType> template<typename T> struct ei_eval<T,Dense> { typedef typename ei_plain_matrix_type<T>::type type; -// typedef typename T::PlainMatrixType type; +// typedef typename T::PlainObject type; // typedef T::Matrix<typename ei_traits<T>::Scalar, // ei_traits<T>::RowsAtCompileTime, // ei_traits<T>::ColsAtCompileTime, @@ -201,6 +201,18 @@ template<typename T> struct ei_plain_matrix_type_row_major // we should be able to get rid of this one too template<typename T> struct ei_must_nest_by_value { enum { ret = false }; }; +template<class T> +struct ei_is_reference +{ + enum { ret = false }; +}; + +template<class T> +struct ei_is_reference<T&> +{ + enum { ret = true }; +}; + /** * The reference selector for template expressions. The idea is that we don't * need to use references for expressions since they are light weight proxy @@ -234,7 +246,7 @@ struct ei_ref_selector * const Matrix3d&, because the internal logic of ei_nested determined that since a was already a matrix, there was no point * in copying it into another matrix. */ -template<typename T, int n=1, typename PlainMatrixType = typename ei_eval<T>::type> struct ei_nested +template<typename T, int n=1, typename PlainObject = typename ei_eval<T>::type> struct ei_nested { enum { CostEval = (n+1) * int(NumTraits<typename ei_traits<T>::Scalar>::ReadCost), @@ -244,7 +256,7 @@ template<typename T, int n=1, typename PlainMatrixType = typename ei_eval<T>::ty typedef typename ei_meta_if< ( int(ei_traits<T>::Flags) & EvalBeforeNestingBit ) || ( int(CostEval) <= int(CostNoEval) ), - PlainMatrixType, + PlainObject, typename ei_ref_selector<T>::type >::ret type; }; @@ -258,7 +270,7 @@ template<unsigned int Flags> struct ei_are_flags_consistent * overloads for complex types */ template<typename Derived,typename Scalar,typename OtherScalar, bool EnableIt = !ei_is_same_type<Scalar,OtherScalar>::ret > -struct ei_special_scalar_op_base : public AnyMatrixBase<Derived> +struct ei_special_scalar_op_base : public EigenBase<Derived> { // dummy operator* so that the // "using ei_special_scalar_op_base::operator*" compiles @@ -266,7 +278,7 @@ struct ei_special_scalar_op_base : public AnyMatrixBase<Derived> }; template<typename Derived,typename Scalar,typename OtherScalar> -struct ei_special_scalar_op_base<Derived,Scalar,OtherScalar,true> : public AnyMatrixBase<Derived> +struct ei_special_scalar_op_base<Derived,Scalar,OtherScalar,true> : public EigenBase<Derived> { const CwiseUnaryOp<ei_scalar_multiple2_op<Scalar,OtherScalar>, Derived> operator*(const OtherScalar& scalar) const diff --git a/Eigen/src/Eigen2Support/TriangularSolver.h b/Eigen/src/Eigen2Support/TriangularSolver.h index 94b92577e..833dd58d2 100644 --- a/Eigen/src/Eigen2Support/TriangularSolver.h +++ b/Eigen/src/Eigen2Support/TriangularSolver.h @@ -37,7 +37,7 @@ const unsigned int UnitLowerTriangular = UnitLower; template<typename ExpressionType, unsigned int Added, unsigned int Removed> template<typename OtherDerived> -typename ExpressionType::PlainMatrixType +typename ExpressionType::PlainObject Flagged<ExpressionType,Added,Removed>::solveTriangular(const MatrixBase<OtherDerived>& other) const { return m_matrix.template triangularView<Added>.solve(other.derived()); diff --git a/Eigen/src/Eigenvalues/ComplexSchur.h b/Eigen/src/Eigenvalues/ComplexSchur.h index 5deac3247..0fad415a2 100644 --- a/Eigen/src/Eigenvalues/ComplexSchur.h +++ b/Eigen/src/Eigenvalues/ComplexSchur.h @@ -154,6 +154,14 @@ void ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool skipU) m_matT = hess.matrixH(); if(!skipU) m_matU = hess.matrixQ(); + + // Reduce the Hessenberg matrix m_matT to triangular form by QR iteration. + + // The matrix m_matT is divided in three parts. + // Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero. + // Rows il,...,iu is the part we are working on (the active submatrix). + // Rows iu+1,...,end are already brought in triangular form. + int iu = m_matT.cols() - 1; int il; RealScalar d,sd,sf; @@ -164,7 +172,7 @@ void ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool skipU) int iter = 0; while(true) { - //locate the range in which to iterate + // find iu, the bottom row of the active submatrix while(iu > 0) { d = ei_norm1(m_matT.coeff(iu,iu)) + ei_norm1(m_matT.coeff(iu-1,iu-1)); @@ -187,6 +195,7 @@ void ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool skipU) return; } + // find il, the top row of the active submatrix il = iu-1; while(il > 0) { @@ -202,15 +211,16 @@ void ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool skipU) if( il != 0 ) m_matT.coeffRef(il,il-1) = Complex(0); - // compute the shift (the normalization by sf is to avoid under/overflow) + // compute the shift kappa as one of the eigenvalues of the 2x2 + // diagonal block on the bottom of the active submatrix + Matrix<Scalar,2,2> t = m_matT.template block<2,2>(iu-1,iu-1); sf = t.cwiseAbs().sum(); - t /= sf; - - c = t.determinant(); - b = t.diagonal().sum(); + t /= sf; // the normalization by sf is to avoid under/overflow - disc = ei_sqrt(b*b - RealScalar(4)*c); + b = t.coeff(0,0) + t.coeff(1,1); + c = t.coeff(0,0) - t.coeff(1,1); + disc = ei_sqrt(c*c + RealScalar(4)*t.coeff(0,1)*t.coeff(1,0)); r1 = (b+disc)/RealScalar(2); r2 = (b-disc)/RealScalar(2); @@ -224,6 +234,12 @@ void ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool skipU) kappa = sf * r1; else kappa = sf * r2; + + if (iter == 10 || iter == 20) + { + // exceptional shift, taken from http://www.netlib.org/eispack/comqr.f + kappa = ei_abs(ei_real(m_matT.coeff(iu,iu-1))) + ei_abs(ei_real(m_matT.coeff(iu-1,iu-2))); + } // perform the QR step using Givens rotations PlanarRotation<Complex> rot; @@ -246,18 +262,6 @@ void ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool skipU) } } - // FIXME : is it necessary ? - /* - for(int i=0 ; i<n ; i++) - for(int j=0 ; j<n ; j++) - { - if(ei_abs(ei_real(m_matT.coeff(i,j))) < eps) - ei_real_ref(m_matT.coeffRef(i,j)) = 0; - if(ei_imag(ei_abs(m_matT.coeff(i,j))) < eps) - ei_imag_ref(m_matT.coeffRef(i,j)) = 0; - } - */ - m_isInitialized = true; m_matUisUptodate = !skipU; } diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h index 960e9b417..e8e9bea97 100644 --- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h +++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h @@ -276,7 +276,7 @@ inline Matrix<typename NumTraits<typename ei_traits<Derived>::Scalar>::Real, ei_ MatrixBase<Derived>::eigenvalues() const { ei_assert(Flags&SelfAdjoint); - return SelfAdjointEigenSolver<typename Derived::PlainMatrixType>(eval(),false).eigenvalues(); + return SelfAdjointEigenSolver<typename Derived::PlainObject>(eval(),false).eigenvalues(); } template<typename Derived, bool IsSelfAdjoint> @@ -296,7 +296,7 @@ template<typename Derived> struct ei_operatorNorm_selector<Derived, false> static inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real operatorNorm(const MatrixBase<Derived>& m) { - typename Derived::PlainMatrixType m_eval(m); + typename Derived::PlainObject m_eval(m); // FIXME if it is really guaranteed that the eigenvalues are already sorted, // then we don't need to compute a maxCoeff() here, comparing the 1st and last ones is enough. return ei_sqrt( diff --git a/Eigen/src/Geometry/Homogeneous.h b/Eigen/src/Geometry/Homogeneous.h index 76ca66c57..2b8b9cf8e 100644 --- a/Eigen/src/Geometry/Homogeneous.h +++ b/Eigen/src/Geometry/Homogeneous.h @@ -213,9 +213,9 @@ struct ei_traits<ei_homogeneous_left_product_impl<Homogeneous<MatrixType,Vertica typedef Matrix<typename ei_traits<MatrixType>::Scalar, Lhs::RowsAtCompileTime, MatrixType::ColsAtCompileTime, - MatrixType::PlainMatrixType::Options, + MatrixType::PlainObject::Options, Lhs::MaxRowsAtCompileTime, - MatrixType::MaxColsAtCompileTime> ReturnMatrixType; + MatrixType::MaxColsAtCompileTime> ReturnType; }; template<typename MatrixType,typename Lhs> @@ -251,9 +251,9 @@ struct ei_traits<ei_homogeneous_right_product_impl<Homogeneous<MatrixType,Horizo typedef Matrix<typename ei_traits<MatrixType>::Scalar, MatrixType::RowsAtCompileTime, Rhs::ColsAtCompileTime, - MatrixType::PlainMatrixType::Options, + MatrixType::PlainObject::Options, MatrixType::MaxRowsAtCompileTime, - Rhs::MaxColsAtCompileTime> ReturnMatrixType; + Rhs::MaxColsAtCompileTime> ReturnType; }; template<typename MatrixType,typename Rhs> diff --git a/Eigen/src/Geometry/OrthoMethods.h b/Eigen/src/Geometry/OrthoMethods.h index c10b6abf4..265507eb9 100644 --- a/Eigen/src/Geometry/OrthoMethods.h +++ b/Eigen/src/Geometry/OrthoMethods.h @@ -35,7 +35,7 @@ */ template<typename Derived> template<typename OtherDerived> -inline typename MatrixBase<Derived>::PlainMatrixType +inline typename MatrixBase<Derived>::PlainObject MatrixBase<Derived>::cross(const MatrixBase<OtherDerived>& other) const { EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(Derived,3) @@ -79,7 +79,7 @@ struct ei_cross3_impl { */ template<typename Derived> template<typename OtherDerived> -inline typename MatrixBase<Derived>::PlainMatrixType +inline typename MatrixBase<Derived>::PlainObject MatrixBase<Derived>::cross3(const MatrixBase<OtherDerived>& other) const { EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(Derived,4) @@ -210,7 +210,7 @@ struct ei_unitOrthogonal_selector<Derived,2> * \sa cross() */ template<typename Derived> -typename MatrixBase<Derived>::PlainMatrixType +typename MatrixBase<Derived>::PlainObject MatrixBase<Derived>::unitOrthogonal() const { EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) diff --git a/Eigen/src/Geometry/Quaternion.h b/Eigen/src/Geometry/Quaternion.h index 67319d15b..e0dfadd7c 100644 --- a/Eigen/src/Geometry/Quaternion.h +++ b/Eigen/src/Geometry/Quaternion.h @@ -211,6 +211,7 @@ public: template<typename _Scalar> struct ei_traits<Quaternion<_Scalar> > { + typedef Quaternion<_Scalar> PlainObject; typedef _Scalar Scalar; typedef Matrix<_Scalar,4,1> Coefficients; enum{ diff --git a/Eigen/src/Geometry/RotationBase.h b/Eigen/src/Geometry/RotationBase.h index e8bb16f17..9b865889c 100644 --- a/Eigen/src/Geometry/RotationBase.h +++ b/Eigen/src/Geometry/RotationBase.h @@ -74,12 +74,12 @@ class RotationBase */ template<typename OtherDerived> EIGEN_STRONG_INLINE typename ei_rotation_base_generic_product_selector<Derived,OtherDerived,OtherDerived::IsVectorAtCompileTime>::ReturnType - operator*(const AnyMatrixBase<OtherDerived>& e) const + operator*(const EigenBase<OtherDerived>& e) const { return ei_rotation_base_generic_product_selector<Derived,OtherDerived>::run(derived(), e.derived()); } /** \returns the concatenation of a linear transformation \a l with the rotation \a r */ template<typename OtherDerived> friend - inline RotationMatrixType operator*(const AnyMatrixBase<OtherDerived>& l, const Derived& r) + inline RotationMatrixType operator*(const EigenBase<OtherDerived>& l, const Derived& r) { return l.derived() * r.toRotationMatrix(); } /** \returns the concatenation of the rotation \c *this with a transformation \a t */ diff --git a/Eigen/src/Geometry/Transform.h b/Eigen/src/Geometry/Transform.h index 40bc7033a..1240a95bc 100644 --- a/Eigen/src/Geometry/Transform.h +++ b/Eigen/src/Geometry/Transform.h @@ -226,14 +226,14 @@ public: /** Constructs and initializes a transformation from a Dim^2 or a (Dim+1)^2 matrix. */ template<typename OtherDerived> - inline explicit Transform(const AnyMatrixBase<OtherDerived>& other) + inline explicit Transform(const EigenBase<OtherDerived>& other) { ei_transform_construct_from_matrix<OtherDerived,Mode,Dim,HDim>::run(this, other.derived()); } /** Set \c *this from a Dim^2 or (Dim+1)^2 matrix. */ template<typename OtherDerived> - inline Transform& operator=(const AnyMatrixBase<OtherDerived>& other) + inline Transform& operator=(const EigenBase<OtherDerived>& other) { ei_transform_construct_from_matrix<OtherDerived,Mode,Dim,HDim>::run(this, other.derived()); return *this; @@ -310,7 +310,7 @@ public: // note: this function is defined here because some compilers cannot find the respective declaration template<typename OtherDerived> inline const typename ei_transform_right_product_impl<OtherDerived,Mode,_Dim,_Dim+1>::ResultType - operator * (const AnyMatrixBase<OtherDerived> &other) const + operator * (const EigenBase<OtherDerived> &other) const { return ei_transform_right_product_impl<OtherDerived,Mode,Dim,HDim>::run(*this,other.derived()); } /** \returns the product expression of a transformation matrix \a a times a transform \a b @@ -322,11 +322,11 @@ public: */ template<typename OtherDerived> friend inline const typename ei_transform_left_product_impl<OtherDerived,Mode,_Dim,_Dim+1>::ResultType - operator * (const AnyMatrixBase<OtherDerived> &a, const Transform &b) + operator * (const EigenBase<OtherDerived> &a, const Transform &b) { return ei_transform_left_product_impl<OtherDerived,Mode,Dim,HDim>::run(a.derived(),b); } template<typename OtherDerived> - inline Transform& operator*=(const AnyMatrixBase<OtherDerived>& other) { return *this = *this * other; } + inline Transform& operator*=(const EigenBase<OtherDerived>& other) { return *this = *this * other; } /** Contatenates two transformations */ inline const Transform operator * (const Transform& other) const @@ -1021,7 +1021,7 @@ struct ei_transform_construct_from_matrix<Other, AffineCompact,Dim,HDim, HDim,HD }; /********************************************************* -*** Specializations of operator* with a AnyMatrixBase *** +*** Specializations of operator* with a EigenBase *** *********************************************************/ // ei_general_product_return_type is a generalization of ProductReturnType, for all types (including e.g. DiagonalBase...), diff --git a/Eigen/src/Geometry/Translation.h b/Eigen/src/Geometry/Translation.h index 1e3906bde..bde72c738 100644 --- a/Eigen/src/Geometry/Translation.h +++ b/Eigen/src/Geometry/Translation.h @@ -93,7 +93,7 @@ public: /** Concatenates a translation and a linear transformation */ template<typename OtherDerived> - inline AffineTransformType operator* (const AnyMatrixBase<OtherDerived>& linear) const; + inline AffineTransformType operator* (const EigenBase<OtherDerived>& linear) const; /** Concatenates a translation and a rotation */ template<typename Derived> @@ -103,7 +103,7 @@ public: /** \returns the concatenation of a linear transformation \a l with the translation \a t */ // its a nightmare to define a templated friend function outside its declaration template<typename OtherDerived> friend - inline AffineTransformType operator*(const AnyMatrixBase<OtherDerived>& linear, const Translation& t) + inline AffineTransformType operator*(const EigenBase<OtherDerived>& linear, const Translation& t) { AffineTransformType res; res.matrix().setZero(); @@ -182,7 +182,7 @@ Translation<Scalar,Dim>::operator* (const UniformScaling<Scalar>& other) const template<typename Scalar, int Dim> template<typename OtherDerived> inline typename Translation<Scalar,Dim>::AffineTransformType -Translation<Scalar,Dim>::operator* (const AnyMatrixBase<OtherDerived>& linear) const +Translation<Scalar,Dim>::operator* (const EigenBase<OtherDerived>& linear) const { AffineTransformType res; res.matrix().setZero(); diff --git a/Eigen/src/Householder/Householder.h b/Eigen/src/Householder/Householder.h index d86e287fa..9d1383543 100644 --- a/Eigen/src/Householder/Householder.h +++ b/Eigen/src/Householder/Householder.h @@ -99,7 +99,7 @@ void MatrixBase<Derived>::applyHouseholderOnTheLeft( const Scalar& tau, Scalar* workspace) { - Map<Matrix<Scalar, 1, Base::ColsAtCompileTime, PlainMatrixType::Options, 1, Base::MaxColsAtCompileTime> > tmp(workspace,cols()); + Map<Matrix<Scalar, 1, Base::ColsAtCompileTime, PlainObject::Options, 1, Base::MaxColsAtCompileTime> > tmp(workspace,cols()); Block<Derived, EssentialPart::SizeAtCompileTime, Derived::ColsAtCompileTime> bottom(derived(), 1, 0, rows()-1, cols()); tmp.noalias() = essential.adjoint() * bottom; tmp += this->row(0); @@ -114,7 +114,7 @@ void MatrixBase<Derived>::applyHouseholderOnTheRight( const Scalar& tau, Scalar* workspace) { - Map<Matrix<Scalar, Base::RowsAtCompileTime, 1, PlainMatrixType::Options, Base::MaxRowsAtCompileTime, 1> > tmp(workspace,rows()); + Map<Matrix<Scalar, Base::RowsAtCompileTime, 1, PlainObject::Options, Base::MaxRowsAtCompileTime, 1> > tmp(workspace,rows()); Block<Derived, Derived::RowsAtCompileTime, EssentialPart::SizeAtCompileTime> right(derived(), 0, 1, rows(), cols()-1); tmp.noalias() = right * essential.conjugate(); tmp += this->col(0); diff --git a/Eigen/src/Householder/HouseholderSequence.h b/Eigen/src/Householder/HouseholderSequence.h index f8ae4136b..05e883168 100644 --- a/Eigen/src/Householder/HouseholderSequence.h +++ b/Eigen/src/Householder/HouseholderSequence.h @@ -97,7 +97,7 @@ template<typename OtherScalarType, typename MatrixType> struct ei_matrix_type_ti }; template<typename VectorsType, typename CoeffsType, int Side> class HouseholderSequence - : public AnyMatrixBase<HouseholderSequence<VectorsType,CoeffsType,Side> > + : public EigenBase<HouseholderSequence<VectorsType,CoeffsType,Side> > { enum { RowsAtCompileTime = ei_traits<HouseholderSequence>::RowsAtCompileTime, diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h index 72e878223..dea6ec41c 100644 --- a/Eigen/src/LU/FullPivLU.h +++ b/Eigen/src/LU/FullPivLU.h @@ -119,7 +119,7 @@ template<typename _MatrixType> class FullPivLU * diagonal coefficient of U. */ RealScalar maxPivot() const { return m_maxpivot; } - + /** \returns the permutation matrix P * * \sa permutationQ() @@ -251,6 +251,7 @@ template<typename _MatrixType> class FullPivLU { m_usePrescribedThreshold = true; m_prescribedThreshold = threshold; + return *this; } /** Allows to come back to the default behavior, letting Eigen use its default formula for @@ -360,6 +361,8 @@ template<typename _MatrixType> class FullPivLU (*this, MatrixType::Identity(m_lu.rows(), m_lu.cols())); } + MatrixType reconstructedMatrix() const; + inline int rows() const { return m_lu.rows(); } inline int cols() const { return m_lu.cols(); } @@ -403,6 +406,7 @@ FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const MatrixType& matrix) m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case) m_maxpivot = RealScalar(0); + RealScalar cutoff(0); for(int k = 0; k < size; ++k) { @@ -417,8 +421,14 @@ FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const MatrixType& matrix) row_of_biggest_in_corner += k; // correct the values! since they were computed in the corner, col_of_biggest_in_corner += k; // need to add k to them. - // if the pivot (hence the corner) is exactly zero, terminate to avoid generating nan/inf values - if(biggest_in_corner == RealScalar(0)) + // when k==0, biggest_in_corner is the biggest coeff absolute value in the original matrix + if(k == 0) cutoff = biggest_in_corner * NumTraits<Scalar>::epsilon(); + + // if the pivot (hence the corner) is "zero", terminate to avoid generating nan/inf values. + // Notice that using an exact comparison (biggest_in_corner==0) here, as Golub-van Loan do in + // their pseudo-code, results in numerical instability! The cutoff here has been validated + // by running the unit test 'lu' with many repetitions. + if(biggest_in_corner < cutoff) { // before exiting, make sure to initialize the still uninitialized transpositions // in a sane state without destroying what we already have. @@ -479,6 +489,31 @@ typename ei_traits<MatrixType>::Scalar FullPivLU<MatrixType>::determinant() cons return Scalar(m_det_pq) * Scalar(m_lu.diagonal().prod()); } +/** \returns the matrix represented by the decomposition, + * i.e., it returns the product: P^{-1} L U Q^{-1}. + * This function is provided for debug purpose. */ +template<typename MatrixType> +MatrixType FullPivLU<MatrixType>::reconstructedMatrix() const +{ + ei_assert(m_isInitialized && "LU is not initialized."); + const int smalldim = std::min(m_lu.rows(), m_lu.cols()); + // LU + MatrixType res(m_lu.rows(),m_lu.cols()); + // FIXME the .toDenseMatrix() should not be needed... + res = m_lu.corner(TopLeft,m_lu.rows(),smalldim) + .template triangularView<UnitLower>().toDenseMatrix() + * m_lu.corner(TopLeft,smalldim,m_lu.cols()) + .template triangularView<Upper>().toDenseMatrix(); + + // P^{-1}(LU) + res = m_p.inverse() * res; + + // (P^{-1}LU)Q^{-1} + res = res * m_q.inverse(); + + return res; +} + /********* Implementation of kernel() **************************************************/ template<typename _MatrixType> @@ -630,7 +665,7 @@ struct ei_solve_retval<FullPivLU<_MatrixType>, Rhs> return; } - typename Rhs::PlainMatrixType c(rhs().rows(), rhs().cols()); + typename Rhs::PlainObject c(rhs().rows(), rhs().cols()); // Step 1 c = dec().permutationP() * rhs(); @@ -670,10 +705,10 @@ struct ei_solve_retval<FullPivLU<_MatrixType>, Rhs> * \sa class FullPivLU */ template<typename Derived> -inline const FullPivLU<typename MatrixBase<Derived>::PlainMatrixType> +inline const FullPivLU<typename MatrixBase<Derived>::PlainObject> MatrixBase<Derived>::fullPivLu() const { - return FullPivLU<PlainMatrixType>(eval()); + return FullPivLU<PlainObject>(eval()); } #endif // EIGEN_LU_H diff --git a/Eigen/src/LU/Inverse.h b/Eigen/src/LU/Inverse.h index 36392c8d8..e20da70d6 100644 --- a/Eigen/src/LU/Inverse.h +++ b/Eigen/src/LU/Inverse.h @@ -238,7 +238,7 @@ struct ei_compute_inverse_and_det_with_check<MatrixType, ResultType, 4> template<typename MatrixType> struct ei_traits<ei_inverse_impl<MatrixType> > { - typedef typename MatrixType::PlainMatrixType ReturnMatrixType; + typedef typename MatrixType::PlainObject ReturnType; }; template<typename MatrixType> @@ -327,7 +327,7 @@ inline void MatrixBase<Derived>::computeInverseAndDetWithCheck( typedef typename ei_meta_if< RowsAtCompileTime == 2, typename ei_cleantype<typename ei_nested<Derived, 2>::type>::type, - PlainMatrixType + PlainObject >::ret MatrixType; ei_compute_inverse_and_det_with_check<MatrixType, ResultType>::run (derived(), absDeterminantThreshold, inverse, determinant, invertible); diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h index 3925ac1b0..a60596668 100644 --- a/Eigen/src/LU/PartialPivLU.h +++ b/Eigen/src/LU/PartialPivLU.h @@ -165,6 +165,8 @@ template<typename _MatrixType> class PartialPivLU */ typename ei_traits<MatrixType>::Scalar determinant() const; + MatrixType reconstructedMatrix() const; + inline int rows() const { return m_lu.rows(); } inline int cols() const { return m_lu.cols(); } @@ -400,6 +402,23 @@ typename ei_traits<MatrixType>::Scalar PartialPivLU<MatrixType>::determinant() c return Scalar(m_det_p) * m_lu.diagonal().prod(); } +/** \returns the matrix represented by the decomposition, + * i.e., it returns the product: P^{-1} L U. + * This function is provided for debug purpose. */ +template<typename MatrixType> +MatrixType PartialPivLU<MatrixType>::reconstructedMatrix() const +{ + ei_assert(m_isInitialized && "LU is not initialized."); + // LU + MatrixType res = m_lu.template triangularView<UnitLower>().toDenseMatrix() + * m_lu.template triangularView<Upper>(); + + // P^{-1}(LU) + res = m_p.inverse() * res; + + return res; +} + /***** Implementation of solve() *****************************************************/ template<typename _MatrixType, typename Rhs> @@ -442,10 +461,10 @@ struct ei_solve_retval<PartialPivLU<_MatrixType>, Rhs> * \sa class PartialPivLU */ template<typename Derived> -inline const PartialPivLU<typename MatrixBase<Derived>::PlainMatrixType> +inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject> MatrixBase<Derived>::partialPivLu() const { - return PartialPivLU<PlainMatrixType>(eval()); + return PartialPivLU<PlainObject>(eval()); } /** \lu_module @@ -457,10 +476,10 @@ MatrixBase<Derived>::partialPivLu() const * \sa class PartialPivLU */ template<typename Derived> -inline const PartialPivLU<typename MatrixBase<Derived>::PlainMatrixType> +inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject> MatrixBase<Derived>::lu() const { - return PartialPivLU<PlainMatrixType>(eval()); + return PartialPivLU<PlainObject>(eval()); } #endif // EIGEN_PARTIALLU_H diff --git a/Eigen/src/QR/ColPivHouseholderQR.h b/Eigen/src/QR/ColPivHouseholderQR.h index bf28605dd..1219f1918 100644 --- a/Eigen/src/QR/ColPivHouseholderQR.h +++ b/Eigen/src/QR/ColPivHouseholderQR.h @@ -441,7 +441,7 @@ struct ei_solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs> return; } - typename Rhs::PlainMatrixType c(rhs()); + typename Rhs::PlainObject c(rhs()); // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T c.applyOnTheLeft(householderSequence( @@ -458,7 +458,7 @@ struct ei_solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs> .solveInPlace(c.corner(TopLeft, nonzero_pivots, c.cols())); - typename Rhs::PlainMatrixType d(c); + typename Rhs::PlainObject d(c); d.corner(TopLeft, nonzero_pivots, c.cols()) = dec().matrixQR() .corner(TopLeft, nonzero_pivots, nonzero_pivots) @@ -486,10 +486,10 @@ typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHousehol * \sa class ColPivHouseholderQR */ template<typename Derived> -const ColPivHouseholderQR<typename MatrixBase<Derived>::PlainMatrixType> +const ColPivHouseholderQR<typename MatrixBase<Derived>::PlainObject> MatrixBase<Derived>::colPivHouseholderQr() const { - return ColPivHouseholderQR<PlainMatrixType>(eval()); + return ColPivHouseholderQR<PlainObject>(eval()); } diff --git a/Eigen/src/QR/FullPivHouseholderQR.h b/Eigen/src/QR/FullPivHouseholderQR.h index 1ec60aeaf..07be47f47 100644 --- a/Eigen/src/QR/FullPivHouseholderQR.h +++ b/Eigen/src/QR/FullPivHouseholderQR.h @@ -352,7 +352,7 @@ struct ei_solve_retval<FullPivHouseholderQR<_MatrixType>, Rhs> return; } - typename Rhs::PlainMatrixType c(rhs()); + typename Rhs::PlainObject c(rhs()); Matrix<Scalar,1,Rhs::ColsAtCompileTime> temp(rhs().cols()); for (int k = 0; k < dec().rank(); ++k) @@ -413,10 +413,10 @@ typename FullPivHouseholderQR<MatrixType>::MatrixQType FullPivHouseholderQR<Matr * \sa class FullPivHouseholderQR */ template<typename Derived> -const FullPivHouseholderQR<typename MatrixBase<Derived>::PlainMatrixType> +const FullPivHouseholderQR<typename MatrixBase<Derived>::PlainObject> MatrixBase<Derived>::fullPivHouseholderQr() const { - return FullPivHouseholderQR<PlainMatrixType>(eval()); + return FullPivHouseholderQR<PlainObject>(eval()); } #endif // EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H diff --git a/Eigen/src/QR/HouseholderQR.h b/Eigen/src/QR/HouseholderQR.h index b6b07ea63..4709e4b77 100644 --- a/Eigen/src/QR/HouseholderQR.h +++ b/Eigen/src/QR/HouseholderQR.h @@ -221,7 +221,7 @@ struct ei_solve_retval<HouseholderQR<_MatrixType>, Rhs> const int rank = std::min(rows, cols); ei_assert(rhs().rows() == rows); - typename Rhs::PlainMatrixType c(rhs()); + typename Rhs::PlainObject c(rhs()); // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T c.applyOnTheLeft(householderSequence( @@ -246,10 +246,10 @@ struct ei_solve_retval<HouseholderQR<_MatrixType>, Rhs> * \sa class HouseholderQR */ template<typename Derived> -const HouseholderQR<typename MatrixBase<Derived>::PlainMatrixType> +const HouseholderQR<typename MatrixBase<Derived>::PlainObject> MatrixBase<Derived>::householderQr() const { - return HouseholderQR<PlainMatrixType>(eval()); + return HouseholderQR<PlainObject>(eval()); } diff --git a/Eigen/src/SVD/SVD.h b/Eigen/src/SVD/SVD.h index c308ff3ee..fa3b82ce2 100644 --- a/Eigen/src/SVD/SVD.h +++ b/Eigen/src/SVD/SVD.h @@ -555,10 +555,10 @@ void SVD<MatrixType>::computeScalingRotation(ScalingType *scaling, RotationType * \returns the SVD decomposition of \c *this */ template<typename Derived> -inline SVD<typename MatrixBase<Derived>::PlainMatrixType> +inline SVD<typename MatrixBase<Derived>::PlainObject> MatrixBase<Derived>::svd() const { - return SVD<PlainMatrixType>(derived()); + return SVD<PlainObject>(derived()); } #endif // EIGEN_SVD_H diff --git a/Eigen/src/SVD/UpperBidiagonalization.h b/Eigen/src/SVD/UpperBidiagonalization.h index a2cb44733..0a63b4265 100644 --- a/Eigen/src/SVD/UpperBidiagonalization.h +++ b/Eigen/src/SVD/UpperBidiagonalization.h @@ -141,10 +141,10 @@ UpperBidiagonalization<_MatrixType>& UpperBidiagonalization<_MatrixType>::comput * \sa class Bidiagonalization */ template<typename Derived> -const UpperBidiagonalization<typename MatrixBase<Derived>::PlainMatrixType> +const UpperBidiagonalization<typename MatrixBase<Derived>::PlainObject> MatrixBase<Derived>::bidiagonalization() const { - return UpperBidiagonalization<PlainMatrixType>(eval()); + return UpperBidiagonalization<PlainObject>(eval()); } #endif diff --git a/Eigen/src/Sparse/SparseMatrixBase.h b/Eigen/src/Sparse/SparseMatrixBase.h index 205325939..cf1a5d7bf 100644 --- a/Eigen/src/Sparse/SparseMatrixBase.h +++ b/Eigen/src/Sparse/SparseMatrixBase.h @@ -36,7 +36,7 @@ * * */ -template<typename Derived> class SparseMatrixBase : public AnyMatrixBase<Derived> +template<typename Derived> class SparseMatrixBase : public EigenBase<Derived> { public: @@ -109,7 +109,7 @@ template<typename Derived> class SparseMatrixBase : public AnyMatrixBase<Derived Transpose<Derived> >::ret AdjointReturnType; - typedef SparseMatrix<Scalar, Flags&RowMajorBit ? RowMajor : ColMajor> PlainMatrixType; + typedef SparseMatrix<Scalar, Flags&RowMajorBit ? RowMajor : ColMajor> PlainObject; #define EIGEN_CURRENT_STORAGE_BASE_CLASS Eigen::SparseMatrixBase #include "../plugins/CommonCwiseUnaryOps.h" @@ -396,7 +396,7 @@ template<typename Derived> class SparseMatrixBase : public AnyMatrixBase<Derived template<typename OtherDerived> Scalar dot(const SparseMatrixBase<OtherDerived>& other) const; RealScalar squaredNorm() const; RealScalar norm() const; -// const PlainMatrixType normalized() const; +// const PlainObject normalized() const; // void normalize(); Transpose<Derived> transpose() { return derived(); } diff --git a/Eigen/src/Sparse/SparseSelfAdjointView.h b/Eigen/src/Sparse/SparseSelfAdjointView.h index 039e5c725..f3d4fbcbd 100644 --- a/Eigen/src/Sparse/SparseSelfAdjointView.h +++ b/Eigen/src/Sparse/SparseSelfAdjointView.h @@ -93,8 +93,8 @@ template<typename MatrixType, unsigned int UpLo> class SparseSelfAdjointView template<typename DerivedU> SparseSelfAdjointView& rankUpdate(const MatrixBase<DerivedU>& u, Scalar alpha = Scalar(1)); - // const SparseLLT<PlainMatrixType, UpLo> llt() const; - // const SparseLDLT<PlainMatrixType, UpLo> ldlt() const; + // const SparseLLT<PlainObject, UpLo> llt() const; + // const SparseLDLT<PlainObject, UpLo> ldlt() const; protected: diff --git a/Eigen/src/misc/Image.h b/Eigen/src/misc/Image.h index 2f39d6b9d..1d63d8143 100644 --- a/Eigen/src/misc/Image.h +++ b/Eigen/src/misc/Image.h @@ -40,7 +40,7 @@ struct ei_traits<ei_image_retval_base<DecompositionType> > MatrixType::Options, MatrixType::MaxRowsAtCompileTime, // the image matrix will consist of columns from the original matrix, MatrixType::MaxColsAtCompileTime // so it has the same number of rows and at most as many columns. - > ReturnMatrixType; + > ReturnType; }; template<typename _DecompositionType> struct ei_image_retval_base diff --git a/Eigen/src/misc/Kernel.h b/Eigen/src/misc/Kernel.h index 908c408e9..497b42eab 100644 --- a/Eigen/src/misc/Kernel.h +++ b/Eigen/src/misc/Kernel.h @@ -42,7 +42,7 @@ struct ei_traits<ei_kernel_retval_base<DecompositionType> > MatrixType::MaxColsAtCompileTime, // see explanation for 2nd template parameter MatrixType::MaxColsAtCompileTime // the kernel is a subspace of the domain space, // whose dimension is the number of columns of the original matrix - > ReturnMatrixType; + > ReturnType; }; template<typename _DecompositionType> struct ei_kernel_retval_base diff --git a/Eigen/src/misc/Solve.h b/Eigen/src/misc/Solve.h index 4ab0775fc..028716aa2 100644 --- a/Eigen/src/misc/Solve.h +++ b/Eigen/src/misc/Solve.h @@ -35,9 +35,9 @@ struct ei_traits<ei_solve_retval_base<DecompositionType, Rhs> > typedef Matrix<typename Rhs::Scalar, MatrixType::ColsAtCompileTime, Rhs::ColsAtCompileTime, - Rhs::PlainMatrixType::Options, + Rhs::PlainObject::Options, MatrixType::MaxColsAtCompileTime, - Rhs::MaxColsAtCompileTime> ReturnMatrixType; + Rhs::MaxColsAtCompileTime> ReturnType; }; template<typename _DecompositionType, typename Rhs> struct ei_solve_retval_base |