diff options
58 files changed, 591 insertions, 563 deletions
diff --git a/Eigen/Core b/Eigen/Core index 67d55fee3..1f71dfdb5 100644 --- a/Eigen/Core +++ b/Eigen/Core @@ -23,10 +23,6 @@ // License and a copy of the GNU General Public License along with // Eigen. If not, see <http://www.gnu.org/licenses/>. -#ifdef EIGEN2_SUPPORT -#include "Eigen2Support" -#endif - #ifndef EIGEN_CORE_H #define EIGEN_CORE_H @@ -225,5 +221,9 @@ struct Dense {}; #include "src/Core/util/EnableMSVCWarnings.h" +#ifdef EIGEN2_SUPPORT +#include "Eigen2Support" +#endif + #endif // EIGEN_CORE_H /* vim: set filetype=cpp et sw=2 ts=2 ai: */ diff --git a/Eigen/Eigen2Support b/Eigen/Eigen2Support index b95b51dec..ccba4ff30 100644 --- a/Eigen/Eigen2Support +++ b/Eigen/Eigen2Support @@ -25,14 +25,13 @@ #ifndef EIGEN2SUPPORT_H #define EIGEN2SUPPORT_H -#include "src/Core/util/DisableMSVCWarnings.h" +#include "Array" -#ifndef EIGEN2_SUPPORT -#define EIGEN2_SUPPORT +#if (!defined(EIGEN2_SUPPORT)) || (!defined(EIGEN_CORE_H)) +#error Eigen2 support must be enabled by defining EIGEN2_SUPPORT before any other Eigen header #endif -#include "Core" -#include "Array" +#include "src/Core/util/DisableMSVCWarnings.h" namespace Eigen { @@ -40,8 +39,9 @@ namespace Eigen { * This module provides a couple of deprecated functions improving the compatibility with Eigen2. * * \code - * #include <Eigen/Eigen2Support> + * #define EIGEN2_SUPPORT * \endcode + * */ #include "src/Eigen2Support/Flagged.h" diff --git a/Eigen/src/Array/Reverse.h b/Eigen/src/Array/Reverse.h index a8d8310a9..7d2c34816 100644 --- a/Eigen/src/Array/Reverse.h +++ b/Eigen/src/Array/Reverse.h @@ -59,9 +59,7 @@ struct ei_traits<Reverse<MatrixType, Direction> > LinearAccess = ( (Direction==BothDirections) && (int(_MatrixTypeNested::Flags)&PacketAccessBit) ) ? LinearAccessBit : 0, - Flags = (int(_MatrixTypeNested::Flags) & (HereditaryBits | PacketAccessBit | LinearAccess)) - | (int(_MatrixTypeNested::Flags)&UpperTriangularBit ? LowerTriangularBit : 0) - | (int(_MatrixTypeNested::Flags)&LowerTriangularBit ? UpperTriangularBit : 0), + Flags = (int(_MatrixTypeNested::Flags) & (HereditaryBits | PacketAccessBit | LinearAccess)), CoeffReadCost = _MatrixTypeNested::CoeffReadCost }; diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index 4fb6d3d2c..8b13d8e2e 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -80,7 +80,7 @@ template<typename _MatrixType> class LDLT } /** \returns the lower triangular matrix L */ - inline TriangularView<MatrixType, UnitLowerTriangular> matrixL(void) const + inline TriangularView<MatrixType, UnitLower> matrixL(void) const { ei_assert(m_isInitialized && "LDLT is not initialized."); return m_matrix; @@ -301,13 +301,13 @@ bool LDLT<MatrixType>::solveInPlace(MatrixBase<Derived> &bAndX) const // y = L^-1 z //matrixL().solveInPlace(bAndX); - m_matrix.template triangularView<UnitLowerTriangular>().solveInPlace(bAndX); + m_matrix.template triangularView<UnitLower>().solveInPlace(bAndX); // w = D^-1 y bAndX = m_matrix.diagonal().asDiagonal().inverse() * bAndX; // u = L^-T w - m_matrix.adjoint().template triangularView<UnitUpperTriangular>().solveInPlace(bAndX); + m_matrix.adjoint().template triangularView<UnitUpper>().solveInPlace(bAndX); // x = P^T u for (int i = size-1; i >= 0; --i) bAndX.row(m_transpositions.coeff(i)).swap(bAndX.row(i)); diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h index 02645b23f..8a149a316 100644 --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h @@ -148,7 +148,7 @@ template<typename _MatrixType, int _UpLo> class LLT // forward declaration (defined at the end of this file) template<int UpLo> struct ei_llt_inplace; -template<> struct ei_llt_inplace<LowerTriangular> +template<> struct ei_llt_inplace<Lower> { template<typename MatrixType> static bool unblocked(MatrixType& mat) @@ -198,47 +198,47 @@ template<> struct ei_llt_inplace<LowerTriangular> Block<MatrixType,Dynamic,Dynamic> A22(m,k+bs,k+bs,rs,rs); if(!unblocked(A11)) return false; - if(rs>0) A11.adjoint().template triangularView<UpperTriangular>().template solveInPlace<OnTheRight>(A21); - if(rs>0) A22.template selfadjointView<LowerTriangular>().rankUpdate(A21,-1); // bottleneck + if(rs>0) A11.adjoint().template triangularView<Upper>().template solveInPlace<OnTheRight>(A21); + if(rs>0) A22.template selfadjointView<Lower>().rankUpdate(A21,-1); // bottleneck } return true; } }; -template<> struct ei_llt_inplace<UpperTriangular> +template<> struct ei_llt_inplace<Upper> { template<typename MatrixType> static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat) { Transpose<MatrixType> matt(mat); - return ei_llt_inplace<LowerTriangular>::unblocked(matt); + return ei_llt_inplace<Lower>::unblocked(matt); } template<typename MatrixType> static EIGEN_STRONG_INLINE bool blocked(MatrixType& mat) { Transpose<MatrixType> matt(mat); - return ei_llt_inplace<LowerTriangular>::blocked(matt); + return ei_llt_inplace<Lower>::blocked(matt); } }; -template<typename MatrixType> struct LLT_Traits<MatrixType,LowerTriangular> +template<typename MatrixType> struct LLT_Traits<MatrixType,Lower> { - typedef TriangularView<MatrixType, LowerTriangular> MatrixL; - typedef TriangularView<typename MatrixType::AdjointReturnType, UpperTriangular> MatrixU; + typedef TriangularView<MatrixType, Lower> MatrixL; + typedef TriangularView<typename MatrixType::AdjointReturnType, Upper> MatrixU; inline static MatrixL getL(const MatrixType& m) { return m; } inline static MatrixU getU(const MatrixType& m) { return m.adjoint(); } static bool inplace_decomposition(MatrixType& m) - { return ei_llt_inplace<LowerTriangular>::blocked(m); } + { return ei_llt_inplace<Lower>::blocked(m); } }; -template<typename MatrixType> struct LLT_Traits<MatrixType,UpperTriangular> +template<typename MatrixType> struct LLT_Traits<MatrixType,Upper> { - typedef TriangularView<typename MatrixType::AdjointReturnType, LowerTriangular> MatrixL; - typedef TriangularView<MatrixType, UpperTriangular> MatrixU; + typedef TriangularView<typename MatrixType::AdjointReturnType, Lower> MatrixL; + typedef TriangularView<MatrixType, Upper> MatrixU; inline static MatrixL getL(const MatrixType& m) { return m.adjoint(); } inline static MatrixU getU(const MatrixType& m) { return m; } static bool inplace_decomposition(MatrixType& m) - { return ei_llt_inplace<UpperTriangular>::blocked(m); } + { return ei_llt_inplace<Upper>::blocked(m); } }; /** Computes / recomputes the Cholesky decomposition A = LL^* = U^*U of \a matrix diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 4c30f30ad..7569dab58 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -405,13 +405,6 @@ template<typename Derived> class MatrixBase template<int Size> const VectorBlock<Derived,Size> start() const; template<int Size> VectorBlock<Derived,Size> end(); template<int Size> const VectorBlock<Derived,Size> end() const; - - template<typename OtherDerived> - typename ei_plain_matrix_type_column_major<OtherDerived>::type - solveTriangular(const MatrixBase<OtherDerived>& other) const; - - template<typename OtherDerived> - void solveTriangularInPlace(const MatrixBase<OtherDerived>& other) const; #endif protected: diff --git a/Eigen/src/Core/SelfAdjointView.h b/Eigen/src/Core/SelfAdjointView.h index 9ab8cb2c1..6d01ee495 100644 --- a/Eigen/src/Core/SelfAdjointView.h +++ b/Eigen/src/Core/SelfAdjointView.h @@ -31,7 +31,7 @@ * \brief Expression of a selfadjoint matrix from a triangular part of a dense matrix * * \param MatrixType the type of the dense matrix storing the coefficients - * \param TriangularPart can be either \c LowerTriangular or \c UpperTriangular + * \param TriangularPart can be either \c Lower or \c Upper * * This class is an expression of a sefladjoint matrix from a triangular part of a matrix * with given dense storage of the coefficients. It is the return type of MatrixBase::selfadjointView() @@ -46,7 +46,7 @@ struct ei_traits<SelfAdjointView<MatrixType, TriangularPart> > : ei_traits<Matri typedef typename ei_unref<MatrixTypeNested>::type _MatrixTypeNested; typedef MatrixType ExpressionType; enum { - Mode = TriangularPart | SelfAdjointBit, + Mode = TriangularPart | SelfAdjoint, Flags = _MatrixTypeNested::Flags & (HereditaryBits) & (~(PacketAccessBit | DirectAccessBit | LinearAccessBit)), // FIXME these flags should be preserved CoeffReadCost = _MatrixTypeNested::CoeffReadCost diff --git a/Eigen/src/Core/SolveTriangular.h b/Eigen/src/Core/SolveTriangular.h index 9dc019d17..73cf77387 100644 --- a/Eigen/src/Core/SolveTriangular.h +++ b/Eigen/src/Core/SolveTriangular.h @@ -57,7 +57,7 @@ struct ei_triangular_solver_selector<Lhs,Rhs,OnTheLeft,Mode,NoUnrolling,RowMajor typedef ei_blas_traits<Lhs> LhsProductTraits; typedef typename LhsProductTraits::ExtractType ActualLhsType; enum { - IsLowerTriangular = ((Mode&LowerTriangularBit)==LowerTriangularBit) + IsLower = ((Mode&Lower)==Lower) }; static void run(const Lhs& lhs, Rhs& other) { @@ -65,20 +65,20 @@ struct ei_triangular_solver_selector<Lhs,Rhs,OnTheLeft,Mode,NoUnrolling,RowMajor ActualLhsType actualLhs = LhsProductTraits::extract(lhs); const int size = lhs.cols(); - for(int pi=IsLowerTriangular ? 0 : size; - IsLowerTriangular ? pi<size : pi>0; - IsLowerTriangular ? pi+=PanelWidth : pi-=PanelWidth) + for(int pi=IsLower ? 0 : size; + IsLower ? pi<size : pi>0; + IsLower ? pi+=PanelWidth : pi-=PanelWidth) { - int actualPanelWidth = std::min(IsLowerTriangular ? size - pi : pi, PanelWidth); + int actualPanelWidth = std::min(IsLower ? size - pi : pi, PanelWidth); - int r = IsLowerTriangular ? pi : size - pi; // remaining size + int r = IsLower ? pi : size - pi; // remaining size if (r > 0) { // let's directly call the low level product function because: // 1 - it is faster to compile // 2 - it is slighlty faster at runtime - int startRow = IsLowerTriangular ? pi : pi-actualPanelWidth; - int startCol = IsLowerTriangular ? 0 : pi; + int startRow = IsLower ? pi : pi-actualPanelWidth; + int startCol = IsLower ? 0 : pi; VectorBlock<Rhs,Dynamic> target(other,startRow,actualPanelWidth); ei_cache_friendly_product_rowmajor_times_vector<LhsProductTraits::NeedToConjugate,false>( @@ -89,12 +89,12 @@ struct ei_triangular_solver_selector<Lhs,Rhs,OnTheLeft,Mode,NoUnrolling,RowMajor for(int k=0; k<actualPanelWidth; ++k) { - int i = IsLowerTriangular ? pi+k : pi-k-1; - int s = IsLowerTriangular ? pi : i+1; + int i = IsLower ? pi+k : pi-k-1; + int s = IsLower ? pi : i+1; if (k>0) other.coeffRef(i) -= (lhs.row(i).segment(s,k).transpose().cwiseProduct(other.segment(s,k))).sum(); - if(!(Mode & UnitDiagBit)) + if(!(Mode & UnitDiag)) other.coeffRef(i) /= lhs.coeff(i,i); } } @@ -111,7 +111,7 @@ struct ei_triangular_solver_selector<Lhs,Rhs,OnTheLeft,Mode,NoUnrolling,ColMajor typedef typename LhsProductTraits::ExtractType ActualLhsType; enum { PacketSize = ei_packet_traits<Scalar>::size, - IsLowerTriangular = ((Mode&LowerTriangularBit)==LowerTriangularBit) + IsLower = ((Mode&Lower)==Lower) }; static void run(const Lhs& lhs, Rhs& other) @@ -120,26 +120,26 @@ struct ei_triangular_solver_selector<Lhs,Rhs,OnTheLeft,Mode,NoUnrolling,ColMajor ActualLhsType actualLhs = LhsProductTraits::extract(lhs); const int size = lhs.cols(); - for(int pi=IsLowerTriangular ? 0 : size; - IsLowerTriangular ? pi<size : pi>0; - IsLowerTriangular ? pi+=PanelWidth : pi-=PanelWidth) + for(int pi=IsLower ? 0 : size; + IsLower ? pi<size : pi>0; + IsLower ? pi+=PanelWidth : pi-=PanelWidth) { - int actualPanelWidth = std::min(IsLowerTriangular ? size - pi : pi, PanelWidth); - int startBlock = IsLowerTriangular ? pi : pi-actualPanelWidth; - int endBlock = IsLowerTriangular ? pi + actualPanelWidth : 0; + int actualPanelWidth = std::min(IsLower ? size - pi : pi, PanelWidth); + int startBlock = IsLower ? pi : pi-actualPanelWidth; + int endBlock = IsLower ? pi + actualPanelWidth : 0; for(int k=0; k<actualPanelWidth; ++k) { - int i = IsLowerTriangular ? pi+k : pi-k-1; - if(!(Mode & UnitDiagBit)) + int i = IsLower ? pi+k : pi-k-1; + if(!(Mode & UnitDiag)) other.coeffRef(i) /= lhs.coeff(i,i); int r = actualPanelWidth - k - 1; // remaining size - int s = IsLowerTriangular ? i+1 : i-r; + int s = IsLower ? i+1 : i-r; if (r>0) other.segment(s,r) -= other.coeffRef(i) * Block<Lhs,Dynamic,1>(lhs, s, i, r, 1); } - int r = IsLowerTriangular ? size - endBlock : startBlock; // remaining size + int r = IsLower ? size - endBlock : startBlock; // remaining size if (r > 0) { // let's directly call the low level product function because: @@ -198,16 +198,16 @@ struct ei_triangular_solver_unroller; template<typename Lhs, typename Rhs, int Mode, int Index, int Size> struct ei_triangular_solver_unroller<Lhs,Rhs,Mode,Index,Size,false> { enum { - IsLowerTriangular = ((Mode&LowerTriangularBit)==LowerTriangularBit), - I = IsLowerTriangular ? Index : Size - Index - 1, - S = IsLowerTriangular ? 0 : I+1 + IsLower = ((Mode&Lower)==Lower), + I = IsLower ? Index : Size - Index - 1, + S = IsLower ? 0 : I+1 }; static void run(const Lhs& lhs, Rhs& rhs) { if (Index>0) rhs.coeffRef(I) -= ((lhs.row(I).template segment<Index>(S).transpose()).cwiseProduct(rhs.template segment<Index>(S))).sum(); - if(!(Mode & UnitDiagBit)) + if(!(Mode & UnitDiag)) rhs.coeffRef(I) /= lhs.coeff(I,I); ei_triangular_solver_unroller<Lhs,Rhs,Mode,Index+1,Size>::run(lhs,rhs); @@ -245,8 +245,8 @@ void TriangularView<MatrixType,Mode>::solveInPlace(const MatrixBase<OtherDerived OtherDerived& other = _other.const_cast_derived(); ei_assert(cols() == rows()); ei_assert( (Side==OnTheLeft && cols() == other.rows()) || (Side==OnTheRight && cols() == other.cols()) ); - ei_assert(!(Mode & ZeroDiagBit)); - ei_assert(Mode & (UpperTriangularBit|LowerTriangularBit)); + ei_assert(!(Mode & ZeroDiag)); + ei_assert(Mode & (Upper|Lower)); enum { copy = ei_traits<OtherDerived>::Flags & RowMajorBit && OtherDerived::IsVectorAtCompileTime }; typedef typename ei_meta_if<copy, diff --git a/Eigen/src/Core/Transpose.h b/Eigen/src/Core/Transpose.h index 80b056527..35b8b2ed3 100644 --- a/Eigen/src/Core/Transpose.h +++ b/Eigen/src/Core/Transpose.h @@ -49,10 +49,7 @@ struct ei_traits<Transpose<MatrixType> > : ei_traits<MatrixType> ColsAtCompileTime = MatrixType::RowsAtCompileTime, MaxRowsAtCompileTime = MatrixType::MaxColsAtCompileTime, MaxColsAtCompileTime = MatrixType::MaxRowsAtCompileTime, - Flags = ((int(_MatrixTypeNested::Flags) ^ RowMajorBit) - & ~(LowerTriangularBit | UpperTriangularBit)) - | (int(_MatrixTypeNested::Flags)&UpperTriangularBit ? LowerTriangularBit : 0) - | (int(_MatrixTypeNested::Flags)&LowerTriangularBit ? UpperTriangularBit : 0), + Flags = (int(_MatrixTypeNested::Flags) ^ RowMajorBit), CoeffReadCost = _MatrixTypeNested::CoeffReadCost }; }; @@ -229,7 +226,7 @@ struct ei_inplace_transpose_selector; template<typename MatrixType> struct ei_inplace_transpose_selector<MatrixType,true> { // square matrix static void run(MatrixType& m) { - m.template triangularView<StrictlyUpperTriangular>().swap(m.transpose()); + m.template triangularView<StrictlyUpper>().swap(m.transpose()); } }; @@ -237,7 +234,7 @@ template<typename MatrixType> struct ei_inplace_transpose_selector<MatrixType,false> { // non square matrix static void run(MatrixType& m) { if (m.rows()==m.cols()) - m.template triangularView<StrictlyUpperTriangular>().swap(m.transpose()); + m.template triangularView<StrictlyUpper>().swap(m.transpose()); else m = m.transpose().eval(); } diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h index 62d800fef..8411b0546 100644 --- a/Eigen/src/Core/TriangularMatrix.h +++ b/Eigen/src/Core/TriangularMatrix.h @@ -46,7 +46,7 @@ template<typename Derived> class TriangularBase : public AnyMatrixBase<Derived> }; typedef typename ei_traits<Derived>::Scalar Scalar; - inline TriangularBase() { ei_assert(ei_are_flags_consistent<Mode>::ret); } + inline TriangularBase() { ei_assert(!((Mode&UnitDiag) && (Mode&ZeroDiag))); } inline int rows() const { return derived().rows(); } inline int cols() const { return derived().cols(); } @@ -89,10 +89,10 @@ template<typename Derived> class TriangularBase : public AnyMatrixBase<Derived> void check_coordinates(int row, int col) { ei_assert(col>=0 && col<cols() && row>=0 && row<rows()); - ei_assert( (Mode==UpperTriangular && col>=row) - || (Mode==LowerTriangular && col<=row) - || ((Mode==StrictlyUpperTriangular || Mode==UnitUpperTriangular) && col>row) - || ((Mode==StrictlyLowerTriangular || Mode==UnitLowerTriangular) && col<row)); + ei_assert( (Mode==Upper && col>=row) + || (Mode==Lower && col<=row) + || ((Mode==StrictlyUpper || Mode==UnitUpper) && col>row) + || ((Mode==StrictlyLower || Mode==UnitLower) && col<row)); } #ifdef EIGEN_INTERNAL_DEBUGGING @@ -111,10 +111,10 @@ template<typename Derived> class TriangularBase : public AnyMatrixBase<Derived> * \brief Base class for triangular part in a matrix * * \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 UpperTriangular, - * LowerTriangular, UpperSelfadjoint, or LowerSelfadjoint. This is in fact a bit field; - * it must have either UpperBit or LowerBit, and additionnaly it may have either - * TraingularBit or SelfadjointBit. + * \param Mode the kind of triangular matrix expression to construct. Can be Upper, + * Lower, UpperSelfadjoint, or LowerSelfadjoint. This is in fact a bit field; + * it must have either Upper or Lower, and additionnaly it may have either + * UnitDiag or Selfadjoint. * * This class represents a triangular part of a matrix, not necessarily square. Strictly speaking, for rectangular * matrices one should speak ok "trapezoid" parts. This class is the return type @@ -154,9 +154,10 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView enum { Mode = _Mode, - TransposeMode = (Mode & UpperTriangularBit ? LowerTriangularBit : 0) - | (Mode & LowerTriangularBit ? UpperTriangularBit : 0) - | (Mode & (ZeroDiagBit | UnitDiagBit)) + TransposeMode = (Mode & Upper ? Lower : 0) + | (Mode & Lower ? Upper : 0) + | (Mode & (UnitDiag)) + | (Mode & (ZeroDiag)) }; inline TriangularView(const MatrixType& matrix) : m_matrix(matrix) @@ -283,12 +284,12 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView const SelfAdjointView<_MatrixTypeNested,Mode> selfadjointView() const { - EIGEN_STATIC_ASSERT((Mode&UnitDiagBit)==0,PROGRAMMING_ERROR); + EIGEN_STATIC_ASSERT((Mode&UnitDiag)==0,PROGRAMMING_ERROR); return SelfAdjointView<_MatrixTypeNested,Mode>(m_matrix); } SelfAdjointView<_MatrixTypeNested,Mode> selfadjointView() { - EIGEN_STATIC_ASSERT((Mode&UnitDiagBit)==0,PROGRAMMING_ERROR); + EIGEN_STATIC_ASSERT((Mode&UnitDiag)==0,PROGRAMMING_ERROR); return SelfAdjointView<_MatrixTypeNested,Mode>(m_matrix); } @@ -304,6 +305,16 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView TriangularView<SwapWrapper<MatrixType>,Mode>(const_cast<MatrixType&>(m_matrix)).lazyAssign(other.derived()); } + Scalar determinant() const + { + if (Mode & UnitDiag) + return 1; + else if (Mode & ZeroDiag) + return 0; + else + return m_matrix.diagonal().prod(); + } + protected: const MatrixTypeNested m_matrix; @@ -325,19 +336,19 @@ struct ei_triangular_assignment_selector { ei_triangular_assignment_selector<Derived1, Derived2, Mode, UnrollCount-1, ClearOpposite>::run(dst, src); - ei_assert( Mode == UpperTriangular || Mode == LowerTriangular - || Mode == StrictlyUpperTriangular || Mode == StrictlyLowerTriangular - || Mode == UnitUpperTriangular || Mode == UnitLowerTriangular); - if((Mode == UpperTriangular && row <= col) - || (Mode == LowerTriangular && row >= col) - || (Mode == StrictlyUpperTriangular && row < col) - || (Mode == StrictlyLowerTriangular && row > col) - || (Mode == UnitUpperTriangular && row < col) - || (Mode == UnitLowerTriangular && row > col)) + ei_assert( Mode == Upper || Mode == Lower + || Mode == StrictlyUpper || Mode == StrictlyLower + || Mode == UnitUpper || Mode == UnitLower); + if((Mode == Upper && row <= col) + || (Mode == Lower && row >= col) + || (Mode == StrictlyUpper && row < col) + || (Mode == StrictlyLower && row > col) + || (Mode == UnitUpper && row < col) + || (Mode == UnitLower && row > col)) dst.copyCoeff(row, col, src); else if(ClearOpposite) { - if (Mode&UnitDiagBit && row==col) + if (Mode&UnitDiag && row==col) dst.coeffRef(row, col) = 1; else dst.coeffRef(row, col) = 0; @@ -353,7 +364,7 @@ struct ei_triangular_assignment_selector<Derived1, Derived2, Mode, 0, ClearOppos }; template<typename Derived1, typename Derived2, bool ClearOpposite> -struct ei_triangular_assignment_selector<Derived1, Derived2, UpperTriangular, Dynamic, ClearOpposite> +struct ei_triangular_assignment_selector<Derived1, Derived2, Upper, Dynamic, ClearOpposite> { inline static void run(Derived1 &dst, const Derived2 &src) { @@ -370,7 +381,7 @@ struct ei_triangular_assignment_selector<Derived1, Derived2, UpperTriangular, Dy }; template<typename Derived1, typename Derived2, bool ClearOpposite> -struct ei_triangular_assignment_selector<Derived1, Derived2, LowerTriangular, Dynamic, ClearOpposite> +struct ei_triangular_assignment_selector<Derived1, Derived2, Lower, Dynamic, ClearOpposite> { inline static void run(Derived1 &dst, const Derived2 &src) { @@ -387,7 +398,7 @@ struct ei_triangular_assignment_selector<Derived1, Derived2, LowerTriangular, Dy }; template<typename Derived1, typename Derived2, bool ClearOpposite> -struct ei_triangular_assignment_selector<Derived1, Derived2, StrictlyUpperTriangular, Dynamic, ClearOpposite> +struct ei_triangular_assignment_selector<Derived1, Derived2, StrictlyUpper, Dynamic, ClearOpposite> { inline static void run(Derived1 &dst, const Derived2 &src) { @@ -404,7 +415,7 @@ struct ei_triangular_assignment_selector<Derived1, Derived2, StrictlyUpperTriang }; template<typename Derived1, typename Derived2, bool ClearOpposite> -struct ei_triangular_assignment_selector<Derived1, Derived2, StrictlyLowerTriangular, Dynamic, ClearOpposite> +struct ei_triangular_assignment_selector<Derived1, Derived2, StrictlyLower, Dynamic, ClearOpposite> { inline static void run(Derived1 &dst, const Derived2 &src) { @@ -421,7 +432,7 @@ struct ei_triangular_assignment_selector<Derived1, Derived2, StrictlyLowerTriang }; template<typename Derived1, typename Derived2, bool ClearOpposite> -struct ei_triangular_assignment_selector<Derived1, Derived2, UnitUpperTriangular, Dynamic, ClearOpposite> +struct ei_triangular_assignment_selector<Derived1, Derived2, UnitUpper, Dynamic, ClearOpposite> { inline static void run(Derived1 &dst, const Derived2 &src) { @@ -440,7 +451,7 @@ struct ei_triangular_assignment_selector<Derived1, Derived2, UnitUpperTriangular } }; template<typename Derived1, typename Derived2, bool ClearOpposite> -struct ei_triangular_assignment_selector<Derived1, Derived2, UnitLowerTriangular, Dynamic, ClearOpposite> +struct ei_triangular_assignment_selector<Derived1, Derived2, UnitLower, Dynamic, ClearOpposite> { inline static void run(Derived1 &dst, const Derived2 &src) { @@ -590,8 +601,8 @@ EIGEN_DEPRECATED TriangularView<Derived, Mode> MatrixBase<Derived>::part() /** \nonstableyet * \returns an expression of a triangular view extracted from the current matrix * - * The parameter \a Mode can have the following values: \c UpperTriangular, \c StrictlyUpperTriangular, \c UnitUpperTriangular, - * \c LowerTriangular, \c StrictlyLowerTriangular, \c UnitLowerTriangular. + * The parameter \a Mode can have the following values: \c Upper, \c StrictlyUpper, \c UnitUpper, + * \c Lower, \c StrictlyLower, \c UnitLower. * * Example: \include MatrixBase_extract.cpp * Output: \verbinclude MatrixBase_extract.out @@ -616,22 +627,22 @@ const TriangularView<Derived, Mode> MatrixBase<Derived>::triangularView() const /** \returns true if *this is approximately equal to an upper triangular matrix, * within the precision given by \a prec. * - * \sa isLowerTriangular(), extract(), part(), marked() + * \sa isLower(), extract(), part(), marked() */ template<typename Derived> bool MatrixBase<Derived>::isUpperTriangular(RealScalar prec) const { - RealScalar maxAbsOnUpperTriangularPart = static_cast<RealScalar>(-1); + RealScalar maxAbsOnUpperPart = static_cast<RealScalar>(-1); for(int j = 0; j < cols(); ++j) { int maxi = std::min(j, rows()-1); for(int i = 0; i <= maxi; ++i) { RealScalar absValue = ei_abs(coeff(i,j)); - if(absValue > maxAbsOnUpperTriangularPart) maxAbsOnUpperTriangularPart = absValue; + if(absValue > maxAbsOnUpperPart) maxAbsOnUpperPart = absValue; } } - RealScalar threshold = maxAbsOnUpperTriangularPart * prec; + RealScalar threshold = maxAbsOnUpperPart * prec; for(int j = 0; j < cols(); ++j) for(int i = j+1; i < rows(); ++i) if(ei_abs(coeff(i, j)) > threshold) return false; @@ -641,19 +652,19 @@ bool MatrixBase<Derived>::isUpperTriangular(RealScalar prec) const /** \returns true if *this is approximately equal to a lower triangular matrix, * within the precision given by \a prec. * - * \sa isUpperTriangular(), extract(), part(), marked() + * \sa isUpper(), extract(), part(), marked() */ template<typename Derived> bool MatrixBase<Derived>::isLowerTriangular(RealScalar prec) const { - RealScalar maxAbsOnLowerTriangularPart = static_cast<RealScalar>(-1); + RealScalar maxAbsOnLowerPart = static_cast<RealScalar>(-1); for(int j = 0; j < cols(); ++j) for(int i = j; i < rows(); ++i) { RealScalar absValue = ei_abs(coeff(i,j)); - if(absValue > maxAbsOnLowerTriangularPart) maxAbsOnLowerTriangularPart = absValue; + if(absValue > maxAbsOnLowerPart) maxAbsOnLowerPart = absValue; } - RealScalar threshold = maxAbsOnLowerTriangularPart * prec; + RealScalar threshold = maxAbsOnLowerPart * prec; for(int j = 1; j < cols(); ++j) { int maxi = std::min(j, rows()-1); diff --git a/Eigen/src/Core/products/SelfadjointMatrixMatrix.h b/Eigen/src/Core/products/SelfadjointMatrixMatrix.h index 35efa752e..ac072e189 100644 --- a/Eigen/src/Core/products/SelfadjointMatrixMatrix.h +++ b/Eigen/src/Core/products/SelfadjointMatrixMatrix.h @@ -368,10 +368,10 @@ struct SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,RhsMode,false> SelfadjointProductMatrix(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) {} enum { - LhsUpLo = LhsMode&(UpperTriangularBit|LowerTriangularBit), - LhsIsSelfAdjoint = (LhsMode&SelfAdjointBit)==SelfAdjointBit, - RhsUpLo = RhsMode&(UpperTriangularBit|LowerTriangularBit), - RhsIsSelfAdjoint = (RhsMode&SelfAdjointBit)==SelfAdjointBit + LhsUpLo = LhsMode&(Upper|Lower), + LhsIsSelfAdjoint = (LhsMode&SelfAdjoint)==SelfAdjoint, + RhsUpLo = RhsMode&(Upper|Lower), + RhsIsSelfAdjoint = (RhsMode&SelfAdjoint)==SelfAdjoint }; template<typename Dest> void scaleAndAddTo(Dest& dst, Scalar alpha) const @@ -385,12 +385,12 @@ struct SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,RhsMode,false> * RhsBlasTraits::extractScalarFactor(m_rhs); ei_product_selfadjoint_matrix<Scalar, - EIGEN_LOGICAL_XOR(LhsUpLo==UpperTriangular, + EIGEN_LOGICAL_XOR(LhsUpLo==Upper, ei_traits<Lhs>::Flags &RowMajorBit) ? RowMajor : ColMajor, LhsIsSelfAdjoint, - NumTraits<Scalar>::IsComplex && EIGEN_LOGICAL_XOR(LhsUpLo==UpperTriangular,bool(LhsBlasTraits::NeedToConjugate)), - EIGEN_LOGICAL_XOR(RhsUpLo==UpperTriangular, + NumTraits<Scalar>::IsComplex && EIGEN_LOGICAL_XOR(LhsUpLo==Upper,bool(LhsBlasTraits::NeedToConjugate)), + EIGEN_LOGICAL_XOR(RhsUpLo==Upper, ei_traits<Rhs>::Flags &RowMajorBit) ? RowMajor : ColMajor, RhsIsSelfAdjoint, - NumTraits<Scalar>::IsComplex && EIGEN_LOGICAL_XOR(RhsUpLo==UpperTriangular,bool(RhsBlasTraits::NeedToConjugate)), + NumTraits<Scalar>::IsComplex && EIGEN_LOGICAL_XOR(RhsUpLo==Upper,bool(RhsBlasTraits::NeedToConjugate)), ei_traits<Dest>::Flags&RowMajorBit ? RowMajor : ColMajor> ::run( lhs.rows(), rhs.cols(), // sizes diff --git a/Eigen/src/Core/products/SelfadjointMatrixVector.h b/Eigen/src/Core/products/SelfadjointMatrixVector.h index 32b7f220e..1c48208b3 100644 --- a/Eigen/src/Core/products/SelfadjointMatrixVector.h +++ b/Eigen/src/Core/products/SelfadjointMatrixVector.h @@ -42,7 +42,7 @@ static EIGEN_DONT_INLINE void ei_product_selfadjoint_vector( enum { IsRowMajor = StorageOrder==RowMajor ? 1 : 0, - IsLower = UpLo == LowerTriangularBit ? 1 : 0, + IsLower = UpLo == Lower ? 1 : 0, FirstTriangular = IsRowMajor == IsLower }; @@ -170,7 +170,7 @@ struct SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,0,true> EIGEN_PRODUCT_PUBLIC_INTERFACE(SelfadjointProductMatrix) enum { - LhsUpLo = LhsMode&(UpperTriangularBit|LowerTriangularBit) + LhsUpLo = LhsMode&(Upper|Lower) }; SelfadjointProductMatrix(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) {} diff --git a/Eigen/src/Core/products/SelfadjointProduct.h b/Eigen/src/Core/products/SelfadjointProduct.h index f9cdda637..8967f62be 100644 --- a/Eigen/src/Core/products/SelfadjointProduct.h +++ b/Eigen/src/Core/products/SelfadjointProduct.h @@ -47,7 +47,7 @@ struct ei_selfadjoint_product<Scalar,MatStorageOrder, RowMajor, AAT, UpLo> { static EIGEN_STRONG_INLINE void run(int size, int depth, const Scalar* mat, int matStride, Scalar* res, int resStride, Scalar alpha) { - ei_selfadjoint_product<Scalar, MatStorageOrder, ColMajor, !AAT, UpLo==LowerTriangular?UpperTriangular:LowerTriangular> + ei_selfadjoint_product<Scalar, MatStorageOrder, ColMajor, !AAT, UpLo==Lower?Upper:Lower> ::run(size, depth, mat, matStride, res, resStride, alpha); } }; @@ -100,13 +100,13 @@ struct ei_selfadjoint_product<Scalar,MatStorageOrder, ColMajor, AAT, UpLo> // 1 - before the diagonal => processed with gebp or skipped // 2 - the actual_mc x actual_mc symmetric block => processed with a special kernel // 3 - after the diagonal => processed with gebp or skipped - if (UpLo==LowerTriangular) + if (UpLo==Lower) gebp_kernel(res+i2, resStride, blockA, blockB, actual_mc, actual_kc, std::min(size,i2)); ei_sybb_kernel<Scalar, Blocking::mr, Blocking::nr, Conj, UpLo>() (res+resStride*i2 + i2, resStride, blockA, blockB + actual_kc*Blocking::PacketSize*i2, actual_mc, actual_kc); - if (UpLo==UpperTriangular) + if (UpLo==Upper) { int j2 = i2+actual_mc; gebp_kernel(res+resStride*j2+i2, resStride, blockA, blockB+actual_kc*Blocking::PacketSize*j2, actual_mc, actual_kc, std::max(0,size-j2)); @@ -173,7 +173,7 @@ struct ei_sybb_kernel int actualBlockSize = std::min<int>(BlockSize,size - j); const Scalar* actual_b = blockB+j*depth*PacketSize; - if(UpLo==UpperTriangular) + if(UpLo==Upper) gebp_kernel(res+j*resStride, resStride, blockA, actual_b, j, depth, actualBlockSize); // selfadjoint micro block @@ -186,13 +186,13 @@ struct ei_sybb_kernel for(int j1=0; j1<actualBlockSize; ++j1) { Scalar* r = res + (j+j1)*resStride + i; - for(int i1=UpLo==LowerTriangular ? j1 : 0; - UpLo==LowerTriangular ? i1<actualBlockSize : i1<=j1; ++i1) + for(int i1=UpLo==Lower ? j1 : 0; + UpLo==Lower ? i1<actualBlockSize : i1<=j1; ++i1) r[i1] += buffer(i1,j1); } } - if(UpLo==LowerTriangular) + if(UpLo==Lower) { int i = j+actualBlockSize; gebp_kernel(res+j*resStride+i, resStride, blockA+depth*i, actual_b, size-i, depth, actualBlockSize); diff --git a/Eigen/src/Core/products/SelfadjointRank2Update.h b/Eigen/src/Core/products/SelfadjointRank2Update.h index 1c0e503e6..c95de1146 100644 --- a/Eigen/src/Core/products/SelfadjointRank2Update.h +++ b/Eigen/src/Core/products/SelfadjointRank2Update.h @@ -33,7 +33,7 @@ template<typename Scalar, typename UType, typename VType, int UpLo> struct ei_selfadjoint_rank2_update_selector; template<typename Scalar, typename UType, typename VType> -struct ei_selfadjoint_rank2_update_selector<Scalar,UType,VType,LowerTriangular> +struct ei_selfadjoint_rank2_update_selector<Scalar,UType,VType,Lower> { static void run(Scalar* mat, int stride, const UType& u, const VType& v, Scalar alpha) { @@ -48,7 +48,7 @@ struct ei_selfadjoint_rank2_update_selector<Scalar,UType,VType,LowerTriangular> }; template<typename Scalar, typename UType, typename VType> -struct ei_selfadjoint_rank2_update_selector<Scalar,UType,VType,UpperTriangular> +struct ei_selfadjoint_rank2_update_selector<Scalar,UType,VType,Upper> { static void run(Scalar* mat, int stride, const UType& u, const VType& v, Scalar alpha) { @@ -87,7 +87,7 @@ SelfAdjointView<MatrixType,UpLo>& SelfAdjointView<MatrixType,UpLo> ei_selfadjoint_rank2_update_selector<Scalar, typename ei_cleantype<typename ei_conj_expr_if<IsRowMajor ^ UBlasTraits::NeedToConjugate,_ActualUType>::ret>::type, typename ei_cleantype<typename ei_conj_expr_if<IsRowMajor ^ VBlasTraits::NeedToConjugate,_ActualVType>::ret>::type, - (IsRowMajor ? (UpLo==UpperTriangular ? LowerTriangular : UpperTriangular) : UpLo)> + (IsRowMajor ? (UpLo==Upper ? Lower : Upper) : UpLo)> ::run(const_cast<Scalar*>(_expression().data()),_expression().stride(),actualU,actualV,actualAlpha); return *this; diff --git a/Eigen/src/Core/products/TriangularMatrixMatrix.h b/Eigen/src/Core/products/TriangularMatrixMatrix.h index 701ccb644..37617a915 100644 --- a/Eigen/src/Core/products/TriangularMatrixMatrix.h +++ b/Eigen/src/Core/products/TriangularMatrixMatrix.h @@ -75,7 +75,7 @@ struct ei_product_triangular_matrix_matrix<Scalar,Mode,LhsIsTriangular, Scalar alpha) { ei_product_triangular_matrix_matrix<Scalar, - (Mode&UnitDiagBit) | (Mode&UpperTriangular) ? LowerTriangular : UpperTriangular, + (Mode&UnitDiag) | (Mode&Upper) ? Lower : Upper, (!LhsIsTriangular), RhsStorageOrder==RowMajor ? ColMajor : RowMajor, ConjugateRhs, @@ -113,7 +113,7 @@ struct ei_product_triangular_matrix_matrix<Scalar,Mode,true, typedef ei_product_blocking_traits<Scalar> Blocking; enum { SmallPanelWidth = EIGEN_ENUM_MAX(Blocking::mr,Blocking::nr), - IsLowerTriangular = (Mode&LowerTriangular) == LowerTriangular + IsLower = (Mode&Lower) == Lower }; int kc = std::min<int>(Blocking::Max_kc/4,size); // cache block size along the K direction @@ -130,12 +130,12 @@ struct ei_product_triangular_matrix_matrix<Scalar,Mode,true, ei_gemm_pack_lhs<Scalar,Blocking::mr,LhsStorageOrder> pack_lhs; ei_gemm_pack_rhs<Scalar,Blocking::nr,RhsStorageOrder> pack_rhs; - for(int k2=IsLowerTriangular ? size : 0; - IsLowerTriangular ? k2>0 : k2<size; - IsLowerTriangular ? k2-=kc : k2+=kc) + for(int k2=IsLower ? size : 0; + IsLower ? k2>0 : k2<size; + IsLower ? k2-=kc : k2+=kc) { - const int actual_kc = std::min(IsLowerTriangular ? k2 : size-k2, kc); - int actual_k2 = IsLowerTriangular ? k2-actual_kc : k2; + const int actual_kc = std::min(IsLower ? k2 : size-k2, kc); + int actual_k2 = IsLower ? k2-actual_kc : k2; pack_rhs(blockB, &rhs(actual_k2,0), rhsStride, alpha, actual_kc, cols); @@ -149,7 +149,7 @@ struct ei_product_triangular_matrix_matrix<Scalar,Mode,true, for (int k1=0; k1<actual_kc; k1+=SmallPanelWidth) { int actualPanelWidth = std::min<int>(actual_kc-k1, SmallPanelWidth); - int lengthTarget = IsLowerTriangular ? actual_kc-k1-actualPanelWidth : k1; + int lengthTarget = IsLower ? actual_kc-k1-actualPanelWidth : k1; int startBlock = actual_k2+k1; int blockBOffset = k1; @@ -158,9 +158,9 @@ struct ei_product_triangular_matrix_matrix<Scalar,Mode,true, // To this end we do an extra triangular copy to small temporary buffer for (int k=0;k<actualPanelWidth;++k) { - if (!(Mode&UnitDiagBit)) + if (!(Mode&UnitDiag)) triangularBuffer.coeffRef(k,k) = lhs(startBlock+k,startBlock+k); - for (int i=IsLowerTriangular ? k+1 : 0; IsLowerTriangular ? i<actualPanelWidth : i<k; ++i) + for (int i=IsLower ? k+1 : 0; IsLower ? i<actualPanelWidth : i<k; ++i) triangularBuffer.coeffRef(i,k) = lhs(startBlock+i,startBlock+k); } pack_lhs(blockA, triangularBuffer.data(), triangularBuffer.stride(), actualPanelWidth, actualPanelWidth); @@ -171,7 +171,7 @@ struct ei_product_triangular_matrix_matrix<Scalar,Mode,true, // GEBP with remaining micro panel if (lengthTarget>0) { - int startTarget = IsLowerTriangular ? actual_k2+k1+actualPanelWidth : actual_k2; + int startTarget = IsLower ? actual_k2+k1+actualPanelWidth : actual_k2; pack_lhs(blockA, &lhs(startTarget,startBlock), lhsStride, actualPanelWidth, lengthTarget); @@ -182,8 +182,8 @@ struct ei_product_triangular_matrix_matrix<Scalar,Mode,true, } // the part below the diagonal => GEPP { - int start = IsLowerTriangular ? k2 : 0; - int end = IsLowerTriangular ? size : actual_k2; + int start = IsLower ? k2 : 0; + int end = IsLower ? size : actual_k2; for(int i2=start; i2<end; i2+=mc) { const int actual_mc = std::min(i2+mc,end)-i2; @@ -227,7 +227,7 @@ struct ei_product_triangular_matrix_matrix<Scalar,Mode,false, typedef ei_product_blocking_traits<Scalar> Blocking; enum { SmallPanelWidth = EIGEN_ENUM_MAX(Blocking::mr,Blocking::nr), - IsLowerTriangular = (Mode&LowerTriangular) == LowerTriangular + IsLower = (Mode&Lower) == Lower }; int kc = std::min<int>(Blocking::Max_kc/4,size); // cache block size along the K direction @@ -245,16 +245,16 @@ struct ei_product_triangular_matrix_matrix<Scalar,Mode,false, ei_gemm_pack_rhs<Scalar,Blocking::nr,RhsStorageOrder> pack_rhs; ei_gemm_pack_rhs<Scalar,Blocking::nr,RhsStorageOrder,true> pack_rhs_panel; - for(int k2=IsLowerTriangular ? 0 : size; - IsLowerTriangular ? k2<size : k2>0; - IsLowerTriangular ? k2+=kc : k2-=kc) + for(int k2=IsLower ? 0 : size; + IsLower ? k2<size : k2>0; + IsLower ? k2+=kc : k2-=kc) { - const int actual_kc = std::min(IsLowerTriangular ? size-k2 : k2, kc); - int actual_k2 = IsLowerTriangular ? k2 : k2-actual_kc; - int rs = IsLowerTriangular ? actual_k2 : size - k2; + const int actual_kc = std::min(IsLower ? size-k2 : k2, kc); + int actual_k2 = IsLower ? k2 : k2-actual_kc; + int rs = IsLower ? actual_k2 : size - k2; Scalar* geb = blockB+actual_kc*actual_kc*Blocking::PacketSize; - pack_rhs(geb, &rhs(actual_k2,IsLowerTriangular ? 0 : k2), rhsStride, alpha, actual_kc, rs); + pack_rhs(geb, &rhs(actual_k2,IsLower ? 0 : k2), rhsStride, alpha, actual_kc, rs); // pack the triangular part of the rhs padding the unrolled blocks with zeros { @@ -262,8 +262,8 @@ struct ei_product_triangular_matrix_matrix<Scalar,Mode,false, { int actualPanelWidth = std::min<int>(actual_kc-j2, SmallPanelWidth); int actual_j2 = actual_k2 + j2; - int panelOffset = IsLowerTriangular ? j2+actualPanelWidth : 0; - int panelLength = IsLowerTriangular ? actual_kc-j2-actualPanelWidth : j2; + int panelOffset = IsLower ? j2+actualPanelWidth : 0; + int panelLength = IsLower ? actual_kc-j2-actualPanelWidth : j2; // general part pack_rhs_panel(blockB+j2*actual_kc*Blocking::PacketSize, &rhs(actual_k2+panelOffset, actual_j2), rhsStride, alpha, @@ -273,9 +273,9 @@ struct ei_product_triangular_matrix_matrix<Scalar,Mode,false, // append the triangular part via a temporary buffer for (int j=0;j<actualPanelWidth;++j) { - if (!(Mode&UnitDiagBit)) + if (!(Mode&UnitDiag)) triangularBuffer.coeffRef(j,j) = rhs(actual_j2+j,actual_j2+j); - for (int k=IsLowerTriangular ? j+1 : 0; IsLowerTriangular ? k<actualPanelWidth : k<j; ++k) + for (int k=IsLower ? j+1 : 0; IsLower ? k<actualPanelWidth : k<j; ++k) triangularBuffer.coeffRef(k,j) = rhs(actual_j2+k,actual_j2+j); } @@ -296,8 +296,8 @@ struct ei_product_triangular_matrix_matrix<Scalar,Mode,false, for (int j2=0; j2<actual_kc; j2+=SmallPanelWidth) { int actualPanelWidth = std::min<int>(actual_kc-j2, SmallPanelWidth); - int panelLength = IsLowerTriangular ? actual_kc-j2 : j2+actualPanelWidth; - int blockOffset = IsLowerTriangular ? j2 : 0; + int panelLength = IsLower ? actual_kc-j2 : j2+actualPanelWidth; + int blockOffset = IsLower ? j2 : 0; gebp_kernel(res+i2+(actual_k2+j2)*resStride, resStride, blockA, blockB+j2*actual_kc*Blocking::PacketSize, @@ -306,7 +306,7 @@ struct ei_product_triangular_matrix_matrix<Scalar,Mode,false, blockOffset, blockOffset*Blocking::PacketSize);// offsets } } - gebp_kernel(res+i2+(IsLowerTriangular ? 0 : k2)*resStride, resStride, + gebp_kernel(res+i2+(IsLower ? 0 : k2)*resStride, resStride, blockA, geb, actual_mc, actual_kc, rs); } } diff --git a/Eigen/src/Core/products/TriangularMatrixVector.h b/Eigen/src/Core/products/TriangularMatrixVector.h index fc3188cd8..2cad48eb9 100644 --- a/Eigen/src/Core/products/TriangularMatrixVector.h +++ b/Eigen/src/Core/products/TriangularMatrixVector.h @@ -34,8 +34,8 @@ struct ei_product_triangular_vector_selector<Lhs,Rhs,Result,Mode,ConjLhs,ConjRhs { typedef typename Rhs::Scalar Scalar; enum { - IsLowerTriangular = ((Mode&LowerTriangularBit)==LowerTriangularBit), - HasUnitDiag = (Mode & UnitDiagBit)==UnitDiagBit + IsLower = ((Mode&Lower)==Lower), + HasUnitDiag = (Mode & UnitDiag)==UnitDiag }; static EIGEN_DONT_INLINE void run(const Lhs& lhs, const Rhs& rhs, Result& res, typename ei_traits<Lhs>::Scalar alpha) { @@ -50,17 +50,17 @@ struct ei_product_triangular_vector_selector<Lhs,Rhs,Result,Mode,ConjLhs,ConjRhs for (int k=0; k<actualPanelWidth; ++k) { int i = pi + k; - int s = IsLowerTriangular ? (HasUnitDiag ? i+1 : i ) : pi; - int r = IsLowerTriangular ? actualPanelWidth-k : k+1; + int s = IsLower ? (HasUnitDiag ? i+1 : i ) : pi; + int r = IsLower ? actualPanelWidth-k : k+1; if ((!HasUnitDiag) || (--r)>0) res.segment(s,r) += (alpha * cjRhs.coeff(i)) * cjLhs.col(i).segment(s,r); if (HasUnitDiag) res.coeffRef(i) += alpha * cjRhs.coeff(i); } - int r = IsLowerTriangular ? size - pi - actualPanelWidth : pi; + int r = IsLower ? size - pi - actualPanelWidth : pi; if (r>0) { - int s = IsLowerTriangular ? pi+actualPanelWidth : 0; + int s = IsLower ? pi+actualPanelWidth : 0; ei_cache_friendly_product_colmajor_times_vector<ConjLhs,ConjRhs>( r, &(lhs.const_cast_derived().coeffRef(s,pi)), lhs.stride(), @@ -77,8 +77,8 @@ struct ei_product_triangular_vector_selector<Lhs,Rhs,Result,Mode,ConjLhs,ConjRhs { typedef typename Rhs::Scalar Scalar; enum { - IsLowerTriangular = ((Mode&LowerTriangularBit)==LowerTriangularBit), - HasUnitDiag = (Mode & UnitDiagBit)==UnitDiagBit + IsLower = ((Mode&Lower)==Lower), + HasUnitDiag = (Mode & UnitDiag)==UnitDiag }; static void run(const Lhs& lhs, const Rhs& rhs, Result& res, typename ei_traits<Lhs>::Scalar alpha) { @@ -92,17 +92,17 @@ struct ei_product_triangular_vector_selector<Lhs,Rhs,Result,Mode,ConjLhs,ConjRhs for (int k=0; k<actualPanelWidth; ++k) { int i = pi + k; - int s = IsLowerTriangular ? pi : (HasUnitDiag ? i+1 : i); - int r = IsLowerTriangular ? k+1 : actualPanelWidth-k; + int s = IsLower ? pi : (HasUnitDiag ? i+1 : i); + int r = IsLower ? k+1 : actualPanelWidth-k; if ((!HasUnitDiag) || (--r)>0) res.coeffRef(i) += alpha * (cjLhs.row(i).segment(s,r).cwiseProduct(cjRhs.segment(s,r).transpose())).sum(); if (HasUnitDiag) res.coeffRef(i) += alpha * cjRhs.coeff(i); } - int r = IsLowerTriangular ? pi : size - pi - actualPanelWidth; + int r = IsLower ? pi : size - pi - actualPanelWidth; if (r>0) { - int s = IsLowerTriangular ? 0 : pi + actualPanelWidth; + int s = IsLower ? 0 : pi + actualPanelWidth; Block<Result,Dynamic,1> target(res,pi,0,actualPanelWidth,1); ei_cache_friendly_product_rowmajor_times_vector<ConjLhs,ConjRhs>( &(lhs.const_cast_derived().coeffRef(pi,s)), lhs.stride(), diff --git a/Eigen/src/Core/products/TriangularSolverMatrix.h b/Eigen/src/Core/products/TriangularSolverMatrix.h index e49fac956..23a645d7c 100644 --- a/Eigen/src/Core/products/TriangularSolverMatrix.h +++ b/Eigen/src/Core/products/TriangularSolverMatrix.h @@ -36,7 +36,7 @@ struct ei_triangular_solve_matrix<Scalar,Side,Mode,Conjugate,TriStorageOrder,Row { ei_triangular_solve_matrix< Scalar, Side==OnTheLeft?OnTheRight:OnTheLeft, - (Mode&UnitDiagBit) | ((Mode&UpperTriangular) ? LowerTriangular : UpperTriangular), + (Mode&UnitDiag) | ((Mode&Upper) ? Lower : Upper), NumTraits<Scalar>::IsComplex && Conjugate, TriStorageOrder==RowMajor ? ColMajor : RowMajor, ColMajor> ::run(size, cols, tri, triStride, _other, otherStride); @@ -60,7 +60,7 @@ struct ei_triangular_solve_matrix<Scalar,OnTheLeft,Mode,Conjugate,TriStorageOrde typedef ei_product_blocking_traits<Scalar> Blocking; enum { SmallPanelWidth = EIGEN_ENUM_MAX(Blocking::mr,Blocking::nr), - IsLowerTriangular = (Mode&LowerTriangular) == LowerTriangular + IsLower = (Mode&Lower) == Lower }; int kc = std::min<int>(Blocking::Max_kc/4,size); // cache block size along the K direction @@ -73,11 +73,11 @@ struct ei_triangular_solve_matrix<Scalar,OnTheLeft,Mode,Conjugate,TriStorageOrde ei_gebp_kernel<Scalar, Blocking::mr, Blocking::nr, ei_conj_helper<Conjugate,false> > gebp_kernel; ei_gemm_pack_lhs<Scalar,Blocking::mr,TriStorageOrder> pack_lhs; - for(int k2=IsLowerTriangular ? 0 : size; - IsLowerTriangular ? k2<size : k2>0; - IsLowerTriangular ? k2+=kc : k2-=kc) + for(int k2=IsLower ? 0 : size; + IsLower ? k2<size : k2>0; + IsLower ? k2+=kc : k2-=kc) { - const int actual_kc = std::min(IsLowerTriangular ? size-k2 : k2, kc); + const int actual_kc = std::min(IsLower ? size-k2 : k2, kc); // We have selected and packed a big horizontal panel R1 of rhs. Let B be the packed copy of this panel, // and R2 the remaining part of rhs. The corresponding vertical panel of lhs is split into @@ -101,11 +101,11 @@ struct ei_triangular_solve_matrix<Scalar,OnTheLeft,Mode,Conjugate,TriStorageOrde for (int k=0; k<actualPanelWidth; ++k) { // TODO write a small kernel handling this (can be shared with trsv) - int i = IsLowerTriangular ? k2+k1+k : k2-k1-k-1; - int s = IsLowerTriangular ? k2+k1 : i+1; + int i = IsLower ? k2+k1+k : k2-k1-k-1; + int s = IsLower ? k2+k1 : i+1; int rs = actualPanelWidth - k - 1; // remaining size - Scalar a = (Mode & UnitDiagBit) ? Scalar(1) : Scalar(1)/conj(tri(i,i)); + Scalar a = (Mode & UnitDiag) ? Scalar(1) : Scalar(1)/conj(tri(i,i)); for (int j=0; j<cols; ++j) { if (TriStorageOrder==RowMajor) @@ -120,7 +120,7 @@ struct ei_triangular_solve_matrix<Scalar,OnTheLeft,Mode,Conjugate,TriStorageOrde } else { - int s = IsLowerTriangular ? i+1 : i-rs; + int s = IsLower ? i+1 : i-rs; Scalar b = (other(i,j) *= a); Scalar* r = &other(s,j); const Scalar* l = &tri(s,i); @@ -131,8 +131,8 @@ struct ei_triangular_solve_matrix<Scalar,OnTheLeft,Mode,Conjugate,TriStorageOrde } int lengthTarget = actual_kc-k1-actualPanelWidth; - int startBlock = IsLowerTriangular ? k2+k1 : k2-k1-actualPanelWidth; - int blockBOffset = IsLowerTriangular ? k1 : lengthTarget; + int startBlock = IsLower ? k2+k1 : k2-k1-actualPanelWidth; + int blockBOffset = IsLower ? k1 : lengthTarget; // update the respective rows of B from other ei_gemm_pack_rhs<Scalar, Blocking::nr, ColMajor, true>() @@ -141,7 +141,7 @@ struct ei_triangular_solve_matrix<Scalar,OnTheLeft,Mode,Conjugate,TriStorageOrde // GEBP if (lengthTarget>0) { - int startTarget = IsLowerTriangular ? k2+k1+actualPanelWidth : k2-actual_kc; + int startTarget = IsLower ? k2+k1+actualPanelWidth : k2-actual_kc; pack_lhs(blockA, &tri(startTarget,startBlock), triStride, actualPanelWidth, lengthTarget); @@ -153,14 +153,14 @@ struct ei_triangular_solve_matrix<Scalar,OnTheLeft,Mode,Conjugate,TriStorageOrde // R2 = A2 * B => GEPP { - int start = IsLowerTriangular ? k2+kc : 0; - int end = IsLowerTriangular ? size : k2-kc; + int start = IsLower ? k2+kc : 0; + int end = IsLower ? size : k2-kc; for(int i2=start; i2<end; i2+=mc) { const int actual_mc = std::min(mc,end-i2); if (actual_mc>0) { - pack_lhs(blockA, &tri(i2, IsLowerTriangular ? k2 : k2-kc), triStride, actual_kc, actual_mc); + pack_lhs(blockA, &tri(i2, IsLower ? k2 : k2-kc), triStride, actual_kc, actual_mc); gebp_kernel(_other+i2, otherStride, blockA, blockB, actual_mc, actual_kc, cols); } @@ -191,7 +191,7 @@ struct ei_triangular_solve_matrix<Scalar,OnTheRight,Mode,Conjugate,TriStorageOrd enum { RhsStorageOrder = TriStorageOrder, SmallPanelWidth = EIGEN_ENUM_MAX(Blocking::mr,Blocking::nr), - IsLowerTriangular = (Mode&LowerTriangular) == LowerTriangular + IsLower = (Mode&Lower) == Lower }; int kc = std::min<int>(Blocking::Max_kc/4,size); // cache block size along the K direction @@ -206,15 +206,15 @@ struct ei_triangular_solve_matrix<Scalar,OnTheRight,Mode,Conjugate,TriStorageOrd ei_gemm_pack_rhs<Scalar,Blocking::nr,RhsStorageOrder,true> pack_rhs_panel; ei_gemm_pack_lhs<Scalar, Blocking::mr, ColMajor, false, true> pack_lhs_panel; - for(int k2=IsLowerTriangular ? size : 0; - IsLowerTriangular ? k2>0 : k2<size; - IsLowerTriangular ? k2-=kc : k2+=kc) + for(int k2=IsLower ? size : 0; + IsLower ? k2>0 : k2<size; + IsLower ? k2-=kc : k2+=kc) { - const int actual_kc = std::min(IsLowerTriangular ? k2 : size-k2, kc); - int actual_k2 = IsLowerTriangular ? k2-actual_kc : k2 ; + const int actual_kc = std::min(IsLower ? k2 : size-k2, kc); + int actual_k2 = IsLower ? k2-actual_kc : k2 ; - int startPanel = IsLowerTriangular ? 0 : k2+actual_kc; - int rs = IsLowerTriangular ? actual_k2 : size - actual_k2 - actual_kc; + int startPanel = IsLower ? 0 : k2+actual_kc; + int rs = IsLower ? actual_k2 : size - actual_k2 - actual_kc; Scalar* geb = blockB+actual_kc*actual_kc*Blocking::PacketSize; if (rs>0) pack_rhs(geb, &rhs(actual_k2,startPanel), triStride, -1, actual_kc, rs); @@ -226,8 +226,8 @@ struct ei_triangular_solve_matrix<Scalar,OnTheRight,Mode,Conjugate,TriStorageOrd { int actualPanelWidth = std::min<int>(actual_kc-j2, SmallPanelWidth); int actual_j2 = actual_k2 + j2; - int panelOffset = IsLowerTriangular ? j2+actualPanelWidth : 0; - int panelLength = IsLowerTriangular ? actual_kc-j2-actualPanelWidth : j2; + int panelOffset = IsLower ? j2+actualPanelWidth : 0; + int panelLength = IsLower ? actual_kc-j2-actualPanelWidth : j2; if (panelLength>0) pack_rhs_panel(blockB+j2*actual_kc*Blocking::PacketSize, @@ -244,17 +244,17 @@ struct ei_triangular_solve_matrix<Scalar,OnTheRight,Mode,Conjugate,TriStorageOrd // triangular solver kernel { // for each small block of the diagonal (=> vertical panels of rhs) - for (int j2 = IsLowerTriangular + for (int j2 = IsLower ? (actual_kc - ((actual_kc%SmallPanelWidth) ? (actual_kc%SmallPanelWidth) : SmallPanelWidth)) : 0; - IsLowerTriangular ? j2>=0 : j2<actual_kc; - IsLowerTriangular ? j2-=SmallPanelWidth : j2+=SmallPanelWidth) + IsLower ? j2>=0 : j2<actual_kc; + IsLower ? j2-=SmallPanelWidth : j2+=SmallPanelWidth) { int actualPanelWidth = std::min<int>(actual_kc-j2, SmallPanelWidth); int absolute_j2 = actual_k2 + j2; - int panelOffset = IsLowerTriangular ? j2+actualPanelWidth : 0; - int panelLength = IsLowerTriangular ? actual_kc - j2 - actualPanelWidth : j2; + int panelOffset = IsLower ? j2+actualPanelWidth : 0; + int panelLength = IsLower ? actual_kc - j2 - actualPanelWidth : j2; // GEBP if(panelLength>0) @@ -269,17 +269,17 @@ struct ei_triangular_solve_matrix<Scalar,OnTheRight,Mode,Conjugate,TriStorageOrd // unblocked triangular solve for (int k=0; k<actualPanelWidth; ++k) { - int j = IsLowerTriangular ? absolute_j2+actualPanelWidth-k-1 : absolute_j2+k; + int j = IsLower ? absolute_j2+actualPanelWidth-k-1 : absolute_j2+k; Scalar* r = &lhs(i2,j); for (int k3=0; k3<k; ++k3) { - Scalar b = conj(rhs(IsLowerTriangular ? j+1+k3 : absolute_j2+k3,j)); - Scalar* a = &lhs(i2,IsLowerTriangular ? j+1+k3 : absolute_j2+k3); + Scalar b = conj(rhs(IsLower ? j+1+k3 : absolute_j2+k3,j)); + Scalar* a = &lhs(i2,IsLower ? j+1+k3 : absolute_j2+k3); for (int i=0; i<actual_mc; ++i) r[i] -= a[i] * b; } - Scalar b = (Mode & UnitDiagBit) ? Scalar(1) : Scalar(1)/conj(rhs(j,j)); + Scalar b = (Mode & UnitDiag) ? Scalar(1) : Scalar(1)/conj(rhs(j,j)); for (int i=0; i<actual_mc; ++i) r[i] *= b; } diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h index 8483772df..01c1aeb2f 100644 --- a/Eigen/src/Core/util/Constants.h +++ b/Eigen/src/Core/util/Constants.h @@ -149,44 +149,18 @@ const unsigned int DirectAccessBit = 0x20; * means the first coefficient packet is guaranteed to be aligned */ const unsigned int AlignedBit = 0x40; -/** \ingroup flags - * - * means all diagonal coefficients are equal to 0 */ -const unsigned int ZeroDiagBit = 0x80; - -/** \ingroup flags - * - * means all diagonal coefficients are equal to 1 */ -const unsigned int UnitDiagBit = 0x100; - -/** \ingroup flags - * - * means the matrix is selfadjoint (M=M*). */ -const unsigned int SelfAdjointBit = 0x200; - -/** \ingroup flags - * - * means the strictly lower triangular part is 0 */ -const unsigned int UpperTriangularBit = 0x400; - -/** \ingroup flags - * - * means the strictly upper triangular part is 0 */ -const unsigned int LowerTriangularBit = 0x800; // list of flags that are inherited by default const unsigned int HereditaryBits = RowMajorBit | EvalBeforeNestingBit | EvalBeforeAssigningBit; -// Possible values for the Mode parameter of part() -const unsigned int UpperTriangular = UpperTriangularBit; -const unsigned int StrictlyUpperTriangular = UpperTriangularBit | ZeroDiagBit; -const unsigned int LowerTriangular = LowerTriangularBit; -const unsigned int StrictlyLowerTriangular = LowerTriangularBit | ZeroDiagBit; -const unsigned int SelfAdjoint = SelfAdjointBit; -const unsigned int UnitUpperTriangular = UpperTriangularBit | UnitDiagBit; -const unsigned int UnitLowerTriangular = LowerTriangularBit | UnitDiagBit; +// Possible values for the Mode parameter of triangularView() +enum { + Lower=0x1, Upper=0x2, UnitDiag=0x4, ZeroDiag=0x8, + UnitLower=UnitDiag|Lower, UnitUpper=UnitDiag|Upper, + StrictlyLower=ZeroDiag|Lower, StrictlyUpper=ZeroDiag|Upper, + SelfAdjoint=0x10}; enum { Unaligned=0, Aligned=1 }; enum { ConditionalJumpCost = 5 }; diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index 09c6c9528..6959a1c1b 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -132,7 +132,7 @@ template<typename MatrixType> class ColPivHouseholderQR; template<typename MatrixType> class FullPivHouseholderQR; template<typename MatrixType> class SVD; template<typename MatrixType, unsigned int Options = 0> class JacobiSVD; -template<typename MatrixType, int UpLo = LowerTriangular> class LLT; +template<typename MatrixType, int UpLo = Lower> class LLT; template<typename MatrixType> class LDLT; template<typename VectorsType, typename CoeffsType> class HouseholderSequence; template<typename Scalar> class PlanarRotation; diff --git a/Eigen/src/Core/util/XprHelper.h b/Eigen/src/Core/util/XprHelper.h index 4bff09252..9b17a2b0e 100644 --- a/Eigen/src/Core/util/XprHelper.h +++ b/Eigen/src/Core/util/XprHelper.h @@ -263,8 +263,7 @@ template<typename T, int n=1, typename PlainMatrixType = typename ei_eval<T>::ty template<unsigned int Flags> struct ei_are_flags_consistent { - enum { ret = !( (Flags&UnitDiagBit && Flags&ZeroDiagBit) ) - }; + enum { ret = true }; }; /** \internal Helper base class to add a scalar multiple operator diff --git a/Eigen/src/Eigen2Support/Flagged.h b/Eigen/src/Eigen2Support/Flagged.h index 94e7f9119..bed110b64 100644 --- a/Eigen/src/Eigen2Support/Flagged.h +++ b/Eigen/src/Eigen2Support/Flagged.h @@ -109,6 +109,12 @@ 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; + + template<typename OtherDerived> + void solveTriangularInPlace(const MatrixBase<OtherDerived>& other) const; + protected: ExpressionTypeNested m_matrix; }; diff --git a/Eigen/src/Eigen2Support/TriangularSolver.h b/Eigen/src/Eigen2Support/TriangularSolver.h index e69de29bb..94b92577e 100644 --- a/Eigen/src/Eigen2Support/TriangularSolver.h +++ b/Eigen/src/Eigen2Support/TriangularSolver.h @@ -0,0 +1,53 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2010 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_SOLVER2_H +#define EIGEN_TRIANGULAR_SOLVER2_H + +const unsigned int UnitDiagBit = UnitDiag; +const unsigned int SelfAdjointBit = SelfAdjoint; +const unsigned int UpperTriangularBit = Upper; +const unsigned int LowerTriangularBit = Lower; + +const unsigned int UpperTriangular = Upper; +const unsigned int LowerTriangular = Lower; +const unsigned int UnitUpperTriangular = UnitUpper; +const unsigned int UnitLowerTriangular = UnitLower; + +template<typename ExpressionType, unsigned int Added, unsigned int Removed> +template<typename OtherDerived> +typename ExpressionType::PlainMatrixType +Flagged<ExpressionType,Added,Removed>::solveTriangular(const MatrixBase<OtherDerived>& other) const +{ + return m_matrix.template triangularView<Added>.solve(other.derived()); +} + +template<typename ExpressionType, unsigned int Added, unsigned int Removed> +template<typename OtherDerived> +void Flagged<ExpressionType,Added,Removed>::solveTriangularInPlace(const MatrixBase<OtherDerived>& other) const +{ + m_matrix.template triangularView<Added>.solveInPlace(other.derived()); +} + +#endif // EIGEN_TRIANGULAR_SOLVER2_H diff --git a/Eigen/src/Eigenvalues/HessenbergDecomposition.h b/Eigen/src/Eigenvalues/HessenbergDecomposition.h index 636b2f4f7..e3f64b807 100644 --- a/Eigen/src/Eigenvalues/HessenbergDecomposition.h +++ b/Eigen/src/Eigenvalues/HessenbergDecomposition.h @@ -199,7 +199,7 @@ HessenbergDecomposition<MatrixType>::matrixH() const int n = m_matrix.rows(); MatrixType matH = m_matrix; if (n>2) - matH.corner(BottomLeft,n-2, n-2).template triangularView<LowerTriangular>().setZero(); + matH.corner(BottomLeft,n-2, n-2).template triangularView<Lower>().setZero(); return matH; } diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h index 336976517..960e9b417 100644 --- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h +++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h @@ -247,11 +247,11 @@ compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors matC.adjointInPlace(); // this version works too: // matC = matC.transpose(); -// cholB.matrixL().conjugate().template marked<LowerTriangular>().solveTriangularInPlace(matC); +// cholB.matrixL().conjugate().template marked<Lower>().solveTriangularInPlace(matC); // matC = matC.transpose(); // FIXME: this should work: (currently it only does for small matrices) // Transpose<MatrixType> trMatC(matC); -// cholB.matrixL().conjugate().eval().template marked<LowerTriangular>().solveTriangularInPlace(trMatC); +// cholB.matrixL().conjugate().eval().template marked<Lower>().solveTriangularInPlace(trMatC); compute(matC, computeEigenvectors); @@ -275,7 +275,7 @@ template<typename Derived> inline Matrix<typename NumTraits<typename ei_traits<Derived>::Scalar>::Real, ei_traits<Derived>::ColsAtCompileTime, 1> MatrixBase<Derived>::eigenvalues() const { - ei_assert(Flags&SelfAdjointBit); + ei_assert(Flags&SelfAdjoint); return SelfAdjointEigenSolver<typename Derived::PlainMatrixType>(eval(),false).eigenvalues(); } @@ -316,7 +316,7 @@ template<typename Derived> inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real MatrixBase<Derived>::operatorNorm() const { - return ei_operatorNorm_selector<Derived, Flags&SelfAdjointBit> + return ei_operatorNorm_selector<Derived, Flags&SelfAdjoint> ::operatorNorm(derived()); } diff --git a/Eigen/src/Eigenvalues/Tridiagonalization.h b/Eigen/src/Eigenvalues/Tridiagonalization.h index e43605b0f..ffd374eb2 100644 --- a/Eigen/src/Eigenvalues/Tridiagonalization.h +++ b/Eigen/src/Eigenvalues/Tridiagonalization.h @@ -172,8 +172,8 @@ Tridiagonalization<MatrixType>::matrixT(void) const matT.corner(TopRight,n-1, n-1).diagonal() = subDiagonal().template cast<Scalar>().conjugate(); if (n>2) { - matT.corner(TopRight,n-2, n-2).template triangularView<UpperTriangular>().setZero(); - matT.corner(BottomLeft,n-2, n-2).template triangularView<LowerTriangular>().setZero(); + matT.corner(TopRight,n-2, n-2).template triangularView<Upper>().setZero(); + matT.corner(BottomLeft,n-2, n-2).template triangularView<Lower>().setZero(); } return matT; } @@ -208,12 +208,12 @@ void Tridiagonalization<MatrixType>::_compute(MatrixType& matA, CoeffVectorType& // i.e., A = H A H' where H = I - h v v' and v = matA.col(i).tail(n-i-1) matA.col(i).coeffRef(i+1) = 1; - hCoeffs.tail(n-i-1) = (matA.corner(BottomRight,remainingSize,remainingSize).template selfadjointView<LowerTriangular>() + hCoeffs.tail(n-i-1) = (matA.corner(BottomRight,remainingSize,remainingSize).template selfadjointView<Lower>() * (ei_conj(h) * matA.col(i).tail(remainingSize))); hCoeffs.tail(n-i-1) += (ei_conj(h)*Scalar(-0.5)*(hCoeffs.tail(remainingSize).dot(matA.col(i).tail(remainingSize)))) * matA.col(i).tail(n-i-1); - matA.corner(BottomRight, remainingSize, remainingSize).template selfadjointView<LowerTriangular>() + matA.corner(BottomRight, remainingSize, remainingSize).template selfadjointView<Lower>() .rankUpdate(matA.col(i).tail(remainingSize), hCoeffs.tail(remainingSize), -1); matA.col(i).coeffRef(i+1) = beta; diff --git a/Eigen/src/LU/Determinant.h b/Eigen/src/LU/Determinant.h index 27ad6abe9..bae01256e 100644 --- a/Eigen/src/LU/Determinant.h +++ b/Eigen/src/LU/Determinant.h @@ -41,14 +41,8 @@ const typename Derived::Scalar ei_bruteforce_det4_helper * (matrix.coeff(m,2) * matrix.coeff(n,3) - matrix.coeff(n,2) * matrix.coeff(m,3)); } -// FIXME update computation of triangular det - -const int TriangularDeterminant = 0; - template<typename Derived, - int DeterminantType = - (Derived::Flags & (UpperTriangularBit | LowerTriangularBit)) - ? TriangularDeterminant : Derived::RowsAtCompileTime + int DeterminantType = Derived::RowsAtCompileTime > struct ei_determinant_impl { static inline typename ei_traits<Derived>::Scalar run(const Derived& m) @@ -57,19 +51,6 @@ template<typename Derived, } }; -template<typename Derived> struct ei_determinant_impl<Derived, TriangularDeterminant> -{ - static inline typename ei_traits<Derived>::Scalar run(const Derived& m) - { - if (Derived::Flags & UnitDiagBit) - return 1; - else if (Derived::Flags & ZeroDiagBit) - return 0; - else - return m.diagonal().redux(ei_scalar_product_op<typename ei_traits<Derived>::Scalar>()); - } -}; - template<typename Derived> struct ei_determinant_impl<Derived, 1> { static inline typename ei_traits<Derived>::Scalar run(const Derived& m) diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h index f62dcc1db..148ddcd23 100644 --- a/Eigen/src/LU/FullPivLU.h +++ b/Eigen/src/LU/FullPivLU.h @@ -541,7 +541,7 @@ struct ei_kernel_retval<FullPivLU<_MatrixType> > m.row(i).tail(cols-i) = dec().matrixLU().row(pivots.coeff(i)).tail(cols-i); } m.block(0, 0, rank(), rank()); - m.block(0, 0, rank(), rank()).template triangularView<StrictlyLowerTriangular>().setZero(); + m.block(0, 0, rank(), rank()).template triangularView<StrictlyLower>().setZero(); for(int i = 0; i < rank(); ++i) m.col(i).swap(m.col(pivots.coeff(i))); @@ -549,7 +549,7 @@ struct ei_kernel_retval<FullPivLU<_MatrixType> > // notice that the math behind this suggests that we should apply this to the // negative of the RHS, but for performance we just put the negative sign elsewhere, see below. m.corner(TopLeft, rank(), rank()) - .template triangularView<UpperTriangular>().solveInPlace( + .template triangularView<Upper>().solveInPlace( m.corner(TopRight, rank(), dimker) ); @@ -638,7 +638,7 @@ struct ei_solve_retval<FullPivLU<_MatrixType>, Rhs> // Step 2 dec().matrixLU() .corner(Eigen::TopLeft,smalldim,smalldim) - .template triangularView<UnitLowerTriangular>() + .template triangularView<UnitLower>() .solveInPlace(c.corner(Eigen::TopLeft, smalldim, c.cols())); if(rows>cols) { @@ -650,7 +650,7 @@ struct ei_solve_retval<FullPivLU<_MatrixType>, Rhs> // Step 3 dec().matrixLU() .corner(TopLeft, nonzero_pivots, nonzero_pivots) - .template triangularView<UpperTriangular>() + .template triangularView<Upper>() .solveInPlace(c.corner(TopLeft, nonzero_pivots, c.cols())); // Step 4 diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h index 6bb5c3ac7..f50aa4535 100644 --- a/Eigen/src/LU/PartialPivLU.h +++ b/Eigen/src/LU/PartialPivLU.h @@ -349,7 +349,7 @@ struct ei_partial_lu_impl A_2.row(i).swap(A_2.row(row_transpositions[i])); // A12 = A11^-1 A12 - A11.template triangularView<UnitLowerTriangular>().solveInPlace(A12); + A11.template triangularView<UnitLower>().solveInPlace(A12); A22.noalias() -= A21 * A12; } @@ -426,10 +426,10 @@ struct ei_solve_retval<PartialPivLU<_MatrixType>, Rhs> dst = dec().permutationP() * rhs(); // Step 2 - dec().matrixLU().template triangularView<UnitLowerTriangular>().solveInPlace(dst); + dec().matrixLU().template triangularView<UnitLower>().solveInPlace(dst); // Step 3 - dec().matrixLU().template triangularView<UpperTriangular>().solveInPlace(dst); + dec().matrixLU().template triangularView<Upper>().solveInPlace(dst); } }; diff --git a/Eigen/src/QR/ColPivHouseholderQR.h b/Eigen/src/QR/ColPivHouseholderQR.h index 8b705de3c..f4edc9324 100644 --- a/Eigen/src/QR/ColPivHouseholderQR.h +++ b/Eigen/src/QR/ColPivHouseholderQR.h @@ -381,7 +381,7 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const m_nonzero_pivots = k; m_hCoeffs.tail(size-k).setZero(); m_qr.corner(BottomRight,rows-k,cols-k) - .template triangularView<StrictlyLowerTriangular>() + .template triangularView<StrictlyLower>() .setZero(); break; } @@ -453,7 +453,7 @@ struct ei_solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs> dec().matrixQR() .corner(TopLeft, nonzero_pivots, nonzero_pivots) - .template triangularView<UpperTriangular>() + .template triangularView<Upper>() .solveInPlace(c.corner(TopLeft, nonzero_pivots, c.cols())); @@ -461,7 +461,7 @@ struct ei_solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs> d.corner(TopLeft, nonzero_pivots, c.cols()) = dec().matrixQR() .corner(TopLeft, nonzero_pivots, nonzero_pivots) - .template triangularView<UpperTriangular>() + .template triangularView<Upper>() * c.corner(TopLeft, nonzero_pivots, c.cols()); for(int i = 0; i < nonzero_pivots; ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i); diff --git a/Eigen/src/QR/FullPivHouseholderQR.h b/Eigen/src/QR/FullPivHouseholderQR.h index 598f44dc3..717ff19f8 100644 --- a/Eigen/src/QR/FullPivHouseholderQR.h +++ b/Eigen/src/QR/FullPivHouseholderQR.h @@ -376,7 +376,7 @@ struct ei_solve_retval<FullPivHouseholderQR<_MatrixType>, Rhs> } dec().matrixQR() .corner(TopLeft, dec().rank(), dec().rank()) - .template triangularView<UpperTriangular>() + .template triangularView<Upper>() .solveInPlace(c.corner(TopLeft, dec().rank(), c.cols())); for(int i = 0; i < dec().rank(); ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i); diff --git a/Eigen/src/QR/HouseholderQR.h b/Eigen/src/QR/HouseholderQR.h index 24aa96e04..3afaed9e0 100644 --- a/Eigen/src/QR/HouseholderQR.h +++ b/Eigen/src/QR/HouseholderQR.h @@ -231,7 +231,7 @@ struct ei_solve_retval<HouseholderQR<_MatrixType>, Rhs> dec().matrixQR() .corner(TopLeft, rank, rank) - .template triangularView<UpperTriangular>() + .template triangularView<Upper>() .solveInPlace(c.corner(TopLeft, rank, c.cols())); dst.corner(TopLeft, rank, c.cols()) = c.corner(TopLeft, rank, c.cols()); diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h index d2cd8e478..55fac3d12 100644 --- a/Eigen/src/SVD/JacobiSVD.h +++ b/Eigen/src/SVD/JacobiSVD.h @@ -234,7 +234,7 @@ struct ei_svd_precondition_if_more_rows_than_cols<MatrixType, Options, true> if(rows > cols) { FullPivHouseholderQR<MatrixType> qr(matrix); - work_matrix = qr.matrixQR().block(0,0,diagSize,diagSize).template triangularView<UpperTriangular>(); + work_matrix = qr.matrixQR().block(0,0,diagSize,diagSize).template triangularView<Upper>(); if(ComputeU) svd.m_matrixU = qr.matrixQ(); if(ComputeV) svd.m_matrixV = qr.colsPermutation(); @@ -278,7 +278,7 @@ struct ei_svd_precondition_if_more_cols_than_rows<MatrixType, Options, true> MatrixOptions,MaxColsAtCompileTime,MaxRowsAtCompileTime> TransposeTypeWithSameStorageOrder; FullPivHouseholderQR<TransposeTypeWithSameStorageOrder> qr(matrix.adjoint()); - work_matrix = qr.matrixQR().block(0,0,diagSize,diagSize).template triangularView<UpperTriangular>().adjoint(); + work_matrix = qr.matrixQR().block(0,0,diagSize,diagSize).template triangularView<Upper>().adjoint(); if(ComputeV) svd.m_matrixV = qr.matrixQ(); if(ComputeU) svd.m_matrixU = qr.colsPermutation(); return true; diff --git a/Eigen/src/Sparse/CholmodSupport.h b/Eigen/src/Sparse/CholmodSupport.h index f02374d7c..fd33b1507 100644 --- a/Eigen/src/Sparse/CholmodSupport.h +++ b/Eigen/src/Sparse/CholmodSupport.h @@ -76,9 +76,9 @@ cholmod_sparse SparseMatrixBase<Derived>::asCholmodMatrix() if (Derived::Flags & SelfAdjoint) { - if (Derived::Flags & UpperTriangular) + if (Derived::Flags & Upper) res.stype = 1; - else if (Derived::Flags & LowerTriangular) + else if (Derived::Flags & Lower) res.stype = -1; else res.stype = 0; diff --git a/Eigen/src/Sparse/RandomSetter.h b/Eigen/src/Sparse/RandomSetter.h index b34ca19a8..76f24cf0e 100644 --- a/Eigen/src/Sparse/RandomSetter.h +++ b/Eigen/src/Sparse/RandomSetter.h @@ -179,8 +179,8 @@ class RandomSetter SwapStorage = 1 - MapTraits<ScalarWrapper>::IsSorted, TargetRowMajor = (SparseMatrixType::Flags & RowMajorBit) ? 1 : 0, SetterRowMajor = SwapStorage ? 1-TargetRowMajor : TargetRowMajor, - IsUpperTriangular = SparseMatrixType::Flags & UpperTriangularBit, - IsLowerTriangular = SparseMatrixType::Flags & LowerTriangularBit + IsUpper = SparseMatrixType::Flags & Upper, + IsLower = SparseMatrixType::Flags & Lower }; public: @@ -303,8 +303,8 @@ class RandomSetter /** \returns a reference to the coefficient at given coordinates \a row, \a col */ Scalar& operator() (int row, int col) { - ei_assert(((!IsUpperTriangular) || (row<=col)) && "Invalid access to an upper triangular matrix"); - ei_assert(((!IsLowerTriangular) || (col<=row)) && "Invalid access to an upper triangular matrix"); + ei_assert(((!IsUpper) || (row<=col)) && "Invalid access to an upper triangular matrix"); + ei_assert(((!IsLower) || (col<=row)) && "Invalid access to an upper triangular matrix"); const int outer = SetterRowMajor ? row : col; const int inner = SetterRowMajor ? col : row; const int outerMajor = outer >> OuterPacketBits; // index of the packet/map diff --git a/Eigen/src/Sparse/SparseLDLT.h b/Eigen/src/Sparse/SparseLDLT.h index 625357e2b..7dfcf329c 100644 --- a/Eigen/src/Sparse/SparseLDLT.h +++ b/Eigen/src/Sparse/SparseLDLT.h @@ -333,12 +333,12 @@ bool SparseLDLT<MatrixType, Backend>::solveInPlace(MatrixBase<Derived> &b) const return false; if (m_matrix.nonZeros()>0) // otherwise L==I - m_matrix.template triangularView<UnitLowerTriangular>().solveInPlace(b); + m_matrix.template triangularView<UnitLower>().solveInPlace(b); b = b.cwiseQuotient(m_diag); // FIXME should be .adjoint() but it fails to compile... if (m_matrix.nonZeros()>0) // otherwise L==I - m_matrix.transpose().template triangularView<UnitUpperTriangular>().solveInPlace(b); + m_matrix.transpose().template triangularView<UnitUpper>().solveInPlace(b); return true; } diff --git a/Eigen/src/Sparse/SparseLLT.h b/Eigen/src/Sparse/SparseLLT.h index 68ae43022..9018fe0df 100644 --- a/Eigen/src/Sparse/SparseLLT.h +++ b/Eigen/src/Sparse/SparseLLT.h @@ -41,7 +41,7 @@ class SparseLLT protected: typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; - typedef SparseMatrix<Scalar,LowerTriangular> CholMatrixType; + typedef SparseMatrix<Scalar> CholMatrixType; enum { SupernodalFactorIsDirty = 0x10000, @@ -193,15 +193,15 @@ bool SparseLLT<MatrixType, Backend>::solveInPlace(MatrixBase<Derived> &b) const const int size = m_matrix.rows(); ei_assert(size==b.rows()); - m_matrix.template triangularView<LowerTriangular>().solveInPlace(b); + m_matrix.template triangularView<Lower>().solveInPlace(b); // FIXME should be simply .adjoint() but it fails to compile... if (NumTraits<Scalar>::IsComplex) { CholMatrixType aux = m_matrix.conjugate(); - aux.transpose().template triangularView<UpperTriangular>().solveInPlace(b); + aux.transpose().template triangularView<Upper>().solveInPlace(b); } else - m_matrix.transpose().template triangularView<UpperTriangular>().solveInPlace(b); + m_matrix.transpose().template triangularView<Upper>().solveInPlace(b); return true; } diff --git a/Eigen/src/Sparse/SparseLU.h b/Eigen/src/Sparse/SparseLU.h index 2ec6d0e74..0802b0807 100644 --- a/Eigen/src/Sparse/SparseLU.h +++ b/Eigen/src/Sparse/SparseLU.h @@ -47,7 +47,7 @@ class SparseLU protected: typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; - typedef SparseMatrix<Scalar,LowerTriangular> LUMatrixType; + typedef SparseMatrix<Scalar> LUMatrixType; enum { MatrixLUIsDirty = 0x10000 diff --git a/Eigen/src/Sparse/SparseMatrixBase.h b/Eigen/src/Sparse/SparseMatrixBase.h index dc66113f2..9e426c2b9 100644 --- a/Eigen/src/Sparse/SparseMatrixBase.h +++ b/Eigen/src/Sparse/SparseMatrixBase.h @@ -532,8 +532,8 @@ template<typename Derived> class SparseMatrixBase : public AnyMatrixBase<Derived // bool isIdentity(RealScalar prec = dummy_precision<Scalar>()) const; // bool isDiagonal(RealScalar prec = dummy_precision<Scalar>()) const; -// bool isUpperTriangular(RealScalar prec = dummy_precision<Scalar>()) const; -// bool isLowerTriangular(RealScalar prec = dummy_precision<Scalar>()) const; +// bool isUpper(RealScalar prec = dummy_precision<Scalar>()) const; +// bool isLower(RealScalar prec = dummy_precision<Scalar>()) const; // template<typename OtherDerived> // bool isOrthogonal(const MatrixBase<OtherDerived>& other, diff --git a/Eigen/src/Sparse/SparseSelfAdjointView.h b/Eigen/src/Sparse/SparseSelfAdjointView.h index 9f9c5b5e8..039e5c725 100644 --- a/Eigen/src/Sparse/SparseSelfAdjointView.h +++ b/Eigen/src/Sparse/SparseSelfAdjointView.h @@ -31,7 +31,7 @@ * \brief Pseudo expression to manipulate a triangular sparse matrix as a selfadjoint matrix. * * \param MatrixType the type of the dense matrix storing the coefficients - * \param UpLo can be either \c LowerTriangular or \c UpperTriangular + * \param UpLo can be either \c Lower or \c Upper * * This class is an expression of a sefladjoint matrix from a triangular part of a matrix * with given dense storage of the coefficients. It is the return type of MatrixBase::selfadjointView() @@ -168,9 +168,9 @@ class SparseSelfAdjointTimeDenseProduct enum { LhsIsRowMajor = (_Lhs::Flags&RowMajorBit)==RowMajorBit, ProcessFirstHalf = - ((UpLo&(UpperTriangularBit|LowerTriangularBit))==(UpperTriangularBit|LowerTriangularBit)) - || ( (UpLo&UpperTriangularBit) && !LhsIsRowMajor) - || ( (UpLo&LowerTriangularBit) && LhsIsRowMajor), + ((UpLo&(Upper|Lower))==(Upper|Lower)) + || ( (UpLo&Upper) && !LhsIsRowMajor) + || ( (UpLo&Lower) && LhsIsRowMajor), ProcessSecondHalf = !ProcessFirstHalf }; for (int j=0; j<m_lhs.outerSize(); ++j) diff --git a/Eigen/src/Sparse/SparseTriangularView.h b/Eigen/src/Sparse/SparseTriangularView.h index b5eb3d6bd..e713220b9 100644 --- a/Eigen/src/Sparse/SparseTriangularView.h +++ b/Eigen/src/Sparse/SparseTriangularView.h @@ -33,8 +33,8 @@ struct ei_traits<SparseTriangularView<MatrixType,Mode> > template<typename MatrixType, int Mode> class SparseTriangularView : public SparseMatrixBase<SparseTriangularView<MatrixType,Mode> > { - enum { SkipFirst = (Mode==LowerTriangular && !(MatrixType::Flags&RowMajorBit)) - || (Mode==UpperTriangular && (MatrixType::Flags&RowMajorBit)) }; + enum { SkipFirst = (Mode==Lower && !(MatrixType::Flags&RowMajorBit)) + || (Mode==Upper && (MatrixType::Flags&RowMajorBit)) }; public: class InnerIterator; diff --git a/Eigen/src/Sparse/SuperLUSupport.h b/Eigen/src/Sparse/SuperLUSupport.h index 708f177e8..1a765c75b 100644 --- a/Eigen/src/Sparse/SuperLUSupport.h +++ b/Eigen/src/Sparse/SuperLUSupport.h @@ -193,9 +193,9 @@ struct SluMatrix : SuperMatrix res.setScalarType<typename MatrixType::Scalar>(); // FIXME the following is not very accurate - if (MatrixType::Flags & UpperTriangular) + if (MatrixType::Flags & Upper) res.Mtype = SLU_TRU; - if (MatrixType::Flags & LowerTriangular) + if (MatrixType::Flags & Lower) res.Mtype = SLU_TRL; if (MatrixType::Flags & SelfAdjoint) ei_assert(false && "SelfAdjoint matrix shape not supported by SuperLU"); @@ -251,9 +251,9 @@ struct SluMatrixMapHelper<SparseMatrixBase<Derived> > res.setScalarType<typename MatrixType::Scalar>(); // FIXME the following is not very accurate - if (MatrixType::Flags & UpperTriangular) + if (MatrixType::Flags & Upper) res.Mtype = SLU_TRU; - if (MatrixType::Flags & LowerTriangular) + if (MatrixType::Flags & Lower) res.Mtype = SLU_TRL; if (MatrixType::Flags & SelfAdjoint) ei_assert(false && "SelfAdjoint matrix shape not supported by SuperLU"); @@ -298,8 +298,8 @@ class SparseLU<MatrixType,SuperLU> : public SparseLU<MatrixType> typedef Matrix<Scalar,Dynamic,1> Vector; typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType; typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType; - typedef SparseMatrix<Scalar,LowerTriangular|UnitDiagBit> LMatrixType; - typedef SparseMatrix<Scalar,UpperTriangular> UMatrixType; + typedef SparseMatrix<Scalar,Lower|UnitDiag> LMatrixType; + typedef SparseMatrix<Scalar,Upper> UMatrixType; using Base::m_flags; using Base::m_status; diff --git a/Eigen/src/Sparse/TaucsSupport.h b/Eigen/src/Sparse/TaucsSupport.h index cbe3d0036..f2a9d19dc 100644 --- a/Eigen/src/Sparse/TaucsSupport.h +++ b/Eigen/src/Sparse/TaucsSupport.h @@ -50,13 +50,13 @@ taucs_ccs_matrix SparseMatrixBase<Derived>::asTaucsMatrix() ei_assert(false && "Scalar type not supported by TAUCS"); } - if (Flags & UpperTriangular) + if (Flags & Upper) res.flags |= TAUCS_UPPER; - if (Flags & LowerTriangular) + if (Flags & Lower) res.flags |= TAUCS_LOWER; if (Flags & SelfAdjoint) res.flags |= (NumTraits<Scalar>::IsComplex ? TAUCS_HERMITIAN : TAUCS_SYMMETRIC); - else if ((Flags & UpperTriangular) || (Flags & LowerTriangular)) + else if ((Flags & Upper) || (Flags & Lower)) res.flags |= TAUCS_TRIANGULAR; return res; diff --git a/Eigen/src/Sparse/TriangularSolver.h b/Eigen/src/Sparse/TriangularSolver.h index 5b761e6ad..039424bf0 100644 --- a/Eigen/src/Sparse/TriangularSolver.h +++ b/Eigen/src/Sparse/TriangularSolver.h @@ -26,17 +26,17 @@ #define EIGEN_SPARSETRIANGULARSOLVER_H template<typename Lhs, typename Rhs, int Mode, - int UpLo = (Mode & LowerTriangularBit) - ? LowerTriangular - : (Mode & UpperTriangularBit) - ? UpperTriangular + int UpLo = (Mode & Lower) + ? Lower + : (Mode & Upper) + ? Upper : -1, int StorageOrder = int(ei_traits<Lhs>::Flags) & RowMajorBit> struct ei_sparse_solve_triangular_selector; // forward substitution, row-major template<typename Lhs, typename Rhs, int Mode> -struct ei_sparse_solve_triangular_selector<Lhs,Rhs,Mode,LowerTriangular,RowMajor> +struct ei_sparse_solve_triangular_selector<Lhs,Rhs,Mode,Lower,RowMajor> { typedef typename Rhs::Scalar Scalar; static void run(const Lhs& lhs, Rhs& other) @@ -56,7 +56,7 @@ struct ei_sparse_solve_triangular_selector<Lhs,Rhs,Mode,LowerTriangular,RowMajor break; tmp -= lastVal * other.coeff(lastIndex,col); } - if (Mode & UnitDiagBit) + if (Mode & UnitDiag) other.coeffRef(i,col) = tmp; else { @@ -70,7 +70,7 @@ struct ei_sparse_solve_triangular_selector<Lhs,Rhs,Mode,LowerTriangular,RowMajor // backward substitution, row-major template<typename Lhs, typename Rhs, int Mode> -struct ei_sparse_solve_triangular_selector<Lhs,Rhs,Mode,UpperTriangular,RowMajor> +struct ei_sparse_solve_triangular_selector<Lhs,Rhs,Mode,Upper,RowMajor> { typedef typename Rhs::Scalar Scalar; static void run(const Lhs& lhs, Rhs& other) @@ -88,7 +88,7 @@ struct ei_sparse_solve_triangular_selector<Lhs,Rhs,Mode,UpperTriangular,RowMajor tmp -= it.value() * other.coeff(it.index(),col); } - if (Mode & UnitDiagBit) + if (Mode & UnitDiag) other.coeffRef(i,col) = tmp; else { @@ -103,7 +103,7 @@ struct ei_sparse_solve_triangular_selector<Lhs,Rhs,Mode,UpperTriangular,RowMajor // forward substitution, col-major template<typename Lhs, typename Rhs, int Mode> -struct ei_sparse_solve_triangular_selector<Lhs,Rhs,Mode,LowerTriangular,ColMajor> +struct ei_sparse_solve_triangular_selector<Lhs,Rhs,Mode,Lower,ColMajor> { typedef typename Rhs::Scalar Scalar; static void run(const Lhs& lhs, Rhs& other) @@ -116,7 +116,7 @@ struct ei_sparse_solve_triangular_selector<Lhs,Rhs,Mode,LowerTriangular,ColMajor if (tmp!=Scalar(0)) // optimization when other is actually sparse { typename Lhs::InnerIterator it(lhs, i); - if(!(Mode & UnitDiagBit)) + if(!(Mode & UnitDiag)) { ei_assert(it.index()==i); tmp /= it.value(); @@ -133,7 +133,7 @@ struct ei_sparse_solve_triangular_selector<Lhs,Rhs,Mode,LowerTriangular,ColMajor // backward substitution, col-major template<typename Lhs, typename Rhs, int Mode> -struct ei_sparse_solve_triangular_selector<Lhs,Rhs,Mode,UpperTriangular,ColMajor> +struct ei_sparse_solve_triangular_selector<Lhs,Rhs,Mode,Upper,ColMajor> { typedef typename Rhs::Scalar Scalar; static void run(const Lhs& lhs, Rhs& other) @@ -145,7 +145,7 @@ struct ei_sparse_solve_triangular_selector<Lhs,Rhs,Mode,UpperTriangular,ColMajor Scalar& tmp = other.coeffRef(i,col); if (tmp!=Scalar(0)) // optimization when other is actually sparse { - if(!(Mode & UnitDiagBit)) + if(!(Mode & UnitDiag)) { // FIXME lhs.coeff(i,i) might not be always efficient while it must simply be the // last element of the column ! @@ -166,8 +166,8 @@ void SparseTriangularView<ExpressionType,Mode>::solveInPlace(MatrixBase<OtherDer { ei_assert(m_matrix.cols() == m_matrix.rows()); ei_assert(m_matrix.cols() == other.rows()); - ei_assert(!(Mode & ZeroDiagBit)); - ei_assert(Mode & (UpperTriangularBit|LowerTriangularBit)); + ei_assert(!(Mode & ZeroDiag)); + ei_assert(Mode & (Upper|Lower)); enum { copy = ei_traits<OtherDerived>::Flags & RowMajorBit }; @@ -194,10 +194,10 @@ SparseTriangularView<ExpressionType,Mode>::solve(const MatrixBase<OtherDerived>& // pure sparse path template<typename Lhs, typename Rhs, int Mode, - int UpLo = (Mode & LowerTriangularBit) - ? LowerTriangular - : (Mode & UpperTriangularBit) - ? UpperTriangular + int UpLo = (Mode & Lower) + ? Lower + : (Mode & Upper) + ? Upper : -1, int StorageOrder = int(Lhs::Flags) & (RowMajorBit)> struct ei_sparse_solve_triangular_sparse_selector; @@ -209,7 +209,7 @@ struct ei_sparse_solve_triangular_sparse_selector<Lhs,Rhs,Mode,UpLo,ColMajor> typedef typename Rhs::Scalar Scalar; static void run(const Lhs& lhs, Rhs& other) { - const bool IsLowerTriangular = (UpLo==LowerTriangular); + const bool IsLower = (UpLo==Lower); AmbiVector<Scalar> tempVector(other.rows()*2); tempVector.setBounds(0,other.rows()); @@ -227,9 +227,9 @@ struct ei_sparse_solve_triangular_sparse_selector<Lhs,Rhs,Mode,UpLo,ColMajor> tempVector.coeffRef(rhsIt.index()) = rhsIt.value(); } - for(int i=IsLowerTriangular?0:lhs.cols()-1; - IsLowerTriangular?i<lhs.cols():i>=0; - i+=IsLowerTriangular?1:-1) + for(int i=IsLower?0:lhs.cols()-1; + IsLower?i<lhs.cols():i>=0; + i+=IsLower?1:-1) { tempVector.restart(); Scalar& ci = tempVector.coeffRef(i); @@ -237,9 +237,9 @@ struct ei_sparse_solve_triangular_sparse_selector<Lhs,Rhs,Mode,UpLo,ColMajor> { // find typename Lhs::InnerIterator it(lhs, i); - if(!(Mode & UnitDiagBit)) + if(!(Mode & UnitDiag)) { - if (IsLowerTriangular) + if (IsLower) { ei_assert(it.index()==i); ci /= it.value(); @@ -248,7 +248,7 @@ struct ei_sparse_solve_triangular_sparse_selector<Lhs,Rhs,Mode,UpLo,ColMajor> ci /= lhs.coeff(i,i); } tempVector.restart(); - if (IsLowerTriangular) + if (IsLower) { if (it.index()==i) ++it; @@ -287,8 +287,8 @@ void SparseTriangularView<ExpressionType,Mode>::solveInPlace(SparseMatrixBase<Ot { ei_assert(m_matrix.cols() == m_matrix.rows()); ei_assert(m_matrix.cols() == other.rows()); - ei_assert(!(Mode & ZeroDiagBit)); - ei_assert(Mode & (UpperTriangularBit|LowerTriangularBit)); + ei_assert(!(Mode & ZeroDiag)); + ei_assert(Mode & (Upper|Lower)); // enum { copy = ei_traits<OtherDerived>::Flags & RowMajorBit }; @@ -311,7 +311,7 @@ template<typename Derived> template<typename OtherDerived> void SparseMatrixBase<Derived>::solveTriangularInPlace(MatrixBase<OtherDerived>& other) const { - this->template triangular<Flags&(UpperTriangularBit|LowerTriangularBit)>().solveInPlace(other); + this->template triangular<Flags&(Upper|Lower)>().solveInPlace(other); } /** \deprecated */ diff --git a/Eigen/src/Sparse/UmfPackSupport.h b/Eigen/src/Sparse/UmfPackSupport.h index 1c7a7aaa8..cde924964 100644 --- a/Eigen/src/Sparse/UmfPackSupport.h +++ b/Eigen/src/Sparse/UmfPackSupport.h @@ -127,8 +127,8 @@ class SparseLU<MatrixType,UmfPack> : public SparseLU<MatrixType> typedef Matrix<Scalar,Dynamic,1> Vector; typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType; typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType; - typedef SparseMatrix<Scalar,LowerTriangular|UnitDiagBit> LMatrixType; - typedef SparseMatrix<Scalar,UpperTriangular> UMatrixType; + typedef SparseMatrix<Scalar,Lower|UnitDiagBit> LMatrixType; + typedef SparseMatrix<Scalar,Upper> UMatrixType; using Base::m_flags; using Base::m_status; diff --git a/doc/A05_PortingFrom2To3.dox b/doc/A05_PortingFrom2To3.dox index e2dac4525..3c1952ebc 100644 --- a/doc/A05_PortingFrom2To3.dox +++ b/doc/A05_PortingFrom2To3.dox @@ -29,13 +29,30 @@ vec.head(length) vec.head<length>() vec.tail(length) vec.tail<length>() -\endcode</td><td>This is a simple search and replace change.</td></tr> +\endcode</td><td>Trivial "search and replace".</td></tr> <tr><td colspan="3"></td></tr> <tr><td>\code mat.cwise().XXX()\endcode</td><td>\code mat.array().XXX()\endcode</td><td>See \ref CoefficientWiseOperations. </td></tr> <tr><td colspan="3"></td></tr> <tr><td>\code c = (a * b).lazy();\endcode</td><td>\code c.noalias() = a * b;\endcode</td><td>See \ref LazyVsNoalias. </td></tr> <tr><td colspan="3"></td></tr> <tr><td>\code A.triangularSolveInPlace<XXX>(X);\endcode</td><td>\code A.triangularView<XXX>().solveInPlace(X);\endcode</td><td></td></tr> +<tr><td colspan="3"></td></tr> +<tr><td>\code +UpperTriangular +LowerTriangular +UnitUpperTriangular +UnitLowerTriangular +StrictlyUpperTriangular +StrictlyLowerTriangular +\endcode</td><td>\code +Upper +Lower +UnitUpper +UnitLower +StrictlyUpper +StrictlyLower +\endcode</td> +<td>Trivial "search and replace".</td></tr> </table> \section CoefficientWiseOperations Coefficient wise operations diff --git a/test/bandmatrix.cpp b/test/bandmatrix.cpp index e243dffe5..dc54812b9 100644 --- a/test/bandmatrix.cpp +++ b/test/bandmatrix.cpp @@ -64,8 +64,8 @@ template<typename MatrixType> void bandmatrix(const MatrixType& _m) int a = std::max(0,cols-d-supers); int b = std::max(0,rows-d-subs); if(a>0) dm1.block(0,d+supers,rows,a).setZero(); - dm1.block(0,supers+1,cols-supers-1-a,cols-supers-1-a).template triangularView<UpperTriangular>().setZero(); - dm1.block(subs+1,0,rows-subs-1-b,rows-subs-1-b).template triangularView<LowerTriangular>().setZero(); + dm1.block(0,supers+1,cols-supers-1-a,cols-supers-1-a).template triangularView<Upper>().setZero(); + dm1.block(subs+1,0,rows-subs-1-b,rows-subs-1-b).template triangularView<Lower>().setZero(); if(b>0) dm1.block(d+subs,0,b,cols).setZero(); //std::cerr << m.m_data << "\n\n" << m.toDense() << "\n\n" << dm1 << "\n\n"; VERIFY_IS_APPROX(dm1,m.toDenseMatrix()); diff --git a/test/cholesky.cpp b/test/cholesky.cpp index c658b902c..1bb808d20 100644 --- a/test/cholesky.cpp +++ b/test/cholesky.cpp @@ -58,12 +58,12 @@ template<typename MatrixType> void cholesky(const MatrixType& m) symm += a1 * a1.adjoint(); } - SquareMatrixType symmUp = symm.template triangularView<UpperTriangular>(); - SquareMatrixType symmLo = symm.template triangularView<LowerTriangular>(); + SquareMatrixType symmUp = symm.template triangularView<Upper>(); + SquareMatrixType symmLo = symm.template triangularView<Lower>(); // to test if really Cholesky only uses the upper triangular part, uncomment the following // FIXME: currently that fails !! - //symm.template part<StrictlyLowerTriangular>().setZero(); + //symm.template part<StrictlyLower>().setZero(); #ifdef HAS_GSL // if (ei_is_same_type<RealScalar,double>::ret) @@ -94,7 +94,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m) #endif { - LLT<SquareMatrixType,LowerTriangular> chollo(symmLo); + LLT<SquareMatrixType,Lower> chollo(symmLo); VERIFY_IS_APPROX(symm, chollo.matrixL().toDenseMatrix() * chollo.matrixL().adjoint().toDenseMatrix()); vecX = chollo.solve(vecB); VERIFY_IS_APPROX(symm * vecX, vecB); @@ -102,7 +102,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m) VERIFY_IS_APPROX(symm * matX, matB); // test the upper mode - LLT<SquareMatrixType,UpperTriangular> cholup(symmUp); + LLT<SquareMatrixType,Upper> cholup(symmUp); VERIFY_IS_APPROX(symm, cholup.matrixL().toDenseMatrix() * cholup.matrixL().adjoint().toDenseMatrix()); vecX = cholup.solve(vecB); VERIFY_IS_APPROX(symm * vecX, vecB); diff --git a/test/product_notemporary.cpp b/test/product_notemporary.cpp index 534d23d84..1da91e5d7 100644 --- a/test/product_notemporary.cpp +++ b/test/product_notemporary.cpp @@ -85,31 +85,31 @@ template<typename MatrixType> void product_notemporary(const MatrixType& m) // NOTE this is because the Block expression is not handled yet by our expression analyser VERIFY_EVALUATION_COUNT(( m3.block(r0,r0,r1,r1).noalias() = s1 * m1.block(r0,c0,r1,c1) * (s1*m2).block(c0,r0,c1,r1) ), 1); - VERIFY_EVALUATION_COUNT( m3.noalias() -= (s1 * m1).template triangularView<LowerTriangular>() * m2, 0); - VERIFY_EVALUATION_COUNT( rm3.noalias() = (s1 * m1.adjoint()).template triangularView<UpperTriangular>() * (m2+m2), 1); - VERIFY_EVALUATION_COUNT( rm3.noalias() = (s1 * m1.adjoint()).template triangularView<UnitUpperTriangular>() * m2.adjoint(), 0); + VERIFY_EVALUATION_COUNT( m3.noalias() -= (s1 * m1).template triangularView<Lower>() * m2, 0); + VERIFY_EVALUATION_COUNT( rm3.noalias() = (s1 * m1.adjoint()).template triangularView<Upper>() * (m2+m2), 1); + VERIFY_EVALUATION_COUNT( rm3.noalias() = (s1 * m1.adjoint()).template triangularView<UnitUpper>() * m2.adjoint(), 0); - VERIFY_EVALUATION_COUNT( rm3.col(c0).noalias() = (s1 * m1.adjoint()).template triangularView<UnitUpperTriangular>() * (s2*m2.row(c0)).adjoint(), 0); + VERIFY_EVALUATION_COUNT( rm3.col(c0).noalias() = (s1 * m1.adjoint()).template triangularView<UnitUpper>() * (s2*m2.row(c0)).adjoint(), 0); - VERIFY_EVALUATION_COUNT( m1.template triangularView<LowerTriangular>().solveInPlace(m3), 0); - VERIFY_EVALUATION_COUNT( m1.adjoint().template triangularView<LowerTriangular>().solveInPlace(m3.transpose()), 0); + VERIFY_EVALUATION_COUNT( m1.template triangularView<Lower>().solveInPlace(m3), 0); + VERIFY_EVALUATION_COUNT( m1.adjoint().template triangularView<Lower>().solveInPlace(m3.transpose()), 0); - VERIFY_EVALUATION_COUNT( m3.noalias() -= (s1 * m1).adjoint().template selfadjointView<LowerTriangular>() * (-m2*s3).adjoint(), 0); - VERIFY_EVALUATION_COUNT( m3.noalias() = s2 * m2.adjoint() * (s1 * m1.adjoint()).template selfadjointView<UpperTriangular>(), 0); - VERIFY_EVALUATION_COUNT( rm3.noalias() = (s1 * m1.adjoint()).template selfadjointView<LowerTriangular>() * m2.adjoint(), 0); + VERIFY_EVALUATION_COUNT( m3.noalias() -= (s1 * m1).adjoint().template selfadjointView<Lower>() * (-m2*s3).adjoint(), 0); + VERIFY_EVALUATION_COUNT( m3.noalias() = s2 * m2.adjoint() * (s1 * m1.adjoint()).template selfadjointView<Upper>(), 0); + VERIFY_EVALUATION_COUNT( rm3.noalias() = (s1 * m1.adjoint()).template selfadjointView<Lower>() * m2.adjoint(), 0); - VERIFY_EVALUATION_COUNT( m3.col(c0).noalias() = (s1 * m1).adjoint().template selfadjointView<LowerTriangular>() * (-m2.row(c0)*s3).adjoint(), 0); - VERIFY_EVALUATION_COUNT( m3.col(c0).noalias() -= (s1 * m1).adjoint().template selfadjointView<UpperTriangular>() * (-m2.row(c0)*s3).adjoint(), 0); + VERIFY_EVALUATION_COUNT( m3.col(c0).noalias() = (s1 * m1).adjoint().template selfadjointView<Lower>() * (-m2.row(c0)*s3).adjoint(), 0); + VERIFY_EVALUATION_COUNT( m3.col(c0).noalias() -= (s1 * m1).adjoint().template selfadjointView<Upper>() * (-m2.row(c0)*s3).adjoint(), 0); - VERIFY_EVALUATION_COUNT( m3.block(r0,c0,r1,c1).noalias() += m1.block(r0,r0,r1,r1).template selfadjointView<UpperTriangular>() * (s1*m2.block(r0,c0,r1,c1)), 0); - VERIFY_EVALUATION_COUNT( m3.block(r0,c0,r1,c1).noalias() = m1.block(r0,r0,r1,r1).template selfadjointView<UpperTriangular>() * m2.block(r0,c0,r1,c1), 0); + VERIFY_EVALUATION_COUNT( m3.block(r0,c0,r1,c1).noalias() += m1.block(r0,r0,r1,r1).template selfadjointView<Upper>() * (s1*m2.block(r0,c0,r1,c1)), 0); + VERIFY_EVALUATION_COUNT( m3.block(r0,c0,r1,c1).noalias() = m1.block(r0,r0,r1,r1).template selfadjointView<Upper>() * m2.block(r0,c0,r1,c1), 0); - VERIFY_EVALUATION_COUNT( m3.template selfadjointView<LowerTriangular>().rankUpdate(m2.adjoint()), 0); + VERIFY_EVALUATION_COUNT( m3.template selfadjointView<Lower>().rankUpdate(m2.adjoint()), 0); m3.resize(1,1); - VERIFY_EVALUATION_COUNT( m3.noalias() = m1.block(r0,r0,r1,r1).template selfadjointView<LowerTriangular>() * m2.block(r0,c0,r1,c1), 0); + VERIFY_EVALUATION_COUNT( m3.noalias() = m1.block(r0,r0,r1,r1).template selfadjointView<Lower>() * m2.block(r0,c0,r1,c1), 0); m3.resize(1,1); - VERIFY_EVALUATION_COUNT( m3.noalias() = m1.block(r0,r0,r1,r1).template triangularView<UnitUpperTriangular>() * m2.block(r0,c0,r1,c1), 0); + VERIFY_EVALUATION_COUNT( m3.noalias() = m1.block(r0,r0,r1,r1).template triangularView<UnitUpper>() * m2.block(r0,c0,r1,c1), 0); } void test_product_notemporary() diff --git a/test/product_selfadjoint.cpp b/test/product_selfadjoint.cpp index 2f3833a02..2027fc8e5 100644 --- a/test/product_selfadjoint.cpp +++ b/test/product_selfadjoint.cpp @@ -53,25 +53,25 @@ template<typename MatrixType> void product_selfadjoint(const MatrixType& m) m1 = (m1.adjoint() + m1).eval(); // rank2 update - m2 = m1.template triangularView<LowerTriangular>(); - m2.template selfadjointView<LowerTriangular>().rankUpdate(v1,v2); - VERIFY_IS_APPROX(m2, (m1 + v1 * v2.adjoint()+ v2 * v1.adjoint()).template triangularView<LowerTriangular>().toDenseMatrix()); + m2 = m1.template triangularView<Lower>(); + m2.template selfadjointView<Lower>().rankUpdate(v1,v2); + VERIFY_IS_APPROX(m2, (m1 + v1 * v2.adjoint()+ v2 * v1.adjoint()).template triangularView<Lower>().toDenseMatrix()); - m2 = m1.template triangularView<UpperTriangular>(); - m2.template selfadjointView<UpperTriangular>().rankUpdate(-v1,s2*v2,s3); - VERIFY_IS_APPROX(m2, (m1 + (-s2*s3) * (v1 * v2.adjoint()+ v2 * v1.adjoint())).template triangularView<UpperTriangular>().toDenseMatrix()); + m2 = m1.template triangularView<Upper>(); + m2.template selfadjointView<Upper>().rankUpdate(-v1,s2*v2,s3); + VERIFY_IS_APPROX(m2, (m1 + (-s2*s3) * (v1 * v2.adjoint()+ v2 * v1.adjoint())).template triangularView<Upper>().toDenseMatrix()); - m2 = m1.template triangularView<UpperTriangular>(); - m2.template selfadjointView<UpperTriangular>().rankUpdate(-r1.adjoint(),r2.adjoint()*s3,s1); - VERIFY_IS_APPROX(m2, (m1 + (-s3*s1) * (r1.adjoint() * r2 + r2.adjoint() * r1)).template triangularView<UpperTriangular>().toDenseMatrix()); + m2 = m1.template triangularView<Upper>(); + m2.template selfadjointView<Upper>().rankUpdate(-r1.adjoint(),r2.adjoint()*s3,s1); + VERIFY_IS_APPROX(m2, (m1 + (-s3*s1) * (r1.adjoint() * r2 + r2.adjoint() * r1)).template triangularView<Upper>().toDenseMatrix()); if (rows>1) { - m2 = m1.template triangularView<LowerTriangular>(); - m2.block(1,1,rows-1,cols-1).template selfadjointView<LowerTriangular>().rankUpdate(v1.tail(rows-1),v2.head(cols-1)); + m2 = m1.template triangularView<Lower>(); + m2.block(1,1,rows-1,cols-1).template selfadjointView<Lower>().rankUpdate(v1.tail(rows-1),v2.head(cols-1)); m3 = m1; m3.block(1,1,rows-1,cols-1) += v1.tail(rows-1) * v2.head(cols-1).adjoint()+ v2.head(cols-1) * v1.tail(rows-1).adjoint(); - VERIFY_IS_APPROX(m2, m3.template triangularView<LowerTriangular>().toDenseMatrix()); + VERIFY_IS_APPROX(m2, m3.template triangularView<Lower>().toDenseMatrix()); } } diff --git a/test/product_symm.cpp b/test/product_symm.cpp index a9e055bd3..a0d80080f 100644 --- a/test/product_symm.cpp +++ b/test/product_symm.cpp @@ -28,10 +28,10 @@ template<int OtherSize> struct symm_extra { template<typename M1, typename M2, typename Scalar> static void run(M1& m1, M1& m2, M2& rhs2, M2& rhs22, M2& rhs23, Scalar s1, Scalar s2) { - m2 = m1.template triangularView<LowerTriangular>(); - VERIFY_IS_APPROX(rhs22 = (rhs2) * (m2).template selfadjointView<LowerTriangular>(), + m2 = m1.template triangularView<Lower>(); + VERIFY_IS_APPROX(rhs22 = (rhs2) * (m2).template selfadjointView<Lower>(), rhs23 = (rhs2) * (m1)); - VERIFY_IS_APPROX(rhs22 = (s2*rhs2) * (s1*m2).template selfadjointView<LowerTriangular>(), + VERIFY_IS_APPROX(rhs22 = (s2*rhs2) * (s1*m2).template selfadjointView<Lower>(), rhs23 = (s2*rhs2) * (s1*m1)); } }; @@ -65,38 +65,38 @@ template<typename Scalar, int Size, int OtherSize> void symm(int size = Size, in Scalar s1 = ei_random<Scalar>(), s2 = ei_random<Scalar>(); - m2 = m1.template triangularView<LowerTriangular>(); - VERIFY_IS_APPROX(rhs12 = (s1*m2).template selfadjointView<LowerTriangular>() * (s2*rhs1), + m2 = m1.template triangularView<Lower>(); + VERIFY_IS_APPROX(rhs12 = (s1*m2).template selfadjointView<Lower>() * (s2*rhs1), rhs13 = (s1*m1) * (s2*rhs1)); - m2 = m1.template triangularView<UpperTriangular>(); rhs12.setRandom(); rhs13 = rhs12; - VERIFY_IS_APPROX(rhs12 += (s1*m2).template selfadjointView<UpperTriangular>() * (s2*rhs1), + m2 = m1.template triangularView<Upper>(); rhs12.setRandom(); rhs13 = rhs12; + VERIFY_IS_APPROX(rhs12 += (s1*m2).template selfadjointView<Upper>() * (s2*rhs1), rhs13 += (s1*m1) * (s2*rhs1)); - m2 = m1.template triangularView<LowerTriangular>(); - VERIFY_IS_APPROX(rhs12 = (s1*m2).template selfadjointView<LowerTriangular>() * (s2*rhs2.adjoint()), + m2 = m1.template triangularView<Lower>(); + VERIFY_IS_APPROX(rhs12 = (s1*m2).template selfadjointView<Lower>() * (s2*rhs2.adjoint()), rhs13 = (s1*m1) * (s2*rhs2.adjoint())); - m2 = m1.template triangularView<UpperTriangular>(); - VERIFY_IS_APPROX(rhs12 = (s1*m2).template selfadjointView<UpperTriangular>() * (s2*rhs2.adjoint()), + m2 = m1.template triangularView<Upper>(); + VERIFY_IS_APPROX(rhs12 = (s1*m2).template selfadjointView<Upper>() * (s2*rhs2.adjoint()), rhs13 = (s1*m1) * (s2*rhs2.adjoint())); - m2 = m1.template triangularView<UpperTriangular>(); - VERIFY_IS_APPROX(rhs12 = (s1*m2.adjoint()).template selfadjointView<LowerTriangular>() * (s2*rhs2.adjoint()), + m2 = m1.template triangularView<Upper>(); + VERIFY_IS_APPROX(rhs12 = (s1*m2.adjoint()).template selfadjointView<Lower>() * (s2*rhs2.adjoint()), rhs13 = (s1*m1.adjoint()) * (s2*rhs2.adjoint())); // test row major = <...> - m2 = m1.template triangularView<LowerTriangular>(); rhs12.setRandom(); rhs13 = rhs12; - VERIFY_IS_APPROX(rhs12 -= (s1*m2).template selfadjointView<LowerTriangular>() * (s2*rhs3), + m2 = m1.template triangularView<Lower>(); rhs12.setRandom(); rhs13 = rhs12; + VERIFY_IS_APPROX(rhs12 -= (s1*m2).template selfadjointView<Lower>() * (s2*rhs3), rhs13 -= (s1*m1) * (s2 * rhs3)); - m2 = m1.template triangularView<UpperTriangular>(); - VERIFY_IS_APPROX(rhs12 = (s1*m2.adjoint()).template selfadjointView<LowerTriangular>() * (s2*rhs3).conjugate(), + m2 = m1.template triangularView<Upper>(); + VERIFY_IS_APPROX(rhs12 = (s1*m2.adjoint()).template selfadjointView<Lower>() * (s2*rhs3).conjugate(), rhs13 = (s1*m1.adjoint()) * (s2*rhs3).conjugate()); - m2 = m1.template triangularView<UpperTriangular>(); rhs13 = rhs12; - VERIFY_IS_APPROX(rhs12.noalias() += s1 * ((m2.adjoint()).template selfadjointView<LowerTriangular>() * (s2*rhs3).conjugate()), + m2 = m1.template triangularView<Upper>(); rhs13 = rhs12; + VERIFY_IS_APPROX(rhs12.noalias() += s1 * ((m2.adjoint()).template selfadjointView<Lower>() * (s2*rhs3).conjugate()), rhs13 += (s1*m1.adjoint()) * (s2*rhs3).conjugate()); // test matrix * selfadjoint diff --git a/test/product_syrk.cpp b/test/product_syrk.cpp index 9f6aec0e2..e597ac88a 100644 --- a/test/product_syrk.cpp +++ b/test/product_syrk.cpp @@ -45,28 +45,28 @@ template<typename MatrixType> void syrk(const MatrixType& m) Scalar s1 = ei_random<Scalar>(); m2.setZero(); - VERIFY_IS_APPROX((m2.template selfadjointView<LowerTriangular>().rankUpdate(rhs2,s1)._expression()), - ((s1 * rhs2 * rhs2.adjoint()).eval().template triangularView<LowerTriangular>().toDenseMatrix())); + VERIFY_IS_APPROX((m2.template selfadjointView<Lower>().rankUpdate(rhs2,s1)._expression()), + ((s1 * rhs2 * rhs2.adjoint()).eval().template triangularView<Lower>().toDenseMatrix())); m2.setZero(); - VERIFY_IS_APPROX(m2.template selfadjointView<UpperTriangular>().rankUpdate(rhs2,s1)._expression(), - (s1 * rhs2 * rhs2.adjoint()).eval().template triangularView<UpperTriangular>().toDenseMatrix()); + VERIFY_IS_APPROX(m2.template selfadjointView<Upper>().rankUpdate(rhs2,s1)._expression(), + (s1 * rhs2 * rhs2.adjoint()).eval().template triangularView<Upper>().toDenseMatrix()); m2.setZero(); - VERIFY_IS_APPROX(m2.template selfadjointView<LowerTriangular>().rankUpdate(rhs1.adjoint(),s1)._expression(), - (s1 * rhs1.adjoint() * rhs1).eval().template triangularView<LowerTriangular>().toDenseMatrix()); + VERIFY_IS_APPROX(m2.template selfadjointView<Lower>().rankUpdate(rhs1.adjoint(),s1)._expression(), + (s1 * rhs1.adjoint() * rhs1).eval().template triangularView<Lower>().toDenseMatrix()); m2.setZero(); - VERIFY_IS_APPROX(m2.template selfadjointView<UpperTriangular>().rankUpdate(rhs1.adjoint(),s1)._expression(), - (s1 * rhs1.adjoint() * rhs1).eval().template triangularView<UpperTriangular>().toDenseMatrix()); + VERIFY_IS_APPROX(m2.template selfadjointView<Upper>().rankUpdate(rhs1.adjoint(),s1)._expression(), + (s1 * rhs1.adjoint() * rhs1).eval().template triangularView<Upper>().toDenseMatrix()); m2.setZero(); - VERIFY_IS_APPROX(m2.template selfadjointView<LowerTriangular>().rankUpdate(rhs3.adjoint(),s1)._expression(), - (s1 * rhs3.adjoint() * rhs3).eval().template triangularView<LowerTriangular>().toDenseMatrix()); + VERIFY_IS_APPROX(m2.template selfadjointView<Lower>().rankUpdate(rhs3.adjoint(),s1)._expression(), + (s1 * rhs3.adjoint() * rhs3).eval().template triangularView<Lower>().toDenseMatrix()); m2.setZero(); - VERIFY_IS_APPROX(m2.template selfadjointView<UpperTriangular>().rankUpdate(rhs3.adjoint(),s1)._expression(), - (s1 * rhs3.adjoint() * rhs3).eval().template triangularView<UpperTriangular>().toDenseMatrix()); + VERIFY_IS_APPROX(m2.template selfadjointView<Upper>().rankUpdate(rhs3.adjoint(),s1)._expression(), + (s1 * rhs3.adjoint() * rhs3).eval().template triangularView<Upper>().toDenseMatrix()); } void test_product_syrk() diff --git a/test/product_trmm.cpp b/test/product_trmm.cpp index 5f92391e6..69e97f7aa 100644 --- a/test/product_trmm.cpp +++ b/test/product_trmm.cpp @@ -36,27 +36,27 @@ template<typename Scalar> void trmm(int size,int othersize) s2 = ei_random<Scalar>(); tri.setRandom(); - loTri = tri.template triangularView<LowerTriangular>(); - upTri = tri.template triangularView<UpperTriangular>(); + loTri = tri.template triangularView<Lower>(); + upTri = tri.template triangularView<Upper>(); ge1.setRandom(); ge2.setRandom(); - VERIFY_IS_APPROX( ge3 = tri.template triangularView<LowerTriangular>() * ge1, loTri * ge1); - VERIFY_IS_APPROX(rge3 = tri.template triangularView<LowerTriangular>() * ge1, loTri * ge1); - VERIFY_IS_APPROX( ge3 = tri.template triangularView<UpperTriangular>() * ge1, upTri * ge1); - VERIFY_IS_APPROX(rge3 = tri.template triangularView<UpperTriangular>() * ge1, upTri * ge1); - VERIFY_IS_APPROX( ge3 = (s1*tri.adjoint()).template triangularView<UpperTriangular>() * (s2*ge1), s1*loTri.adjoint() * (s2*ge1)); - VERIFY_IS_APPROX(rge3 = tri.adjoint().template triangularView<UpperTriangular>() * ge1, loTri.adjoint() * ge1); - VERIFY_IS_APPROX( ge3 = tri.adjoint().template triangularView<LowerTriangular>() * ge1, upTri.adjoint() * ge1); - VERIFY_IS_APPROX(rge3 = tri.adjoint().template triangularView<LowerTriangular>() * ge1, upTri.adjoint() * ge1); - VERIFY_IS_APPROX( ge3 = tri.template triangularView<LowerTriangular>() * ge2.adjoint(), loTri * ge2.adjoint()); - VERIFY_IS_APPROX(rge3 = tri.template triangularView<LowerTriangular>() * ge2.adjoint(), loTri * ge2.adjoint()); - VERIFY_IS_APPROX( ge3 = tri.template triangularView<UpperTriangular>() * ge2.adjoint(), upTri * ge2.adjoint()); - VERIFY_IS_APPROX(rge3 = tri.template triangularView<UpperTriangular>() * ge2.adjoint(), upTri * ge2.adjoint()); - VERIFY_IS_APPROX( ge3 = (s1*tri).adjoint().template triangularView<UpperTriangular>() * ge2.adjoint(), ei_conj(s1) * loTri.adjoint() * ge2.adjoint()); - VERIFY_IS_APPROX(rge3 = tri.adjoint().template triangularView<UpperTriangular>() * ge2.adjoint(), loTri.adjoint() * ge2.adjoint()); - VERIFY_IS_APPROX( ge3 = tri.adjoint().template triangularView<LowerTriangular>() * ge2.adjoint(), upTri.adjoint() * ge2.adjoint()); - VERIFY_IS_APPROX(rge3 = tri.adjoint().template triangularView<LowerTriangular>() * ge2.adjoint(), upTri.adjoint() * ge2.adjoint()); + VERIFY_IS_APPROX( ge3 = tri.template triangularView<Lower>() * ge1, loTri * ge1); + VERIFY_IS_APPROX(rge3 = tri.template triangularView<Lower>() * ge1, loTri * ge1); + VERIFY_IS_APPROX( ge3 = tri.template triangularView<Upper>() * ge1, upTri * ge1); + VERIFY_IS_APPROX(rge3 = tri.template triangularView<Upper>() * ge1, upTri * ge1); + VERIFY_IS_APPROX( ge3 = (s1*tri.adjoint()).template triangularView<Upper>() * (s2*ge1), s1*loTri.adjoint() * (s2*ge1)); + VERIFY_IS_APPROX(rge3 = tri.adjoint().template triangularView<Upper>() * ge1, loTri.adjoint() * ge1); + VERIFY_IS_APPROX( ge3 = tri.adjoint().template triangularView<Lower>() * ge1, upTri.adjoint() * ge1); + VERIFY_IS_APPROX(rge3 = tri.adjoint().template triangularView<Lower>() * ge1, upTri.adjoint() * ge1); + VERIFY_IS_APPROX( ge3 = tri.template triangularView<Lower>() * ge2.adjoint(), loTri * ge2.adjoint()); + VERIFY_IS_APPROX(rge3 = tri.template triangularView<Lower>() * ge2.adjoint(), loTri * ge2.adjoint()); + VERIFY_IS_APPROX( ge3 = tri.template triangularView<Upper>() * ge2.adjoint(), upTri * ge2.adjoint()); + VERIFY_IS_APPROX(rge3 = tri.template triangularView<Upper>() * ge2.adjoint(), upTri * ge2.adjoint()); + VERIFY_IS_APPROX( ge3 = (s1*tri).adjoint().template triangularView<Upper>() * ge2.adjoint(), ei_conj(s1) * loTri.adjoint() * ge2.adjoint()); + VERIFY_IS_APPROX(rge3 = tri.adjoint().template triangularView<Upper>() * ge2.adjoint(), loTri.adjoint() * ge2.adjoint()); + VERIFY_IS_APPROX( ge3 = tri.adjoint().template triangularView<Lower>() * ge2.adjoint(), upTri.adjoint() * ge2.adjoint()); + VERIFY_IS_APPROX(rge3 = tri.adjoint().template triangularView<Lower>() * ge2.adjoint(), upTri.adjoint() * ge2.adjoint()); } void test_product_trmm() diff --git a/test/product_trmv.cpp b/test/product_trmv.cpp index 5016a5b1f..f0962557a 100644 --- a/test/product_trmv.cpp +++ b/test/product_trmv.cpp @@ -44,37 +44,37 @@ template<typename MatrixType> void trmv(const MatrixType& m) m1 = MatrixType::Random(rows, cols); // check with a column-major matrix - m3 = m1.template triangularView<Eigen::LowerTriangular>(); - VERIFY((m3 * v1).isApprox(m1.template triangularView<Eigen::LowerTriangular>() * v1, largerEps)); - m3 = m1.template triangularView<Eigen::UpperTriangular>(); - VERIFY((m3 * v1).isApprox(m1.template triangularView<Eigen::UpperTriangular>() * v1, largerEps)); - m3 = m1.template triangularView<Eigen::UnitLowerTriangular>(); - VERIFY((m3 * v1).isApprox(m1.template triangularView<Eigen::UnitLowerTriangular>() * v1, largerEps)); - m3 = m1.template triangularView<Eigen::UnitUpperTriangular>(); - VERIFY((m3 * v1).isApprox(m1.template triangularView<Eigen::UnitUpperTriangular>() * v1, largerEps)); + m3 = m1.template triangularView<Eigen::Lower>(); + VERIFY((m3 * v1).isApprox(m1.template triangularView<Eigen::Lower>() * v1, largerEps)); + m3 = m1.template triangularView<Eigen::Upper>(); + VERIFY((m3 * v1).isApprox(m1.template triangularView<Eigen::Upper>() * v1, largerEps)); + m3 = m1.template triangularView<Eigen::UnitLower>(); + VERIFY((m3 * v1).isApprox(m1.template triangularView<Eigen::UnitLower>() * v1, largerEps)); + m3 = m1.template triangularView<Eigen::UnitUpper>(); + VERIFY((m3 * v1).isApprox(m1.template triangularView<Eigen::UnitUpper>() * v1, largerEps)); // check conjugated and scalar multiple expressions (col-major) - m3 = m1.template triangularView<Eigen::LowerTriangular>(); - VERIFY(((s1*m3).conjugate() * v1).isApprox((s1*m1).conjugate().template triangularView<Eigen::LowerTriangular>() * v1, largerEps)); - m3 = m1.template triangularView<Eigen::UpperTriangular>(); - VERIFY((m3.conjugate() * v1.conjugate()).isApprox(m1.conjugate().template triangularView<Eigen::UpperTriangular>() * v1.conjugate(), largerEps)); + m3 = m1.template triangularView<Eigen::Lower>(); + VERIFY(((s1*m3).conjugate() * v1).isApprox((s1*m1).conjugate().template triangularView<Eigen::Lower>() * v1, largerEps)); + m3 = m1.template triangularView<Eigen::Upper>(); + VERIFY((m3.conjugate() * v1.conjugate()).isApprox(m1.conjugate().template triangularView<Eigen::Upper>() * v1.conjugate(), largerEps)); // check with a row-major matrix - m3 = m1.template triangularView<Eigen::UpperTriangular>(); - VERIFY((m3.transpose() * v1).isApprox(m1.transpose().template triangularView<Eigen::LowerTriangular>() * v1, largerEps)); - m3 = m1.template triangularView<Eigen::LowerTriangular>(); - VERIFY((m3.transpose() * v1).isApprox(m1.transpose().template triangularView<Eigen::UpperTriangular>() * v1, largerEps)); - m3 = m1.template triangularView<Eigen::UnitUpperTriangular>(); - VERIFY((m3.transpose() * v1).isApprox(m1.transpose().template triangularView<Eigen::UnitLowerTriangular>() * v1, largerEps)); - m3 = m1.template triangularView<Eigen::UnitLowerTriangular>(); - VERIFY((m3.transpose() * v1).isApprox(m1.transpose().template triangularView<Eigen::UnitUpperTriangular>() * v1, largerEps)); + m3 = m1.template triangularView<Eigen::Upper>(); + VERIFY((m3.transpose() * v1).isApprox(m1.transpose().template triangularView<Eigen::Lower>() * v1, largerEps)); + m3 = m1.template triangularView<Eigen::Lower>(); + VERIFY((m3.transpose() * v1).isApprox(m1.transpose().template triangularView<Eigen::Upper>() * v1, largerEps)); + m3 = m1.template triangularView<Eigen::UnitUpper>(); + VERIFY((m3.transpose() * v1).isApprox(m1.transpose().template triangularView<Eigen::UnitLower>() * v1, largerEps)); + m3 = m1.template triangularView<Eigen::UnitLower>(); + VERIFY((m3.transpose() * v1).isApprox(m1.transpose().template triangularView<Eigen::UnitUpper>() * v1, largerEps)); // check conjugated and scalar multiple expressions (row-major) - m3 = m1.template triangularView<Eigen::UpperTriangular>(); - VERIFY((m3.adjoint() * v1).isApprox(m1.adjoint().template triangularView<Eigen::LowerTriangular>() * v1, largerEps)); - m3 = m1.template triangularView<Eigen::LowerTriangular>(); - VERIFY((m3.adjoint() * (s1*v1.conjugate())).isApprox(m1.adjoint().template triangularView<Eigen::UpperTriangular>() * (s1*v1.conjugate()), largerEps)); - m3 = m1.template triangularView<Eigen::UnitUpperTriangular>(); + m3 = m1.template triangularView<Eigen::Upper>(); + VERIFY((m3.adjoint() * v1).isApprox(m1.adjoint().template triangularView<Eigen::Lower>() * v1, largerEps)); + m3 = m1.template triangularView<Eigen::Lower>(); + VERIFY((m3.adjoint() * (s1*v1.conjugate())).isApprox(m1.adjoint().template triangularView<Eigen::Upper>() * (s1*v1.conjugate()), largerEps)); + m3 = m1.template triangularView<Eigen::UnitUpper>(); // TODO check with sub-matrices } diff --git a/test/product_trsolve.cpp b/test/product_trsolve.cpp index 4477a29d1..6e916230e 100644 --- a/test/product_trsolve.cpp +++ b/test/product_trsolve.cpp @@ -49,28 +49,28 @@ template<typename Scalar,int Size, int Cols> void trsolve(int size=Size,int cols cmLhs.setRandom(); cmLhs *= static_cast<RealScalar>(0.1); cmLhs.diagonal().array() += static_cast<RealScalar>(1); rmLhs.setRandom(); rmLhs *= static_cast<RealScalar>(0.1); rmLhs.diagonal().array() += static_cast<RealScalar>(1); - VERIFY_TRSM(cmLhs.conjugate().template triangularView<LowerTriangular>(), cmRhs); - VERIFY_TRSM(cmLhs .template triangularView<UpperTriangular>(), cmRhs); - VERIFY_TRSM(cmLhs .template triangularView<LowerTriangular>(), rmRhs); - VERIFY_TRSM(cmLhs.conjugate().template triangularView<UpperTriangular>(), rmRhs); + VERIFY_TRSM(cmLhs.conjugate().template triangularView<Lower>(), cmRhs); + VERIFY_TRSM(cmLhs .template triangularView<Upper>(), cmRhs); + VERIFY_TRSM(cmLhs .template triangularView<Lower>(), rmRhs); + VERIFY_TRSM(cmLhs.conjugate().template triangularView<Upper>(), rmRhs); - VERIFY_TRSM(cmLhs.conjugate().template triangularView<UnitLowerTriangular>(), cmRhs); - VERIFY_TRSM(cmLhs .template triangularView<UnitUpperTriangular>(), rmRhs); + VERIFY_TRSM(cmLhs.conjugate().template triangularView<UnitLower>(), cmRhs); + VERIFY_TRSM(cmLhs .template triangularView<UnitUpper>(), rmRhs); - VERIFY_TRSM(rmLhs .template triangularView<LowerTriangular>(), cmRhs); - VERIFY_TRSM(rmLhs.conjugate().template triangularView<UnitUpperTriangular>(), rmRhs); + VERIFY_TRSM(rmLhs .template triangularView<Lower>(), cmRhs); + VERIFY_TRSM(rmLhs.conjugate().template triangularView<UnitUpper>(), rmRhs); - VERIFY_TRSM_ONTHERIGHT(cmLhs.conjugate().template triangularView<LowerTriangular>(), cmRhs); - VERIFY_TRSM_ONTHERIGHT(cmLhs .template triangularView<UpperTriangular>(), cmRhs); - VERIFY_TRSM_ONTHERIGHT(cmLhs .template triangularView<LowerTriangular>(), rmRhs); - VERIFY_TRSM_ONTHERIGHT(cmLhs.conjugate().template triangularView<UpperTriangular>(), rmRhs); + VERIFY_TRSM_ONTHERIGHT(cmLhs.conjugate().template triangularView<Lower>(), cmRhs); + VERIFY_TRSM_ONTHERIGHT(cmLhs .template triangularView<Upper>(), cmRhs); + VERIFY_TRSM_ONTHERIGHT(cmLhs .template triangularView<Lower>(), rmRhs); + VERIFY_TRSM_ONTHERIGHT(cmLhs.conjugate().template triangularView<Upper>(), rmRhs); - VERIFY_TRSM_ONTHERIGHT(cmLhs.conjugate().template triangularView<UnitLowerTriangular>(), cmRhs); - VERIFY_TRSM_ONTHERIGHT(cmLhs .template triangularView<UnitUpperTriangular>(), rmRhs); + VERIFY_TRSM_ONTHERIGHT(cmLhs.conjugate().template triangularView<UnitLower>(), cmRhs); + VERIFY_TRSM_ONTHERIGHT(cmLhs .template triangularView<UnitUpper>(), rmRhs); - VERIFY_TRSM_ONTHERIGHT(rmLhs .template triangularView<LowerTriangular>(), cmRhs); - VERIFY_TRSM_ONTHERIGHT(rmLhs.conjugate().template triangularView<UnitUpperTriangular>(), rmRhs); + VERIFY_TRSM_ONTHERIGHT(rmLhs .template triangularView<Lower>(), cmRhs); + VERIFY_TRSM_ONTHERIGHT(rmLhs.conjugate().template triangularView<UnitUpper>(), rmRhs); } void test_product_trsolve() diff --git a/test/qr_colpivoting.cpp b/test/qr_colpivoting.cpp index 8b56cd296..16eb27c52 100644 --- a/test/qr_colpivoting.cpp +++ b/test/qr_colpivoting.cpp @@ -47,7 +47,7 @@ template<typename MatrixType> void qr() MatrixQType q = qr.householderQ(); VERIFY_IS_UNITARY(q); - MatrixType r = qr.matrixQR().template triangularView<UpperTriangular>(); + MatrixType r = qr.matrixQR().template triangularView<Upper>(); MatrixType c = q * r * qr.colsPermutation().inverse(); VERIFY_IS_APPROX(m1, c); @@ -72,7 +72,7 @@ template<typename MatrixType, int Cols2> void qr_fixedsize() VERIFY(!qr.isInvertible()); VERIFY(!qr.isSurjective()); - Matrix<Scalar,Rows,Cols> r = qr.matrixQR().template triangularView<UpperTriangular>(); + Matrix<Scalar,Rows,Cols> r = qr.matrixQR().template triangularView<Upper>(); Matrix<Scalar,Rows,Cols> c = qr.householderQ() * r * qr.colsPermutation().inverse(); VERIFY_IS_APPROX(m1, c); diff --git a/test/sparse_product.cpp b/test/sparse_product.cpp index 00c2cdf74..04e089784 100644 --- a/test/sparse_product.cpp +++ b/test/sparse_product.cpp @@ -115,9 +115,9 @@ template<typename SparseMatrixType> void sparse_product(const SparseMatrixType& VERIFY_IS_APPROX(mS, refS); VERIFY_IS_APPROX(x=mS*b, refX=refS*b); - VERIFY_IS_APPROX(x=mUp.template selfadjointView<UpperTriangular>()*b, refX=refS*b); - VERIFY_IS_APPROX(x=mLo.template selfadjointView<LowerTriangular>()*b, refX=refS*b); - VERIFY_IS_APPROX(x=mS.template selfadjointView<UpperTriangular|LowerTriangular>()*b, refX=refS*b); + VERIFY_IS_APPROX(x=mUp.template selfadjointView<Upper>()*b, refX=refS*b); + VERIFY_IS_APPROX(x=mLo.template selfadjointView<Lower>()*b, refX=refS*b); + VERIFY_IS_APPROX(x=mS.template selfadjointView<Upper|Lower>()*b, refX=refS*b); } } diff --git a/test/sparse_solvers.cpp b/test/sparse_solvers.cpp index 059747c3d..fab2ab56e 100644 --- a/test/sparse_solvers.cpp +++ b/test/sparse_solvers.cpp @@ -65,23 +65,23 @@ template<typename Scalar> void sparse_solvers(int rows, int cols) // lower - dense initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords); - VERIFY_IS_APPROX(refMat2.template triangularView<LowerTriangular>().solve(vec2), - m2.template triangularView<LowerTriangular>().solve(vec3)); + VERIFY_IS_APPROX(refMat2.template triangularView<Lower>().solve(vec2), + m2.template triangularView<Lower>().solve(vec3)); // upper - dense initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, &zeroCoords, &nonzeroCoords); - VERIFY_IS_APPROX(refMat2.template triangularView<UpperTriangular>().solve(vec2), - m2.template triangularView<UpperTriangular>().solve(vec3)); + VERIFY_IS_APPROX(refMat2.template triangularView<Upper>().solve(vec2), + m2.template triangularView<Upper>().solve(vec3)); // lower - transpose initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords); - VERIFY_IS_APPROX(refMat2.transpose().template triangularView<UpperTriangular>().solve(vec2), - m2.transpose().template triangularView<UpperTriangular>().solve(vec3)); + VERIFY_IS_APPROX(refMat2.transpose().template triangularView<Upper>().solve(vec2), + m2.transpose().template triangularView<Upper>().solve(vec3)); // upper - transpose initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, &zeroCoords, &nonzeroCoords); - VERIFY_IS_APPROX(refMat2.transpose().template triangularView<LowerTriangular>().solve(vec2), - m2.transpose().template triangularView<LowerTriangular>().solve(vec3)); + VERIFY_IS_APPROX(refMat2.transpose().template triangularView<Lower>().solve(vec2), + m2.transpose().template triangularView<Lower>().solve(vec3)); SparseMatrix<Scalar> matB(rows, rows); DenseMatrix refMatB = DenseMatrix::Zero(rows, rows); @@ -89,21 +89,21 @@ template<typename Scalar> void sparse_solvers(int rows, int cols) // lower - sparse initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular); initSparse<Scalar>(density, refMatB, matB); - refMat2.template triangularView<LowerTriangular>().solveInPlace(refMatB); - m2.template triangularView<LowerTriangular>().solveInPlace(matB); + refMat2.template triangularView<Lower>().solveInPlace(refMatB); + m2.template triangularView<Lower>().solveInPlace(matB); VERIFY_IS_APPROX(matB.toDense(), refMatB); // upper - sparse initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular); initSparse<Scalar>(density, refMatB, matB); - refMat2.template triangularView<UpperTriangular>().solveInPlace(refMatB); - m2.template triangularView<UpperTriangular>().solveInPlace(matB); + refMat2.template triangularView<Upper>().solveInPlace(refMatB); + m2.template triangularView<Upper>().solveInPlace(matB); VERIFY_IS_APPROX(matB, refMatB); // test deprecated API initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords); - VERIFY_IS_APPROX(refMat2.template triangularView<LowerTriangular>().solve(vec2), - m2.template triangularView<LowerTriangular>().solve(vec3)); + VERIFY_IS_APPROX(refMat2.template triangularView<Lower>().solve(vec2), + m2.template triangularView<Lower>().solve(vec3)); } // test LLT @@ -118,29 +118,28 @@ template<typename Scalar> void sparse_solvers(int rows, int cols) initSPD(density, refMat2, m2); refX = refMat2.llt().solve(b); - typedef SparseMatrix<Scalar,LowerTriangular|SelfAdjoint> SparseSelfAdjointMatrix; if (!NumTraits<Scalar>::IsComplex) { x = b; - SparseLLT<SparseSelfAdjointMatrix> (m2).solveInPlace(x); + SparseLLT<SparseMatrix<Scalar> > (m2).solveInPlace(x); VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: default"); } #ifdef EIGEN_CHOLMOD_SUPPORT x = b; - SparseLLT<SparseSelfAdjointMatrix,Cholmod>(m2).solveInPlace(x); + SparseLLT<SparseMatrix<Scalar> ,Cholmod>(m2).solveInPlace(x); VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: cholmod"); #endif #ifdef EIGEN_TAUCS_SUPPORT x = b; - SparseLLT<SparseSelfAdjointMatrix,Taucs>(m2,IncompleteFactorization).solveInPlace(x); + SparseLLT<SparseMatrix<Scalar> ,Taucs>(m2,IncompleteFactorization).solveInPlace(x); VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (IncompleteFactorization)"); // TODO fix TAUCS with complexes x = b; - SparseLLT<SparseSelfAdjointMatrix,Taucs>(m2,SupernodalMultifrontal).solveInPlace(x); + SparseLLT<SparseMatrix<Scalar> ,Taucs>(m2,SupernodalMultifrontal).solveInPlace(x); VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (SupernodalMultifrontal)"); x = b; - SparseLLT<SparseSelfAdjointMatrix,Taucs>(m2,SupernodalLeftLooking).solveInPlace(x); + SparseLLT<SparseMatrix<Scalar> ,Taucs>(m2,SupernodalLeftLooking).solveInPlace(x); VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (SupernodalLeftLooking)"); #endif } @@ -161,7 +160,7 @@ template<typename Scalar> void sparse_solvers(int rows, int cols) refMat2.diagonal() *= 0.5; refX = refMat2.llt().solve(b); // FIXME use LLT to compute the reference because LDLT seems to fail with large matrices - typedef SparseMatrix<Scalar,UpperTriangular|SelfAdjoint> SparseSelfAdjointMatrix; + typedef SparseMatrix<Scalar,Upper|SelfAdjoint> SparseSelfAdjointMatrix; x = b; SparseLDLT<SparseSelfAdjointMatrix> ldlt(m2); if (ldlt.succeeded()) diff --git a/test/triangular.cpp b/test/triangular.cpp index 0de9f5841..2903e247c 100644 --- a/test/triangular.cpp +++ b/test/triangular.cpp @@ -53,14 +53,14 @@ template<typename MatrixType> void triangular_square(const MatrixType& m) v2 = VectorType::Random(rows), vzero = VectorType::Zero(rows); - MatrixType m1up = m1.template triangularView<UpperTriangular>(); - MatrixType m2up = m2.template triangularView<UpperTriangular>(); + MatrixType m1up = m1.template triangularView<Upper>(); + MatrixType m2up = m2.template triangularView<Upper>(); if (rows*cols>1) { - VERIFY(m1up.isUpperTriangular()); - VERIFY(m2up.transpose().isLowerTriangular()); - VERIFY(!m2.isLowerTriangular()); + VERIFY(m1up.isUpper()); + VERIFY(m2up.transpose().isLower()); + VERIFY(!m2.isLower()); } // VERIFY_IS_APPROX(m1up.transpose() * m2, m1.upper().transpose().lower() * m2); @@ -68,20 +68,20 @@ template<typename MatrixType> void triangular_square(const MatrixType& m) // test overloaded operator+= r1.setZero(); r2.setZero(); - r1.template triangularView<UpperTriangular>() += m1; + r1.template triangularView<Upper>() += m1; r2 += m1up; VERIFY_IS_APPROX(r1,r2); // test overloaded operator= m1.setZero(); - m1.template triangularView<UpperTriangular>() = m2.transpose() + m2; + m1.template triangularView<Upper>() = m2.transpose() + m2; m3 = m2.transpose() + m2; - VERIFY_IS_APPROX(m3.template triangularView<LowerTriangular>().transpose().toDenseMatrix(), m1); + VERIFY_IS_APPROX(m3.template triangularView<Lower>().transpose().toDenseMatrix(), m1); // test overloaded operator= m1.setZero(); - m1.template triangularView<LowerTriangular>() = m2.transpose() + m2; - VERIFY_IS_APPROX(m3.template triangularView<LowerTriangular>().toDenseMatrix(), m1); + m1.template triangularView<Lower>() = m2.transpose() + m2; + VERIFY_IS_APPROX(m3.template triangularView<Lower>().toDenseMatrix(), m1); m1 = MatrixType::Random(rows, cols); for (int i=0; i<rows; ++i) @@ -89,49 +89,49 @@ template<typename MatrixType> void triangular_square(const MatrixType& m) Transpose<MatrixType> trm4(m4); // test back and forward subsitution with a vector as the rhs - m3 = m1.template triangularView<UpperTriangular>(); - VERIFY(v2.isApprox(m3.adjoint() * (m1.adjoint().template triangularView<LowerTriangular>().solve(v2)), largerEps)); - m3 = m1.template triangularView<LowerTriangular>(); - VERIFY(v2.isApprox(m3.transpose() * (m1.transpose().template triangularView<UpperTriangular>().solve(v2)), largerEps)); - m3 = m1.template triangularView<UpperTriangular>(); - VERIFY(v2.isApprox(m3 * (m1.template triangularView<UpperTriangular>().solve(v2)), largerEps)); - m3 = m1.template triangularView<LowerTriangular>(); - VERIFY(v2.isApprox(m3.conjugate() * (m1.conjugate().template triangularView<LowerTriangular>().solve(v2)), largerEps)); + m3 = m1.template triangularView<Upper>(); + VERIFY(v2.isApprox(m3.adjoint() * (m1.adjoint().template triangularView<Lower>().solve(v2)), largerEps)); + m3 = m1.template triangularView<Lower>(); + VERIFY(v2.isApprox(m3.transpose() * (m1.transpose().template triangularView<Upper>().solve(v2)), largerEps)); + m3 = m1.template triangularView<Upper>(); + VERIFY(v2.isApprox(m3 * (m1.template triangularView<Upper>().solve(v2)), largerEps)); + m3 = m1.template triangularView<Lower>(); + VERIFY(v2.isApprox(m3.conjugate() * (m1.conjugate().template triangularView<Lower>().solve(v2)), largerEps)); // test back and forward subsitution with a matrix as the rhs - m3 = m1.template triangularView<UpperTriangular>(); - VERIFY(m2.isApprox(m3.adjoint() * (m1.adjoint().template triangularView<LowerTriangular>().solve(m2)), largerEps)); - m3 = m1.template triangularView<LowerTriangular>(); - VERIFY(m2.isApprox(m3.transpose() * (m1.transpose().template triangularView<UpperTriangular>().solve(m2)), largerEps)); - m3 = m1.template triangularView<UpperTriangular>(); - VERIFY(m2.isApprox(m3 * (m1.template triangularView<UpperTriangular>().solve(m2)), largerEps)); - m3 = m1.template triangularView<LowerTriangular>(); - VERIFY(m2.isApprox(m3.conjugate() * (m1.conjugate().template triangularView<LowerTriangular>().solve(m2)), largerEps)); + m3 = m1.template triangularView<Upper>(); + VERIFY(m2.isApprox(m3.adjoint() * (m1.adjoint().template triangularView<Lower>().solve(m2)), largerEps)); + m3 = m1.template triangularView<Lower>(); + VERIFY(m2.isApprox(m3.transpose() * (m1.transpose().template triangularView<Upper>().solve(m2)), largerEps)); + m3 = m1.template triangularView<Upper>(); + VERIFY(m2.isApprox(m3 * (m1.template triangularView<Upper>().solve(m2)), largerEps)); + m3 = m1.template triangularView<Lower>(); + VERIFY(m2.isApprox(m3.conjugate() * (m1.conjugate().template triangularView<Lower>().solve(m2)), largerEps)); // check M * inv(L) using in place API m4 = m3; - m3.transpose().template triangularView<Eigen::UpperTriangular>().solveInPlace(trm4); + m3.transpose().template triangularView<Eigen::Upper>().solveInPlace(trm4); VERIFY(m4.cwiseAbs().isIdentity(test_precision<RealScalar>())); // check M * inv(U) using in place API - m3 = m1.template triangularView<UpperTriangular>(); + m3 = m1.template triangularView<Upper>(); m4 = m3; - m3.transpose().template triangularView<Eigen::LowerTriangular>().solveInPlace(trm4); + m3.transpose().template triangularView<Eigen::Lower>().solveInPlace(trm4); VERIFY(m4.cwiseAbs().isIdentity(test_precision<RealScalar>())); // check solve with unit diagonal - m3 = m1.template triangularView<UnitUpperTriangular>(); - VERIFY(m2.isApprox(m3 * (m1.template triangularView<UnitUpperTriangular>().solve(m2)), largerEps)); + m3 = m1.template triangularView<UnitUpper>(); + VERIFY(m2.isApprox(m3 * (m1.template triangularView<UnitUpper>().solve(m2)), largerEps)); -// VERIFY(( m1.template triangularView<UpperTriangular>() -// * m2.template triangularView<UpperTriangular>()).isUpperTriangular()); +// VERIFY(( m1.template triangularView<Upper>() +// * m2.template triangularView<Upper>()).isUpper()); // test swap m1.setOnes(); m2.setZero(); - m2.template triangularView<UpperTriangular>().swap(m1); + m2.template triangularView<Upper>().swap(m1); m3.setZero(); - m3.template triangularView<UpperTriangular>().setOnes(); + m3.template triangularView<Upper>().setOnes(); VERIFY_IS_APPROX(m2,m3); } @@ -165,69 +165,69 @@ template<typename MatrixType> void triangular_rect(const MatrixType& m) v2 = VectorType::Random(rows), vzero = VectorType::Zero(rows); - MatrixType m1up = m1.template triangularView<UpperTriangular>(); - MatrixType m2up = m2.template triangularView<UpperTriangular>(); + MatrixType m1up = m1.template triangularView<Upper>(); + MatrixType m2up = m2.template triangularView<Upper>(); if (rows*cols>1) { - VERIFY(m1up.isUpperTriangular()); - VERIFY(m2up.transpose().isLowerTriangular()); - VERIFY(!m2.isLowerTriangular()); + VERIFY(m1up.isUpper()); + VERIFY(m2up.transpose().isLower()); + VERIFY(!m2.isLower()); } // test overloaded operator+= r1.setZero(); r2.setZero(); - r1.template triangularView<UpperTriangular>() += m1; + r1.template triangularView<Upper>() += m1; r2 += m1up; VERIFY_IS_APPROX(r1,r2); // test overloaded operator= m1.setZero(); - m1.template triangularView<UpperTriangular>() = 3 * m2; + m1.template triangularView<Upper>() = 3 * m2; m3 = 3 * m2; - VERIFY_IS_APPROX(m3.template triangularView<UpperTriangular>().toDenseMatrix(), m1); + VERIFY_IS_APPROX(m3.template triangularView<Upper>().toDenseMatrix(), m1); m1.setZero(); - m1.template triangularView<LowerTriangular>() = 3 * m2; - VERIFY_IS_APPROX(m3.template triangularView<LowerTriangular>().toDenseMatrix(), m1); + m1.template triangularView<Lower>() = 3 * m2; + VERIFY_IS_APPROX(m3.template triangularView<Lower>().toDenseMatrix(), m1); m1.setZero(); - m1.template triangularView<StrictlyUpperTriangular>() = 3 * m2; - VERIFY_IS_APPROX(m3.template triangularView<StrictlyUpperTriangular>().toDenseMatrix(), m1); + m1.template triangularView<StrictlyUpper>() = 3 * m2; + VERIFY_IS_APPROX(m3.template triangularView<StrictlyUpper>().toDenseMatrix(), m1); m1.setZero(); - m1.template triangularView<StrictlyLowerTriangular>() = 3 * m2; - VERIFY_IS_APPROX(m3.template triangularView<StrictlyLowerTriangular>().toDenseMatrix(), m1); + m1.template triangularView<StrictlyLower>() = 3 * m2; + VERIFY_IS_APPROX(m3.template triangularView<StrictlyLower>().toDenseMatrix(), m1); m1.setRandom(); - m2 = m1.template triangularView<UpperTriangular>(); - VERIFY(m2.isUpperTriangular()); - VERIFY(!m2.isLowerTriangular()); - m2 = m1.template triangularView<StrictlyUpperTriangular>(); - VERIFY(m2.isUpperTriangular()); + m2 = m1.template triangularView<Upper>(); + VERIFY(m2.isUpper()); + VERIFY(!m2.isLower()); + m2 = m1.template triangularView<StrictlyUpper>(); + VERIFY(m2.isUpper()); VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1))); - m2 = m1.template triangularView<UnitUpperTriangular>(); - VERIFY(m2.isUpperTriangular()); + m2 = m1.template triangularView<UnitUpper>(); + VERIFY(m2.isUpper()); m2.diagonal().array() -= Scalar(1); VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1))); - m2 = m1.template triangularView<LowerTriangular>(); - VERIFY(m2.isLowerTriangular()); - VERIFY(!m2.isUpperTriangular()); - m2 = m1.template triangularView<StrictlyLowerTriangular>(); - VERIFY(m2.isLowerTriangular()); + m2 = m1.template triangularView<Lower>(); + VERIFY(m2.isLower()); + VERIFY(!m2.isUpper()); + m2 = m1.template triangularView<StrictlyLower>(); + VERIFY(m2.isLower()); VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1))); - m2 = m1.template triangularView<UnitLowerTriangular>(); - VERIFY(m2.isLowerTriangular()); + m2 = m1.template triangularView<UnitLower>(); + VERIFY(m2.isLower()); m2.diagonal().array() -= Scalar(1); VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1))); // test swap m1.setOnes(); m2.setZero(); - m2.template triangularView<UpperTriangular>().swap(m1); + m2.template triangularView<Upper>().swap(m1); m3.setZero(); - m3.template triangularView<UpperTriangular>().setOnes(); + m3.template triangularView<Upper>().setOnes(); VERIFY_IS_APPROX(m2,m3); } |