diff options
| -rw-r--r-- | CMakeLists.txt | 2 | ||||
| -rw-r--r-- | Eigen/Core | 7 | ||||
| -rw-r--r-- | Eigen/CoreDeclarations | 4 | ||||
| -rw-r--r-- | Eigen/src/Cholesky/Cholesky.h | 6 | ||||
| -rw-r--r-- | Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h | 8 | ||||
| -rw-r--r-- | Eigen/src/Core/Assign.h | 7 | ||||
| -rw-r--r-- | Eigen/src/Core/CwiseBinaryOp.h | 6 | ||||
| -rw-r--r-- | Eigen/src/Core/Flagged.h | 48 | ||||
| -rwxr-xr-x | Eigen/src/Core/InverseProduct.h | 83 | ||||
| -rw-r--r-- | Eigen/src/Core/MatrixBase.h | 51 | ||||
| -rw-r--r-- | Eigen/src/Core/Part.h | 220 | ||||
| -rw-r--r-- | Eigen/src/Core/ProductWIP.h | 127 | ||||
| -rw-r--r-- | Eigen/src/Core/Transpose.h | 2 | ||||
| -rwxr-xr-x | Eigen/src/Core/Triangular.h | 220 | ||||
| -rw-r--r-- | Eigen/src/Core/TriangularAssign.h | 104 | ||||
| -rw-r--r-- | Eigen/src/Core/util/Constants.h | 28 | ||||
| -rw-r--r-- | Eigen/src/Core/util/ForwardDeclarations.h | 5 | ||||
| -rw-r--r-- | Eigen/src/Core/util/Macros.h | 2 | ||||
| -rw-r--r-- | Eigen/src/Core/util/Meta.h | 6 | ||||
| -rw-r--r-- | Eigen/src/LU/Determinant.h | 4 | ||||
| -rw-r--r-- | Eigen/src/QR/QR.h | 2 | ||||
| -rw-r--r-- | test/triangular.cpp | 20 |
22 files changed, 529 insertions, 433 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt index 95a8a2311..b533e3b66 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -7,7 +7,7 @@ SET(CMAKE_INCLUDE_CURRENT_DIR ON) IF(CMAKE_COMPILER_IS_GNUCXX) IF(CMAKE_SYSTEM_NAME MATCHES Linux) - SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -pedantic -Wnon-virtual-dtor -Wno-long-long -ansi -Wundef -Wcast-align -Wchar-subscripts -Wall -W -Wpointer-arith -Wwrite-strings -Wformat-security -fno-exceptions -fno-check-new -fno-common -fstrict-aliasing") + SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O1 -g -pedantic -Wnon-virtual-dtor -Wno-long-long -ansi -Wundef -Wcast-align -Wchar-subscripts -Wall -W -Wpointer-arith -Wwrite-strings -Wformat-security -fno-exceptions -fno-check-new -fno-common -fstrict-aliasing") IF(TEST_OPENMP) SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fopenmp") MESSAGE("Enabling OpenMP in tests/examples") diff --git a/Eigen/Core b/Eigen/Core index e5310d732..7af09300c 100644 --- a/Eigen/Core +++ b/Eigen/Core @@ -35,11 +35,7 @@ namespace Eigen { #include "src/Core/CwiseBinaryOp.h" #include "src/Core/CwiseUnaryOp.h" #include "src/Core/CwiseNullaryOp.h" -#if (defined EIGEN_WIP_PRODUCT) #include "src/Core/ProductWIP.h" -#else -#include "src/Core/Product.h" -#endif #include "src/Core/Block.h" #include "src/Core/Minor.h" #include "src/Core/Transpose.h" @@ -54,7 +50,8 @@ namespace Eigen { #include "src/Core/Swap.h" #include "src/Core/CommaInitializer.h" #include "src/Core/Triangular.h" -#include "src/Core/TriangularAssign.h" +#include "src/Core/Part.h" +#include "src/Core/InverseProduct.h" } // namespace Eigen diff --git a/Eigen/CoreDeclarations b/Eigen/CoreDeclarations index d141a7c08..40bbb011e 100644 --- a/Eigen/CoreDeclarations +++ b/Eigen/CoreDeclarations @@ -1,10 +1,8 @@ #ifndef EIGEN_CORE_DECLARATIONS_H #define EIGEN_CORE_DECLARATIONS_H -#define EIGEN_WIP_PRODUCT - #ifdef __GNUC__ -#define EIGEN_GNUC_AT_LEAST(x,y) (__GNUC__>=x && __GNUC_MINOR__>=y) || __GNUC__>x +#define EIGEN_GNUC_AT_LEAST(x,y) ((__GNUC__>=x && __GNUC_MINOR__>=y) || __GNUC__>x) #else #define EIGEN_GNUC_AT_LEAST(x,y) 0 #endif diff --git a/Eigen/src/Cholesky/Cholesky.h b/Eigen/src/Cholesky/Cholesky.h index 325c9c6fc..8d9b6e8d8 100644 --- a/Eigen/src/Cholesky/Cholesky.h +++ b/Eigen/src/Cholesky/Cholesky.h @@ -58,9 +58,9 @@ template<typename MatrixType> class Cholesky compute(matrix); } - Triangular<Lower, MatrixType> matrixL(void) const + Extract<MatrixType, Lower> matrixL(void) const { - return m_matrix.lower(); + return m_matrix; } bool isPositiveDefinite(void) const { return m_isPositiveDefinite; } @@ -118,7 +118,7 @@ typename Derived::Eval Cholesky<MatrixType>::solve(MatrixBase<Derived> &b) const int size = m_matrix.rows(); ei_assert(size==b.size()); - return m_matrix.adjoint().upper().inverseProduct(m_matrix.lower().inverseProduct(b)); + return m_matrix.adjoint().template extract<Upper>().inverseProduct(matrixL().inverseProduct(b)); } diff --git a/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h b/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h index 589dd858d..5ebd04d2c 100644 --- a/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h +++ b/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h @@ -58,9 +58,9 @@ template<typename MatrixType> class CholeskyWithoutSquareRoot } /** \returns the lower triangular matrix L */ - Triangular<Lower|UnitDiagBit, MatrixType > matrixL(void) const + Extract<MatrixType, UnitLower> matrixL(void) const { - return m_matrix.lowerWithUnitDiag(); + return m_matrix; } /** \returns the coefficients of the diagonal matrix D */ @@ -131,9 +131,9 @@ typename Derived::Eval CholeskyWithoutSquareRoot<MatrixType>::solve(MatrixBase<D const int size = m_matrix.rows(); ei_assert(size==vecB.size()); - return m_matrix.adjoint().upperWithUnitDiag() + return m_matrix.adjoint().template extract<UnitUpper>() .inverseProduct( - (m_matrix.lowerWithUnitDiag() + (matrixL() .inverseProduct(vecB)) .cwiseQuotient(m_matrix.diagonal()) ); diff --git a/Eigen/src/Core/Assign.h b/Eigen/src/Core/Assign.h index 1eb88d8ff..56f4c956e 100644 --- a/Eigen/src/Core/Assign.h +++ b/Eigen/src/Core/Assign.h @@ -104,8 +104,7 @@ bool Vectorize = (int(Derived::Flags) & int(OtherDerived::Flags) & VectorizableB && ( (int(Derived::Flags) & int(OtherDerived::Flags) & Like1DArrayBit) ||((int(Derived::Flags)&RowMajorBit) ? int(Derived::ColsAtCompileTime)!=Dynamic && (int(Derived::ColsAtCompileTime)%ei_packet_traits<typename Derived::Scalar>::size==0) - : int(Derived::RowsAtCompileTime)!=Dynamic && (int(Derived::RowsAtCompileTime)%ei_packet_traits<typename Derived::Scalar>::size==0)) ), -bool TriangularAssign = Derived::Flags & (NullLowerBit | NullUpperBit)> + : int(Derived::RowsAtCompileTime)!=Dynamic && (int(Derived::RowsAtCompileTime)%ei_packet_traits<typename Derived::Scalar>::size==0)) )> struct ei_assignment_impl; template<typename Derived> @@ -146,7 +145,7 @@ inline Derived& MatrixBase<Derived> } template <typename Derived, typename OtherDerived> -struct ei_assignment_impl<Derived, OtherDerived, false, false> +struct ei_assignment_impl<Derived, OtherDerived, false> { static void execute(Derived & dst, const OtherDerived & src) { @@ -180,7 +179,7 @@ struct ei_assignment_impl<Derived, OtherDerived, false, false> }; template <typename Derived, typename OtherDerived> -struct ei_assignment_impl<Derived, OtherDerived, true, false> +struct ei_assignment_impl<Derived, OtherDerived, true> { static void execute(Derived & dst, const OtherDerived & src) { diff --git a/Eigen/src/Core/CwiseBinaryOp.h b/Eigen/src/Core/CwiseBinaryOp.h index 3142a7885..c514e2169 100644 --- a/Eigen/src/Core/CwiseBinaryOp.h +++ b/Eigen/src/Core/CwiseBinaryOp.h @@ -65,11 +65,11 @@ struct ei_traits<CwiseBinaryOp<BinaryOp, Lhs, Rhs> > ColsAtCompileTime = Lhs::ColsAtCompileTime, MaxRowsAtCompileTime = Lhs::MaxRowsAtCompileTime, MaxColsAtCompileTime = Lhs::MaxColsAtCompileTime, - Flags = ((int(LhsFlags) | int(RhsFlags)) & ( + Flags = (int(LhsFlags) | int(RhsFlags)) & ( HereditaryBits - | int(LhsFlags) & int(RhsFlags) & Like1DArrayBit + | (int(LhsFlags) & int(RhsFlags) & Like1DArrayBit) | (ei_functor_traits<BinaryOp>::IsVectorizable && ((int(LhsFlags) & RowMajorBit)==(int(RhsFlags) & RowMajorBit)) - ? int(LhsFlags) & int(RhsFlags) & VectorizableBit : 0))), + ? int(LhsFlags) & int(RhsFlags) & VectorizableBit : 0)), CoeffReadCost = LhsCoeffReadCost + RhsCoeffReadCost + ei_functor_traits<BinaryOp>::Cost }; }; diff --git a/Eigen/src/Core/Flagged.h b/Eigen/src/Core/Flagged.h index 527d85b12..9167c8a97 100644 --- a/Eigen/src/Core/Flagged.h +++ b/Eigen/src/Core/Flagged.h @@ -60,50 +60,70 @@ template<typename ExpressionType, unsigned int Added, unsigned int Removed> clas EIGEN_GENERIC_PUBLIC_INTERFACE(Flagged) - inline Flagged(const ExpressionType& matrix) : m_expression(matrix) {} + inline Flagged(const ExpressionType& matrix) : m_matrix(matrix) {} /** \internal */ - inline ExpressionType _expression() const { return m_expression; } + inline ExpressionType _expression() const { return m_matrix; } private: - inline int _rows() const { return m_expression.rows(); } - inline int _cols() const { return m_expression.cols(); } + inline int _rows() const { return m_matrix.rows(); } + inline int _cols() const { return m_matrix.cols(); } + inline int _stride() const { return m_matrix.stride(); } inline const Scalar _coeff(int row, int col) const { - return m_expression.coeff(row, col); + return m_matrix.coeff(row, col); } inline Scalar& _coeffRef(int row, int col) { - return m_expression.const_cast_derived().coeffRef(row, col); + return m_matrix.const_cast_derived().coeffRef(row, col); } template<int LoadMode> inline const PacketScalar _packetCoeff(int row, int col) const { - return m_expression.template packetCoeff<LoadMode>(row, col); + return m_matrix.template packetCoeff<LoadMode>(row, col); } template<int LoadMode> inline void _writePacketCoeff(int row, int col, const PacketScalar& x) { - m_expression.const_cast_derived().template writePacketCoeff<LoadMode>(row, col, x); + m_matrix.const_cast_derived().template writePacketCoeff<LoadMode>(row, col, x); } protected: - const ExpressionType m_expression; + const ExpressionType m_matrix; }; -/** \returns an expression of the temporary version of *this. +/** \returns an expression of *this with added flags */ template<typename Derived> -template<unsigned int Added, unsigned int Removed> -inline const Flagged<Derived, Added, Removed> -MatrixBase<Derived>::flagged() const +template<unsigned int Added> +inline const Flagged<Derived, Added, 0> +MatrixBase<Derived>::marked() const { - return Flagged<Derived, Added, Removed>(derived()); + return derived(); +} + +/** \returns an expression of *this with the following flags removed: + * EvalBeforeNestingBit and EvalBeforeAssigningBit. + */ +template<typename Derived> +inline const Flagged<Derived, 0, EvalBeforeNestingBit | EvalBeforeAssigningBit> +MatrixBase<Derived>::lazy() const +{ + return derived(); +} + +/** \returns an expression of *this with the NestByValueBit flag added. + */ +template<typename Derived> +inline const Flagged<Derived, NestByValueBit, 0> +MatrixBase<Derived>::temporary() const +{ + return derived(); } #endif // EIGEN_FLAGGED_H diff --git a/Eigen/src/Core/InverseProduct.h b/Eigen/src/Core/InverseProduct.h new file mode 100755 index 000000000..057590259 --- /dev/null +++ b/Eigen/src/Core/InverseProduct.h @@ -0,0 +1,83 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#ifndef EIGEN_INVERSEPRODUCT_H +#define EIGEN_INVERSEPRODUCT_H + +/** \returns the product of the inverse of *this with \a other. + * + * This function computes the inverse-matrix matrix product inverse(*this) * \a other + * It works as a forward (resp. backward) substitution if *this is an upper (resp. lower) + * triangular matrix. + */ +template<typename Derived> +template<typename OtherDerived> +typename OtherDerived::Eval MatrixBase<Derived>::inverseProduct(const MatrixBase<OtherDerived>& other) const +{ + assert(cols() == other.rows()); + assert(!(Flags & ZeroDiagBit)); + assert(Flags & (UpperTriangularBit|LowerTriangularBit)); + + typename OtherDerived::Eval res(other.rows(), other.cols()); + + for(int c=0 ; c<other.cols() ; ++c) + { + if(Flags & LowerTriangularBit) + { + // forward substitution + if(Flags & UnitDiagBit) + res.coeffRef(0,c) = other.coeff(0,c); + else + res.coeffRef(0,c) = other.coeff(0,c)/coeff(0, 0); + for(int i=1; i<rows(); ++i) + { + Scalar tmp = other.coeff(i,c) - ((this->row(i).start(i)) * res.col(c).start(i)).coeff(0,0); + if (Flags & UnitDiagBit) + res.coeffRef(i,c) = tmp; + else + res.coeffRef(i,c) = tmp/coeff(i,i); + } + } + else + { + // backward substitution + if(Flags & UnitDiagBit) + res.coeffRef(cols()-1,c) = other.coeff(cols()-1,c); + else + res.coeffRef(cols()-1,c) = other.coeff(cols()-1, c)/coeff(rows()-1, cols()-1); + for(int i=rows()-2 ; i>=0 ; --i) + { + Scalar tmp = other.coeff(i,c) + - ((this->row(i).end(cols()-i-1)) * res.col(c).end(cols()-i-1)).coeff(0,0); + if (Flags & UnitDiagBit) + res.coeffRef(i,c) = tmp; + else + res.coeffRef(i,c) = tmp/coeff(i,i); + } + } + } + return res; +} + +#endif // EIGEN_INVERSEPRODUCT_H diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 50c4edfc8..18525a5d1 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -136,6 +136,12 @@ template<typename Derived> class MatrixBase CoeffReadCost = ei_traits<Derived>::CoeffReadCost }; + MatrixBase() + { + assert(!( (Flags&UnitDiagBit && Flags&ZeroDiagBit) + || (Flags&UpperTriangularBit && Flags&LowerTriangularBit) )); + } + /** This is the "real scalar" type; if the \a Scalar type is already real numbers * (e.g. int, float or double) then \a RealScalar is just the same as \a Scalar. If * \a Scalar is \a std::complex<T> then RealScalar is \a T. @@ -280,6 +286,10 @@ template<typename Derived> class MatrixBase template<typename OtherDerived> Derived& operator*=(const MatrixBase<OtherDerived>& other); + + template<typename OtherDerived> + typename OtherDerived::Eval inverseProduct(const MatrixBase<OtherDerived>& other) const; + //@} /** \name Dot product and related notions @@ -296,7 +306,7 @@ template<typename Derived> class MatrixBase const Transpose<Derived> transpose() const; const Transpose< Flagged<CwiseUnaryOp<ei_scalar_conjugate_op<typename ei_traits<Derived>::Scalar>, Derived> - , TemporaryBit, 0> > + , NestByValueBit, 0> > adjoint() const; //@} @@ -347,6 +357,9 @@ template<typename Derived> class MatrixBase DiagonalCoeffs<Derived> diagonal(); const DiagonalCoeffs<Derived> diagonal() const; + + template<unsigned int PartType> Part<Derived, PartType> part(); + template<unsigned int ExtractType> const Extract<Derived, ExtractType> extract() const; //@} /// \name Generating special matrices @@ -406,6 +419,9 @@ template<typename Derived> class MatrixBase bool isIdentity(RealScalar prec = precision<Scalar>()) const; bool isDiagonal(RealScalar prec = precision<Scalar>()) const; + bool isUpper(RealScalar prec = precision<Scalar>()) const; + bool isLower(RealScalar prec = precision<Scalar>()) const; + template<typename OtherDerived> bool isOrtho(const MatrixBase<OtherDerived>& other, RealScalar prec = precision<Scalar>()) const; @@ -433,24 +449,17 @@ template<typename Derived> class MatrixBase template<typename OtherDerived> void swap(const MatrixBase<OtherDerived>& other); - template<unsigned int Added, unsigned int Removed> - const Flagged<Derived, Added, Removed> flagged() const; - - const Flagged<Derived, 0, EvalBeforeNestingBit | EvalBeforeAssigningBit> lazy() const - { - return derived(); - } - const Flagged<Derived, TemporaryBit, 0> temporary() const - { - return derived(); - } + template<unsigned int Added> + const Flagged<Derived, Added, 0> marked() const; + const Flagged<Derived, 0, EvalBeforeNestingBit | EvalBeforeAssigningBit> lazy() const; + const Flagged<Derived, NestByValueBit, 0> temporary() const; /** \returns number of elements to skip to pass from one row (resp. column) to another * for a row-major (resp. column-major) matrix. * Combined with coeffRef() and the compile times flags, it allows a direct access to the data * of the underlying matrix. */ - int stride(void) const { return derived()._stride(); } + inline int stride(void) const { return derived()._stride(); } //@} /// \name Coefficient-wise operations @@ -553,22 +562,6 @@ template<typename Derived> class MatrixBase { return *static_cast<Derived*>(const_cast<MatrixBase*>(this)); } //@} - /// \name Triangular matrices - //@{ - Triangular<Upper, Derived> upper(void); - const Triangular<Upper, Derived> upper(void) const; - const Triangular<Upper|UnitDiagBit, Derived> upperWithUnitDiag(void) const; - const Triangular<Upper|NullDiagBit, Derived> upperWithNullDiag(void) const; - - Triangular<Lower, Derived> lower(void); - const Triangular<Lower, Derived> lower(void) const; - const Triangular<Lower|UnitDiagBit, Derived> lowerWithUnitDiag(void) const; - const Triangular<Lower|NullDiagBit, Derived> lowerWithNullDiag(void) const; - - bool isUpper(RealScalar prec = precision<Scalar>()) const; - bool isLower(RealScalar prec = precision<Scalar>()) const; - //@} - /** \name LU module * * \code #include <Eigen/LU> \endcode diff --git a/Eigen/src/Core/Part.h b/Eigen/src/Core/Part.h new file mode 100644 index 000000000..74c29f203 --- /dev/null +++ b/Eigen/src/Core/Part.h @@ -0,0 +1,220 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2008 Benoit Jacob <jacob@math.jussieu.fr> +// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#ifndef EIGEN_PART_H +#define EIGEN_PART_H + +template<typename MatrixType, unsigned int Mode> +class Part +{ + public: + Part(MatrixType& matrix) : m_matrix(matrix) {} + template<typename Other> void lazyAssign(const Other& other); + template<typename Other> void operator=(const Other& other); + template<typename Other> void operator+=(const Other& other); + template<typename Other> void operator-=(const Other& other); + void operator*=(const typename ei_traits<MatrixType>::Scalar& other); + void operator/=(const typename ei_traits<MatrixType>::Scalar& other); + void setConstant(const typename ei_traits<MatrixType>::Scalar& value); + void setZero(); + void setOnes(); + void setRandom(); + void setIdentity(); + + private: + MatrixType& m_matrix; +}; + +template<typename MatrixType, unsigned int Mode> +template<typename Other> +inline void Part<MatrixType, Mode>::operator=(const Other& other) +{ + if(Other::Flags & EvalBeforeAssigningBit) + lazyAssign(other.eval()); + else + lazyAssign(other.derived()); +} + +template<typename Derived1, typename Derived2, unsigned int Mode, int UnrollCount> +struct ei_part_assignment_unroller +{ + enum { + col = (UnrollCount-1) / Derived1::RowsAtCompileTime, + row = (UnrollCount-1) % Derived1::RowsAtCompileTime + }; + + inline static void run(Derived1 &dst, const Derived2 &src) + { + ei_part_assignment_unroller<Derived1, Derived2, Mode, UnrollCount-1>::run(dst, src); + + if(Mode == SelfAdjoint) + { + if(row == col) + dst.coeffRef(row, col) = ei_real(src.coeff(row, col)); + else if(row < col) + dst.coeffRef(col, row) = ei_conj(dst.coeffRef(row, col) = src.coeff(row, col)); + } + else + { + if((Mode == Upper && row <= col) + || (Mode == Lower && row >= col) + || (Mode == StrictlyUpper && row < col) + || (Mode == StrictlyLower && row > col)) + dst.coeffRef(row, col) = src.coeff(row, col); + } + } +}; + +template<typename Derived1, typename Derived2, unsigned int Mode> +struct ei_part_assignment_unroller<Derived1, Derived2, Mode, 1> +{ + inline static void run(Derived1 &dst, const Derived2 &src) + { + if(!(Mode & ZeroDiagBit)) + dst.coeffRef(0, 0) = src.coeff(0, 0); + } +}; + +// prevent buggy user code from causing an infinite recursion +template<typename Derived1, typename Derived2, unsigned int Mode> +struct ei_part_assignment_unroller<Derived1, Derived2, Mode, 0> +{ + inline static void run(Derived1 &, const Derived2 &) {} +}; + +template<typename Derived1, typename Derived2, unsigned int Mode> +struct ei_part_assignment_unroller<Derived1, Derived2, Mode, Dynamic> +{ + inline static void run(Derived1 &, const Derived2 &) {} +}; + + +template<typename MatrixType, unsigned int Mode> +template<typename Other> +void Part<MatrixType, Mode>::lazyAssign(const Other& other) +{ + const bool unroll = MatrixType::SizeAtCompileTime * Other::CoeffReadCost / 2 <= EIGEN_UNROLLING_LIMIT; + ei_assert(m_matrix.rows() == other.rows() && m_matrix.cols() == other.cols()); + if(unroll) + { + ei_part_assignment_unroller + <MatrixType, Other, Mode, + unroll ? int(MatrixType::SizeAtCompileTime) : Dynamic + >::run(m_matrix, other.derived()); + } + else + { + switch(Mode) + { + case Upper: + for(int j = 0; j < m_matrix.cols(); j++) + for(int i = 0; i <= j; i++) + m_matrix.coeffRef(i, j) = other.coeff(i, j); + break; + case Lower: + for(int j = 0; j < m_matrix.cols(); j++) + for(int i = j; i < m_matrix.rows(); i++) + m_matrix.coeffRef(i, j) = other.coeff(i, j); + break; + case StrictlyUpper: + for(int j = 0; j < m_matrix.cols(); j++) + for(int i = 0; i < j; i++) + m_matrix.coeffRef(i, j) = other.coeff(i, j); + break; + case StrictlyLower: + for(int j = 0; j < m_matrix.cols(); j++) + for(int i = j+1; i < m_matrix.rows(); i++) + m_matrix.coeffRef(i, j) = other.coeff(i, j); + break; + case SelfAdjoint: + for(int j = 0; j < m_matrix.cols(); j++) + { + for(int i = 0; i < j; i++) + m_matrix.coeffRef(j, i) = ei_conj(m_matrix.coeffRef(i, j) = other.coeff(i, j)); + m_matrix.coeffRef(j, j) = ei_real(other.coeff(j, j)); + } + break; + } + } +} + +template<typename MatrixType, unsigned int Mode> +template<typename Other> inline void Part<MatrixType, Mode>::operator+=(const Other& other) +{ + *this = m_matrix + other; +} + +template<typename MatrixType, unsigned int Mode> +template<typename Other> inline void Part<MatrixType, Mode>::operator-=(const Other& other) +{ + *this = m_matrix - other; +} + +template<typename MatrixType, unsigned int Mode> +inline void Part<MatrixType, Mode>::operator*= +(const typename ei_traits<MatrixType>::Scalar& other) +{ + *this = m_matrix * other; +} + +template<typename MatrixType, unsigned int Mode> +inline void Part<MatrixType, Mode>::operator/= +(const typename ei_traits<MatrixType>::Scalar& other) +{ + *this = m_matrix / other; +} + +template<typename MatrixType, unsigned int Mode> +inline void Part<MatrixType, Mode>::setConstant(const typename ei_traits<MatrixType>::Scalar& value) +{ + *this = MatrixType::constant(m_matrix.rows(), m_matrix.cols(), value); +} + +template<typename MatrixType, unsigned int Mode> +inline void Part<MatrixType, Mode>::setZero() +{ + setConstant((typename ei_traits<MatrixType>::Scalar)(0)); +} + +template<typename MatrixType, unsigned int Mode> +inline void Part<MatrixType, Mode>::setOnes() +{ + setConstant((typename ei_traits<MatrixType>::Scalar)(1)); +} + +template<typename MatrixType, unsigned int Mode> +inline void Part<MatrixType, Mode>::setRandom() +{ + *this = MatrixType::random(m_matrix.rows(), m_matrix.cols()); +} + +template<typename Derived> +template<unsigned int Mode> +inline Part<Derived, Mode> MatrixBase<Derived>::part() +{ + return Part<Derived, Mode>(derived()); +} + +#endif // EIGEN_PART_H diff --git a/Eigen/src/Core/ProductWIP.h b/Eigen/src/Core/ProductWIP.h index 6815b28ab..9e7dc91f1 100644 --- a/Eigen/src/Core/ProductWIP.h +++ b/Eigen/src/Core/ProductWIP.h @@ -167,7 +167,7 @@ template<typename T> class ei_product_eval_to_column_major template<typename T, int n=1> struct ei_product_nested_rhs { typedef typename ei_meta_if< - (ei_is_temporary<T>::ret && !(ei_traits<T>::Flags & RowMajorBit)), + (ei_traits<T>::Flags & NestByValueBit) && !(ei_traits<T>::Flags & RowMajorBit), T, typename ei_meta_if< ((ei_traits<T>::Flags & EvalBeforeNestingBit) @@ -183,7 +183,7 @@ template<typename T, int n=1> struct ei_product_nested_rhs template<typename T, int n=1> struct ei_product_nested_lhs { typedef typename ei_meta_if< - ei_is_temporary<T>::ret, + ei_traits<T>::Flags & NestByValueBit, T, typename ei_meta_if< int(ei_traits<T>::Flags) & EvalBeforeNestingBit @@ -202,7 +202,7 @@ struct ei_traits<Product<Lhs, Rhs, EvalMode> > // the cache friendly product evals lhs once only // FIXME what to do if we chose to dynamically call the normal product from the cache friendly one for small matrices ? typedef typename ei_meta_if<EvalMode==CacheFriendlyProduct, - typename ei_product_nested_lhs<Lhs,0>::type, + typename ei_product_nested_lhs<Lhs,1>::type, typename ei_nested<Lhs,Rhs::ColsAtCompileTime>::type>::ret LhsNested; // NOTE that rhs must be ColumnMajor, so we might need a special nested type calculation @@ -231,9 +231,7 @@ struct ei_traits<Product<Lhs, Rhs, EvalMode> > (_RowMajor ? 0 : RowMajorBit) | ((RowsAtCompileTime == Dynamic || ColsAtCompileTime == Dynamic) ? 0 : LargeBit)), Flags = ((unsigned int)(LhsFlags | RhsFlags) & _LostBits) - #ifndef EIGEN_WIP_PRODUCT_DIRTY - | EvalBeforeAssigningBit //FIXME - #endif + | EvalBeforeAssigningBit | EvalBeforeNestingBit | (_Vectorizable ? VectorizableBit : 0), CoeffReadCost @@ -335,6 +333,9 @@ template<typename Lhs, typename Rhs, int EvalMode> class Product : ei_no_assignm return res; } + template<typename Lhs_, typename Rhs_, int EvalMode_, typename DestDerived_, bool DirectAccess_> + friend struct ei_cache_friendly_selector; + protected: const LhsNested m_lhs; const RhsNested m_rhs; @@ -342,9 +343,6 @@ template<typename Lhs, typename Rhs, int EvalMode> class Product : ei_no_assignm /** \returns the matrix product of \c *this and \a other. * - * \note This function causes an immediate evaluation. If you want to perform a matrix product - * without immediate evaluation, call .lazy() on one of the matrices before taking the product. - * * \sa lazy(), operator*=(const MatrixBase&) */ template<typename Derived> @@ -374,7 +372,7 @@ inline Derived& MatrixBase<Derived>::operator+=(const Flagged<Product<Lhs,Rhs,CacheFriendlyProduct>, 0, EvalBeforeNestingBit | EvalBeforeAssigningBit>& other) { other._expression()._cacheFriendlyEvalAndAdd(const_cast_derived()); - return derived(); + return derived(); } template<typename Derived> @@ -385,54 +383,89 @@ inline Derived& MatrixBase<Derived>::lazyAssign(const Product<Lhs,Rhs,CacheFrien return derived(); } -template<typename Lhs, typename Rhs, int EvalMode> -template<typename DestDerived> -inline void Product<Lhs,Rhs,EvalMode>::_cacheFriendlyEval(DestDerived& res) const +template<typename Lhs, typename Rhs, int EvalMode, typename DestDerived, bool DirectAccess> +struct ei_cache_friendly_selector { - if ( _rows()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD - && _cols()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD - && m_lhs.cols()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD - ) + typedef Product<Lhs,Rhs,EvalMode> Prod; + typedef typename Prod::_LhsNested _LhsNested; + typedef typename Prod::_RhsNested _RhsNested; + typedef typename Prod::Scalar Scalar; + static inline void eval(const Prod& product, DestDerived& res) { - #ifndef EIGEN_WIP_PRODUCT_DIRTY - res.setZero(); - #endif - - ei_cache_friendly_product<Scalar>( - _rows(), _cols(), m_lhs.cols(), - _LhsNested::Flags&RowMajorBit, &(m_lhs.const_cast_derived().coeffRef(0,0)), m_lhs.stride(), - _RhsNested::Flags&RowMajorBit, &(m_rhs.const_cast_derived().coeffRef(0,0)), m_rhs.stride(), - Flags&RowMajorBit, &(res.coeffRef(0,0)), res.stride() - ); + if ( product._rows()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD + && product._cols()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD + && product.m_lhs.cols()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD + ) + { + res.setZero(); + ei_cache_friendly_product<Scalar>( + product._rows(), product._cols(), product.m_lhs.cols(), + _LhsNested::Flags&RowMajorBit, &(product.m_lhs.const_cast_derived().coeffRef(0,0)), product.m_lhs.stride(), + _RhsNested::Flags&RowMajorBit, &(product.m_rhs.const_cast_derived().coeffRef(0,0)), product.m_rhs.stride(), + Prod::Flags&RowMajorBit, &(res.coeffRef(0,0)), res.stride() + ); + } + else + { + res = Product<_LhsNested,_RhsNested,NormalProduct>(product.m_lhs, product.m_rhs).lazy(); + } } - else + + static inline void eval_and_add(const Prod& product, DestDerived& res) { - // lazy product - res = Product<_LhsNested,_RhsNested,NormalProduct>(m_lhs, m_rhs).lazy(); + if ( product._rows()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD + && product._cols()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD + && product.m_lhs.cols()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD + ) + { + ei_cache_friendly_product<Scalar>( + product._rows(), product._cols(), product.m_lhs.cols(), + _LhsNested::Flags&RowMajorBit, &(product.m_lhs.const_cast_derived().coeffRef(0,0)), product.m_lhs.stride(), + _RhsNested::Flags&RowMajorBit, &(product.m_rhs.const_cast_derived().coeffRef(0,0)), product.m_rhs.stride(), + Prod::Flags&RowMajorBit, &(res.coeffRef(0,0)), res.stride() + ); + } + else + { + res += Product<_LhsNested,_RhsNested,NormalProduct>(product.m_lhs, product.m_rhs).lazy(); + } } -} +}; -template<typename Lhs, typename Rhs, int EvalMode> -template<typename DestDerived> -inline void Product<Lhs,Rhs,EvalMode>::_cacheFriendlyEvalAndAdd(DestDerived& res) const +template<typename Lhs, typename Rhs, int EvalMode, typename DestDerived> +struct ei_cache_friendly_selector<Lhs,Rhs,EvalMode,DestDerived,false> { - if ( _rows()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD - && _cols()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD - && m_lhs.cols()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD - ) + typedef Product<Lhs,Rhs,EvalMode> Prod; + typedef typename Prod::_LhsNested _LhsNested; + typedef typename Prod::_RhsNested _RhsNested; + typedef typename Prod::Scalar Scalar; + static inline void eval(const Prod& product, DestDerived& res) { - ei_cache_friendly_product<Scalar>( - _rows(), _cols(), m_lhs.cols(), - _LhsNested::Flags&RowMajorBit, &(m_lhs.const_cast_derived().coeffRef(0,0)), m_lhs.stride(), - _RhsNested::Flags&RowMajorBit, &(m_rhs.const_cast_derived().coeffRef(0,0)), m_rhs.stride(), - Flags&RowMajorBit, &(res.coeffRef(0,0)), res.stride() - ); + res = Product<_LhsNested,_RhsNested,NormalProduct>(product.m_lhs, product.m_rhs).lazy(); } - else + + static inline void eval_and_add(const Prod& product, DestDerived& res) { - // lazy product - res += Product<_LhsNested,_RhsNested,NormalProduct>(m_lhs, m_rhs).lazy(); + res += Product<_LhsNested,_RhsNested,NormalProduct>(product.m_lhs, product.m_rhs).lazy(); } +}; + +template<typename Lhs, typename Rhs, int EvalMode> +template<typename DestDerived> +inline void Product<Lhs,Rhs,EvalMode>::_cacheFriendlyEval(DestDerived& res) const +{ + ei_cache_friendly_selector<Lhs,Rhs,EvalMode,DestDerived, + _LhsNested::Flags&_RhsNested::Flags&DirectAccessBit> + ::eval(*this, res); +} + +template<typename Lhs, typename Rhs, int EvalMode> +template<typename DestDerived> +inline void Product<Lhs,Rhs,EvalMode>::_cacheFriendlyEvalAndAdd(DestDerived& res) const +{ + ei_cache_friendly_selector<Lhs,Rhs,EvalMode,DestDerived, + _LhsNested::Flags&_RhsNested::Flags&DirectAccessBit> + ::eval_and_add(*this, res); } #endif // EIGEN_PRODUCT_H diff --git a/Eigen/src/Core/Transpose.h b/Eigen/src/Core/Transpose.h index 9c19f80f3..fd28f4bab 100644 --- a/Eigen/src/Core/Transpose.h +++ b/Eigen/src/Core/Transpose.h @@ -127,7 +127,7 @@ MatrixBase<Derived>::transpose() const template<typename Derived> inline const Transpose< Flagged<CwiseUnaryOp<ei_scalar_conjugate_op<typename ei_traits<Derived>::Scalar>, Derived > - , TemporaryBit, 0> > + , NestByValueBit, 0> > MatrixBase<Derived>::adjoint() const { return conjugate().temporary().transpose(); diff --git a/Eigen/src/Core/Triangular.h b/Eigen/src/Core/Triangular.h index 12f34cd99..e9323733a 100755 --- a/Eigen/src/Core/Triangular.h +++ b/Eigen/src/Core/Triangular.h @@ -2,6 +2,7 @@ // for linear algebra. Eigen itself is part of the KDE project. // // Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> +// Copyright (C) 2008 Benoit Jacob <jacob@math.jussieu.fr> // // Eigen is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public @@ -22,42 +23,27 @@ // License and a copy of the GNU General Public License along with // Eigen. If not, see <http://www.gnu.org/licenses/>. -#ifndef EIGEN_TRIANGULAR_H -#define EIGEN_TRIANGULAR_H +#ifndef EIGEN_EXTRACT_H +#define EIGEN_EXTRACT_H -/** \class Triangular +/** \class Extract * - * \brief Expression of a triangular matrix from a square matrix + * \brief Expression of a triangular matrix extracted from a given matrix * - * \param Mode or-ed bit field indicating the triangular part (Upper or Lower) we are taking, - * and the property of the diagonal if any (UnitDiagBit or NullDiagBit). * \param MatrixType the type of the object in which we are taking the triangular part + * \param Mode the kind of triangular matrix expression to construct. Can be Upper, StrictlyUpper, + * UnitUpper, Lower, StrictlyLower, UnitLower. This is in fact a bit field; it must have either + * UpperTriangularBit or LowerTriangularBit, and additionnaly it may have either ZeroDiagBit or + * UnitDiagBit. * * This class represents an expression of the upper or lower triangular part of - * a squared matrix. It is the return type of MatrixBase::upper(), MatrixBase::lower(), - * MatrixBase::upperWithUnitDiagBit(), etc., and used to optimize operations involving - * triangular matrices. Most of the time this is the only way it is used. + * a square matrix, possibly with a further assumption on the diagonal. It is the return type + * of MatrixBase::extract() and most of the time this is the only way it is used. * - * Examples of some key features: - * \code - * m1 = (<any expression>).upper(); - * \endcode - * In this example, the strictly lower part of the expression is not evaluated, - * m1 might be resized and the strict lower part of m1 == 0. - * - * \code - * m1.upper() = <any expression>; - * \endcode - * This example diverge from the previous one in the sense that the strictly - * lower part of m1 is left unchanged, and optimal loops are employed. Note that - * m1 might also be resized. - * - * Of course, in both examples \c <any \c expression> has to be a square matrix. - * - * \sa MatrixBase::upper(), MatrixBase::lower(), class TriangularProduct + * \sa MatrixBase::extract() */ -template<int Mode, typename MatrixType> -struct ei_traits<Triangular<Mode, MatrixType> > +template<typename MatrixType, unsigned int Mode> +struct ei_traits<Extract<MatrixType, Mode> > { typedef typename MatrixType::Scalar Scalar; typedef typename ei_nested<MatrixType>::type MatrixTypeNested; @@ -72,105 +58,30 @@ struct ei_traits<Triangular<Mode, MatrixType> > }; }; -template<int Mode, typename MatrixType> class Triangular - : public MatrixBase<Triangular<Mode,MatrixType> > +template<typename MatrixType, unsigned int Mode> class Extract + : public MatrixBase<Extract<MatrixType, Mode> > { public: - EIGEN_GENERIC_PUBLIC_INTERFACE(Triangular) - - inline Triangular(const MatrixType& matrix) - : m_matrix(matrix) - { - assert(!( (Flags&UnitDiagBit) && (Flags&NullDiagBit))); - assert(matrix.rows()==matrix.cols()); - } + EIGEN_GENERIC_PUBLIC_INTERFACE(Extract) - EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Triangular) + inline Extract(const MatrixType& matrix) : m_matrix(matrix) {} - /** Overloaded to keep a Triangular expression */ - inline Triangular<(Upper | Lower) ^ Mode, Flagged<Transpose<MatrixType>,TemporaryBit,0> > transpose() - { - return m_matrix.transpose().temporary(); - } - - /** Overloaded to keep a Triangular expression */ - inline const Triangular<(Upper | Lower) ^ Mode, Flagged<Transpose<MatrixType>,TemporaryBit,0> > transpose() const - { - return m_matrix.transpose().temporary(); - } - - /** \returns the product of the inverse of *this with \a other. - * - * This function computes the inverse-matrix matrix product inverse(*this) * \a other - * It works as a forward (resp. backward) substitution if *this is an upper (resp. lower) - * triangular matrix. - */ - template<typename OtherDerived> - typename OtherDerived::Eval inverseProduct(const MatrixBase<OtherDerived>& other) const - { - assert(_cols() == other.rows()); - assert(!(Flags & NullDiagBit)); - - typename OtherDerived::Eval res(other.rows(), other.cols()); - - for (int c=0 ; c<other.cols() ; ++c) - { - if (Flags & Lower) - { - // forward substitution - if (Flags & UnitDiagBit) - res(0,c) = other(0,c); - else - res(0,c) = other(0,c)/_coeff(0, 0); - for (int i=1 ; i<_rows() ; ++i) - { - Scalar tmp = other(i,c) - ((this->row(i).start(i)) * res.col(c).start(i))(0,0); - if (Flags & UnitDiagBit) - res(i,c) = tmp; - else - res(i,c) = tmp/_coeff(i,i); - } - } - else - { - // backward substitution - if (Flags & UnitDiagBit) - res(_cols()-1,c) = other(_cols()-1,c); - else - res(_cols()-1,c) = other(_cols()-1, c)/_coeff(_rows()-1, _cols()-1); - for (int i=_rows()-2 ; i>=0 ; --i) - { - Scalar tmp = other(i,c) - ((this->row(i).end(_cols()-i-1)) * res.col(c).end(_cols()-i-1))(0,0); - if (Flags & UnitDiagBit) - res(i,c) = tmp; - else - res(i,c) = tmp/_coeff(i,i); - } - } - } - return res; - } + EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Extract) private: inline int _rows() const { return m_matrix.rows(); } inline int _cols() const { return m_matrix.cols(); } - inline Scalar& _coeffRef(int row, int col) - { - ei_assert( ((! (Flags & Lower)) && row<=col) || (Flags & Lower && col<=row)); - return m_matrix.const_cast_derived().coeffRef(row, col); - } - inline Scalar _coeff(int row, int col) const { - if ((Flags & Lower) ? col>row : row>col) - return 0; - if (Flags & UnitDiagBit) - return col==row ? 1 : m_matrix.coeff(row, col); - else if (Flags & NullDiagBit) - return col==row ? 0 : m_matrix.coeff(row, col); + if(Flags & LowerTriangularBit ? col>row : row>col) + return (Scalar)0; + if(Flags & UnitDiagBit) + return col==row ? (Scalar)1 : m_matrix.coeff(row, col); + else if(Flags & ZeroDiagBit) + return col==row ? (Scalar)0 : m_matrix.coeff(row, col); else return m_matrix.coeff(row, col); } @@ -180,86 +91,21 @@ template<int Mode, typename MatrixType> class Triangular const typename MatrixType::Nested m_matrix; }; -/** \returns an expression of a upper triangular matrix - * - * \sa isUpper(), upperWithNullDiagBit(), upperWithNullDiagBit(), lower() - */ -template<typename Derived> -inline Triangular<Upper, Derived> MatrixBase<Derived>::upper(void) -{ - return Triangular<Upper,Derived>(derived()); -} - -/** This is the const version of upper(). */ -template<typename Derived> -inline const Triangular<Upper, Derived> MatrixBase<Derived>::upper(void) const -{ - return Triangular<Upper,Derived>(derived()); -} - -/** \returns an expression of a lower triangular matrix - * - * \sa isLower(), lowerWithUnitDiag(), lowerWithNullDiag(), upper() - */ -template<typename Derived> -inline Triangular<Lower, Derived> MatrixBase<Derived>::lower(void) -{ - return Triangular<Lower,Derived>(derived()); -} - -/** This is the const version of lower().*/ -template<typename Derived> -inline const Triangular<Lower, Derived> MatrixBase<Derived>::lower(void) const -{ - return Triangular<Lower,Derived>(derived()); -} - -/** \returns an expression of a upper triangular matrix with a unit diagonal - * - * \sa upper(), lowerWithUnitDiagBit() - */ -template<typename Derived> -inline const Triangular<Upper|UnitDiagBit, Derived> MatrixBase<Derived>::upperWithUnitDiag(void) const -{ - return Triangular<Upper|UnitDiagBit, Derived>(derived()); -} - -/** \returns an expression of a strictly upper triangular matrix (diagonal==zero) - * FIXME could also be called strictlyUpper() or upperStrict() - * - * \sa upper(), lowerWithNullDiag() - */ -template<typename Derived> -inline const Triangular<Upper|NullDiagBit, Derived> MatrixBase<Derived>::upperWithNullDiag(void) const -{ - return Triangular<Upper|NullDiagBit, Derived>(derived()); -} - -/** \returns an expression of a lower triangular matrix with a unit diagonal - * - * \sa lower(), upperWithUnitDiag() - */ -template<typename Derived> -inline const Triangular<Lower|UnitDiagBit, Derived> MatrixBase<Derived>::lowerWithUnitDiag(void) const -{ - return Triangular<Lower|UnitDiagBit, Derived>(derived()); -} - -/** \returns an expression of a strictly lower triangular matrix (diagonal==zero) - * FIXME could also be called strictlyLower() or lowerStrict() +/** \returns an expression of a triangular matrix extracted from the current matrix * - * \sa lower(), upperWithNullDiag() + * \sa part(), marked() */ template<typename Derived> -inline const Triangular<Lower|NullDiagBit, Derived> MatrixBase<Derived>::lowerWithNullDiag(void) const +template<unsigned int Mode> +const Extract<Derived, Mode> MatrixBase<Derived>::extract() const { - return Triangular<Lower|NullDiagBit, Derived>(derived()); + return derived(); } /** \returns true if *this is approximately equal to an upper triangular matrix, * within the precision given by \a prec. * - * \sa isLower(), upper() + * \sa isLower(), extract(), part(), marked() */ template<typename Derived> bool MatrixBase<Derived>::isUpper(RealScalar prec) const @@ -281,7 +127,7 @@ bool MatrixBase<Derived>::isUpper(RealScalar prec) const /** \returns true if *this is approximately equal to a lower triangular matrix, * within the precision given by \a prec. * - * \sa isUpper(), upper() + * \sa isUpper(), extract(), part(), marked() */ template<typename Derived> bool MatrixBase<Derived>::isLower(RealScalar prec) const @@ -300,4 +146,4 @@ bool MatrixBase<Derived>::isLower(RealScalar prec) const return true; } -#endif // EIGEN_TRIANGULAR_H +#endif // EIGEN_EXTRACT_H diff --git a/Eigen/src/Core/TriangularAssign.h b/Eigen/src/Core/TriangularAssign.h deleted file mode 100644 index 7b1be135e..000000000 --- a/Eigen/src/Core/TriangularAssign.h +++ /dev/null @@ -1,104 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. Eigen itself is part of the KDE project. -// -// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> -// -// Eigen is free software; you can redistribute it and/or -// modify it under the terms of the GNU Lesser General Public -// License as published by the Free Software Foundation; either -// version 3 of the License, or (at your option) any later version. -// -// Alternatively, you can redistribute it and/or -// modify it under the terms of the GNU General Public License as -// published by the Free Software Foundation; either version 2 of -// the License, or (at your option) any later version. -// -// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY -// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the -// GNU General Public License for more details. -// -// You should have received a copy of the GNU Lesser General Public -// License and a copy of the GNU General Public License along with -// Eigen. If not, see <http://www.gnu.org/licenses/>. - -#ifndef EIGEN_TRIANGULAR_ASSIGN_H -#define EIGEN_TRIANGULAR_ASSIGN_H - -template<typename Derived1, typename Derived2, int UnrollCount, int Mode> -struct ei_triangular_assign_unroller -{ - enum { - col = (UnrollCount-1) / Derived1::RowsAtCompileTime, - row = (UnrollCount-1) % Derived1::RowsAtCompileTime - }; - - inline static void run(Derived1 &dst, const Derived2 &src) - { - ei_triangular_assign_unroller<Derived1, Derived2, - (Mode & Lower) ? - ((row==col) ? UnrollCount-1-row : UnrollCount-1) - : ((row==0) ? UnrollCount-1-Derived1::ColsAtCompileTime+col : UnrollCount-1), - Mode>::run(dst, src); - dst.coeffRef(row, col) = src.coeff(row, col); - } -}; - -template<typename Derived1, typename Derived2, int Mode> -struct ei_triangular_assign_unroller<Derived1, Derived2, 1, Mode> -{ - inline static void run(Derived1 &dst, const Derived2 &src) - { - dst.coeffRef(0, 0) = src.coeff(0, 0); - } -}; - -// prevent buggy user code from causing an infinite recursion -template<typename Derived1, typename Derived2, int Mode> -struct ei_triangular_assign_unroller<Derived1, Derived2, 0, Mode> -{ - inline static void run(Derived1 &, const Derived2 &) {} -}; - -template<typename Derived1, typename Derived2, int Mode> -struct ei_triangular_assign_unroller<Derived1, Derived2, Dynamic, Mode> -{ - inline static void run(Derived1 &, const Derived2 &) {} -}; - - -template <typename Derived, typename OtherDerived, bool DummyVectorize> -struct ei_assignment_impl<Derived, OtherDerived, DummyVectorize, true> -{ - static void execute(Derived & dst, const OtherDerived & src) - { - assert(src.rows()==src.cols()); - assert(dst.rows() == src.rows() && dst.cols() == src.cols()); - - const bool unroll = Derived::SizeAtCompileTime * OtherDerived::CoeffReadCost <= EIGEN_UNROLLING_LIMIT; - - if(unroll) - { - ei_triangular_assign_unroller - <Derived, OtherDerived, unroll ? int(Derived::SizeAtCompileTime) : Dynamic, Derived::Flags>::run - (dst.derived(), src.derived()); - } - else - { - if (Derived::Flags & Lower) - { - for(int j = 0; j < dst.cols(); j++) - for(int i = j; i < dst.rows(); i++) - dst.coeffRef(i, j) = src.coeff(i, j); - } - else - { - for(int j = 0; j < dst.cols(); j++) - for(int i = 0; i <= j; i++) - dst.coeffRef(i, j) = src.coeff(i, j); - } - } - } -}; - -#endif // EIGEN_TRIANGULAR_ASSIGN_H diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h index 12efa9530..e599d8a3d 100644 --- a/Eigen/src/Core/util/Constants.h +++ b/Eigen/src/Core/util/Constants.h @@ -39,15 +39,13 @@ const unsigned int VectorizableBit = 0x10; ///< means the expression might be v const unsigned int VectorizableBit = 0x0; #endif const unsigned int Like1DArrayBit = 0x20; ///< means the expression can be seen as 1D vector (used for explicit vectorization) -const unsigned int NullDiagBit = 0x40; ///< means all diagonal coefficients are equal to 0 +const unsigned int ZeroDiagBit = 0x40; ///< means all diagonal coefficients are equal to 0 const unsigned int UnitDiagBit = 0x80; ///< means all diagonal coefficients are equal to 1 -const unsigned int NullLowerBit = 0x200; ///< means the strictly triangular lower part is 0 -const unsigned int NullUpperBit = 0x400; ///< means the strictly triangular upper part is 0 +const unsigned int SelfAdjointBit = 0x100; ///< means the matrix is selfadjoint (M=M*). +const unsigned int UpperTriangularBit = 0x200; ///< means the strictly triangular lower part is 0 +const unsigned int LowerTriangularBit = 0x400; ///< means the strictly triangular upper part is 0 const unsigned int DirectAccessBit = 0x800; ///< means the underlying matrix data can be direclty accessed -const unsigned int TemporaryBit = 0x1000; ///< means the expression should be copied by value when nested - -enum { Upper=NullLowerBit, Lower=NullUpperBit }; -enum { Aligned=0, UnAligned=1 }; +const unsigned int NestByValueBit = 0x1000; ///< means the expression should be copied by value when nested // list of flags that are inherited by default const unsigned int HereditaryBits = RowMajorBit @@ -55,6 +53,22 @@ const unsigned int HereditaryBits = RowMajorBit | EvalBeforeAssigningBit | LargeBit; +// Possible values for the PartType parameter of part() and the ExtractType parameter of extract() +const unsigned int Upper = UpperTriangularBit; +const unsigned int StrictlyUpper = UpperTriangularBit | ZeroDiagBit; +const unsigned int Lower = LowerTriangularBit; +const unsigned int StrictlyLower = LowerTriangularBit | ZeroDiagBit; + +// additional possible values for the PartType parameter of part() +const unsigned int SelfAdjoint = SelfAdjointBit; + +// additional possible values for the ExtractType parameter of extract() +const unsigned int UnitUpper = UpperTriangularBit | UnitDiagBit; +const unsigned int UnitLower = LowerTriangularBit | UnitDiagBit; + + + +enum { Aligned=0, UnAligned=1 }; enum { ConditionalJumpCost = 5 }; enum CornerType { TopLeft, TopRight, BottomLeft, BottomRight }; enum DirectionType { Vertical, Horizontal }; diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index 85892eb3a..0ab40c6a4 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -46,9 +46,10 @@ template<typename Lhs, typename Rhs, int EvalMode=ei_product_eval_mode<Lhs,Rhs>: template<typename CoeffsVectorType> class DiagonalMatrix; template<typename MatrixType> class DiagonalCoeffs; template<typename MatrixType> class Map; -// template<typename Derived> class Eval; template<int Direction, typename UnaryOp, typename MatrixType> class PartialRedux; -template<int Mode, typename MatrixType> class Triangular; +template<typename MatrixType, unsigned int Mode> class Part; +template<typename MatrixType, unsigned int Mode> class Extract; + template<typename Scalar> struct ei_scalar_sum_op; template<typename Scalar> struct ei_scalar_difference_op; diff --git a/Eigen/src/Core/util/Macros.h b/Eigen/src/Core/util/Macros.h index 93694df09..b511da16d 100644 --- a/Eigen/src/Core/util/Macros.h +++ b/Eigen/src/Core/util/Macros.h @@ -135,7 +135,7 @@ typedef typename Eigen::NumTraits<Scalar>::Real RealScalar; \ typedef typename Base::PacketScalar PacketScalar; \ typedef typename Eigen::ei_nested<Derived>::type Nested; \ typedef typename Eigen::ei_eval<Derived>::type Eval; \ -typedef Eigen::Flagged<Derived, Eigen::TemporaryBit, 0> Temporary; \ +typedef Eigen::Flagged<Derived, NestByValueBit, 0> Temporary; \ enum { RowsAtCompileTime = Base::RowsAtCompileTime, \ ColsAtCompileTime = Base::ColsAtCompileTime, \ MaxRowsAtCompileTime = Base::MaxRowsAtCompileTime, \ diff --git a/Eigen/src/Core/util/Meta.h b/Eigen/src/Core/util/Meta.h index 8d9e4c00a..6c0730531 100644 --- a/Eigen/src/Core/util/Meta.h +++ b/Eigen/src/Core/util/Meta.h @@ -192,15 +192,11 @@ template<typename T> struct ei_unref<T&> { typedef T type; }; template<typename T> struct ei_unconst { typedef T type; }; template<typename T> struct ei_unconst<const T> { typedef T type; }; -template<typename T> struct ei_is_temporary -{ - enum { ret = int(ei_traits<T>::Flags) & TemporaryBit }; -}; template<typename T, int n=1> struct ei_nested { typedef typename ei_meta_if< - ei_is_temporary<T>::ret, + ei_traits<T>::Flags & NestByValueBit, T, typename ei_meta_if< int(ei_traits<T>::Flags) & EvalBeforeNestingBit diff --git a/Eigen/src/LU/Determinant.h b/Eigen/src/LU/Determinant.h index 1675efc10..e6a934147 100644 --- a/Eigen/src/LU/Determinant.h +++ b/Eigen/src/LU/Determinant.h @@ -71,11 +71,11 @@ template<typename Derived> typename ei_traits<Derived>::Scalar MatrixBase<Derived>::determinant() const { assert(rows() == cols()); - if (Derived::Flags & (NullLowerBit | NullUpperBit)) + if (Derived::Flags & (UpperTriangularBit | LowerTriangularBit)) { if (Derived::Flags & UnitDiagBit) return 1; - else if (Derived::Flags & NullDiagBit) + else if (Derived::Flags & ZeroDiagBit) return 0; else return derived().diagonal().redux(ei_scalar_product_op<Scalar>()); diff --git a/Eigen/src/QR/QR.h b/Eigen/src/QR/QR.h index d42153eb9..782f43a04 100644 --- a/Eigen/src/QR/QR.h +++ b/Eigen/src/QR/QR.h @@ -111,7 +111,7 @@ template<typename MatrixType> typename QR<MatrixType>::RMatrixType QR<MatrixType>::matrixR(void) const { int cols = m_qr.cols(); - RMatrixType res = m_qr.block(0,0,cols,cols).upperWithNullDiag(); + RMatrixType res = m_qr.block(0,0,cols,cols).strictlyUpper(); res.diagonal() = m_norms; return res; } diff --git a/test/triangular.cpp b/test/triangular.cpp index f384c4d50..01d4ecf84 100644 --- a/test/triangular.cpp +++ b/test/triangular.cpp @@ -47,8 +47,8 @@ template<typename MatrixType> void triangular(const MatrixType& m) v2 = VectorType::random(rows), vzero = VectorType::zero(rows); - MatrixType m1up = m1.upper(); - MatrixType m2up = m2.upper(); + MatrixType m1up = m1.template extract<Eigen::Upper>(); + MatrixType m2up = m2.template extract<Eigen::Upper>(); if (rows*cols>1) { @@ -62,26 +62,26 @@ template<typename MatrixType> void triangular(const MatrixType& m) // test overloaded operator+= r1.setZero(); r2.setZero(); - r1.upper() += m1; + r1.template part<Eigen::Upper>() += m1; r2 += m1up; VERIFY_IS_APPROX(r1,r2); // test overloaded operator= m1.setZero(); - m1.upper() = (m2.transpose() * m2).lazy(); + m1.template part<Eigen::Upper>() = (m2.transpose() * m2).lazy(); m3 = m2.transpose() * m2; - VERIFY_IS_APPROX(m3.lower().transpose(), m1); + VERIFY_IS_APPROX(m3.template extract<Eigen::Lower>().transpose(), m1); // test overloaded operator= m1.setZero(); - m1.lower() = (m2.transpose() * m2).lazy(); - VERIFY_IS_APPROX(m3.lower(), m1); + m1.template part<Eigen::Lower>() = (m2.transpose() * m2).lazy(); + VERIFY_IS_APPROX(m3.template extract<Eigen::Lower>(), m1); // test back and forward subsitution m1 = MatrixType::random(rows, cols); - VERIFY_IS_APPROX(m1.upper() * (m1.upper().inverseProduct(m2)), m2); - VERIFY_IS_APPROX(m1.lower() * (m1.lower().inverseProduct(m2)), m2); - VERIFY((m1.upper() * m2.upper()).isUpper()); + VERIFY_IS_APPROX(m1.template extract<Eigen::Upper>() * (m1.template extract<Eigen::Upper>().inverseProduct(m2)), m2); + VERIFY_IS_APPROX(m1.template extract<Eigen::Lower>() * (m1.template extract<Eigen::Lower>().inverseProduct(m2)), m2); + VERIFY((m1.template extract<Eigen::Upper>() * m2.template extract<Eigen::Upper>()).isUpper()); } |