From 9e00d945439d801d3f4e33ed1ce57545e3310723 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Sat, 20 Dec 2008 13:36:12 +0000 Subject: * the Upper->UpperTriangular change * finally get ei_add_test right --- Eigen/src/Cholesky/Cholesky.h | 4 +- Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h | 6 +-- Eigen/src/Cholesky/LDLT.h | 4 +- Eigen/src/Cholesky/LLT.h | 4 +- Eigen/src/Core/MatrixBase.h | 4 +- Eigen/src/Core/Part.h | 58 ++++++++++---------- Eigen/src/Core/SolveTriangular.h | 74 +++++++++++++------------- Eigen/src/Core/Transpose.h | 4 +- Eigen/src/Core/util/Constants.h | 14 ++--- Eigen/src/LU/LU.h | 10 ++-- Eigen/src/QR/HessenbergDecomposition.h | 2 +- Eigen/src/QR/QR.h | 4 +- Eigen/src/QR/SelfAdjointEigenSolver.h | 8 +-- Eigen/src/QR/Tridiagonalization.h | 12 ++--- Eigen/src/Sparse/CholmodSupport.h | 4 +- Eigen/src/Sparse/SparseLDLT.h | 2 +- Eigen/src/Sparse/SparseLLT.h | 2 +- Eigen/src/Sparse/SparseLU.h | 2 +- Eigen/src/Sparse/SuperLUSupport.h | 8 +-- Eigen/src/Sparse/TaucsSupport.h | 6 +-- Eigen/src/Sparse/TriangularSolver.h | 8 +-- Eigen/src/Sparse/UmfPackSupport.h | 4 +- bench/btl/libs/eigen2/eigen2_interface.hh | 4 +- bench/sparse_cholesky.cpp | 4 +- bench/sparse_trisolver.cpp | 8 +-- doc/QuickStartGuide.dox | 20 +++---- doc/snippets/MatrixBase_extract.cpp | 6 +-- doc/snippets/MatrixBase_marked.cpp | 6 +-- doc/snippets/MatrixBase_part.cpp | 2 +- doc/snippets/class_LU_2.cpp | 2 +- test/CMakeLists.txt | 6 ++- test/sparse_solvers.cpp | 12 ++--- test/triangular.cpp | 46 ++++++++-------- 33 files changed, 181 insertions(+), 179 deletions(-) diff --git a/Eigen/src/Cholesky/Cholesky.h b/Eigen/src/Cholesky/Cholesky.h index bbf24fab3..31451a6b0 100644 --- a/Eigen/src/Cholesky/Cholesky.h +++ b/Eigen/src/Cholesky/Cholesky.h @@ -52,7 +52,7 @@ template class Cholesky } /** \deprecated */ - inline Part matrixL(void) const { return m_matrix; } + inline Part matrixL(void) const { return m_matrix; } /** \deprecated */ inline bool isPositiveDefinite(void) const { return m_isPositiveDefinite; } @@ -148,7 +148,7 @@ bool Cholesky::solveInPlace(MatrixBase &bAndX) const if (!m_isPositiveDefinite) return false; matrixL().solveTriangularInPlace(bAndX); - m_matrix.adjoint().template part().solveTriangularInPlace(bAndX); + m_matrix.adjoint().template part().solveTriangularInPlace(bAndX); return true; } diff --git a/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h b/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h index 0f3e9e07f..641f9be01 100644 --- a/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h +++ b/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h @@ -46,7 +46,7 @@ template class CholeskyWithoutSquareRoot } /** \returns the lower triangular matrix L */ - inline Part matrixL(void) const { return m_matrix; } + inline Part matrixL(void) const { return m_matrix; } /** \returns the coefficients of the diagonal matrix D */ inline DiagonalCoeffs vectorD(void) const { return m_matrix.diagonal(); } @@ -137,7 +137,7 @@ typename Derived::Eval CholeskyWithoutSquareRoot::solve(const Matrix const int size = m_matrix.rows(); ei_assert(size==b.rows()); - return m_matrix.adjoint().template part() + return m_matrix.adjoint().template part() .solveTriangular( ( m_matrix.cwise().inverse().template part() * matrixL().solveTriangular(b)) @@ -167,7 +167,7 @@ bool CholeskyWithoutSquareRoot::solveInPlace(MatrixBase &bA return false; matrixL().solveTriangularInPlace(bAndX); bAndX = (m_matrix.cwise().inverse().template part() * bAndX).lazy(); - m_matrix.adjoint().template part().solveTriangularInPlace(bAndX); + m_matrix.adjoint().template part().solveTriangularInPlace(bAndX); return true; } diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index a275e093f..4729c3621 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -60,7 +60,7 @@ template class LDLT } /** \returns the lower triangular matrix L */ - inline Part matrixL(void) const { return m_matrix; } + inline Part matrixL(void) const { return m_matrix; } /** \returns the coefficients of the diagonal matrix D */ inline DiagonalCoeffs vectorD(void) const { return m_matrix.diagonal(); } @@ -181,7 +181,7 @@ bool LDLT::solveInPlace(MatrixBase &bAndX) const return false; matrixL().solveTriangularInPlace(bAndX); bAndX = (m_matrix.cwise().inverse().template part() * bAndX).lazy(); - m_matrix.adjoint().template part().solveTriangularInPlace(bAndX); + m_matrix.adjoint().template part().solveTriangularInPlace(bAndX); return true; } diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h index 6e412e20c..c4e7fdffd 100644 --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h @@ -67,7 +67,7 @@ template class LLT } /** \returns the lower triangular matrix L */ - inline Part matrixL(void) const { return m_matrix; } + inline Part matrixL(void) const { return m_matrix; } /** \returns true if the matrix is positive definite */ inline bool isPositiveDefinite(void) const { return m_isPositiveDefinite; } @@ -169,7 +169,7 @@ bool LLT::solveInPlace(MatrixBase &bAndX) const if (!m_isPositiveDefinite) return false; matrixL().solveTriangularInPlace(bAndX); - m_matrix.adjoint().template part().solveTriangularInPlace(bAndX); + m_matrix.adjoint().template part().solveTriangularInPlace(bAndX); return true; } diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index a9a08dc28..42d3bece3 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -470,8 +470,8 @@ template class MatrixBase bool isIdentity(RealScalar prec = precision()) const; bool isDiagonal(RealScalar prec = precision()) const; - bool isUpper(RealScalar prec = precision()) const; - bool isLower(RealScalar prec = precision()) const; + bool isUpperTriangular(RealScalar prec = precision()) const; + bool isLowerTriangular(RealScalar prec = precision()) const; template bool isOrthogonal(const MatrixBase& other, diff --git a/Eigen/src/Core/Part.h b/Eigen/src/Core/Part.h index b1a0b3be9..883b92a98 100644 --- a/Eigen/src/Core/Part.h +++ b/Eigen/src/Core/Part.h @@ -31,8 +31,8 @@ * \brief Expression of a triangular matrix extracted from a given 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 Upper, StrictlyUpper, - * UnitUpper, Lower, StrictlyLower, UnitLower. This is in fact a bit field; it must have either + * \param Mode the kind of triangular matrix expression to construct. Can be UpperTriangular, StrictlyUpperTriangular, + * UnitUpperTriangular, LowerTriangular, StrictlyLowerTriangular, UnitLowerTriangular. This is in fact a bit field; it must have either * UpperTriangularBit or LowerTriangularBit, and additionnaly it may have either ZeroDiagBit or * UnitDiagBit. * @@ -104,10 +104,10 @@ template class Part { EIGEN_STATIC_ASSERT(!(Flags & UnitDiagBit), WRITING_TO_TRIANGULAR_PART_WITH_UNIT_DIAGONAL_IS_NOT_SUPPORTED) EIGEN_STATIC_ASSERT(!(Flags & SelfAdjointBit), COEFFICIENT_WRITE_ACCESS_TO_SELFADJOINT_NOT_SUPPORTED) - ei_assert( (Mode==Upper && col>=row) - || (Mode==Lower && col<=row) - || (Mode==StrictlyUpper && col>row) - || (Mode==StrictlyLower && col=row) + || (Mode==LowerTriangular && col<=row) + || (Mode==StrictlyUpperTriangular && col>row) + || (Mode==StrictlyLowerTriangular && col class Part /** \returns an expression of a triangular matrix extracted from the current matrix * - * The parameter \a Mode can have the following values: \c Upper, \c StrictlyUpper, \c UnitUpper, - * \c Lower, \c StrictlyLower, \c UnitLower. + * The parameter \a Mode can have the following values: \c UpperTriangular, \c StrictlyUpperTriangular, \c UnitUpperTriangular, + * \c LowerTriangular, \c StrictlyLowerTriangular, \c UnitLowerTriangular. * * \addexample PartExample \label How to extract a triangular part of an arbitrary matrix * @@ -187,11 +187,11 @@ struct ei_part_assignment_impl } else { - ei_assert(Mode == Upper || Mode == Lower || Mode == StrictlyUpper || Mode == StrictlyLower); - if((Mode == Upper && row <= col) - || (Mode == Lower && row >= col) - || (Mode == StrictlyUpper && row < col) - || (Mode == StrictlyLower && row > col)) + ei_assert(Mode == UpperTriangular || Mode == LowerTriangular || Mode == StrictlyUpperTriangular || Mode == StrictlyLowerTriangular); + if((Mode == UpperTriangular && row <= col) + || (Mode == LowerTriangular && row >= col) + || (Mode == StrictlyUpperTriangular && row < col) + || (Mode == StrictlyLowerTriangular && row > col)) dst.copyCoeff(row, col, src); } } @@ -215,7 +215,7 @@ struct ei_part_assignment_impl }; template -struct ei_part_assignment_impl +struct ei_part_assignment_impl { inline static void run(Derived1 &dst, const Derived2 &src) { @@ -226,7 +226,7 @@ struct ei_part_assignment_impl }; template -struct ei_part_assignment_impl +struct ei_part_assignment_impl { inline static void run(Derived1 &dst, const Derived2 &src) { @@ -237,7 +237,7 @@ struct ei_part_assignment_impl }; template -struct ei_part_assignment_impl +struct ei_part_assignment_impl { inline static void run(Derived1 &dst, const Derived2 &src) { @@ -247,7 +247,7 @@ struct ei_part_assignment_impl } }; template -struct ei_part_assignment_impl +struct ei_part_assignment_impl { inline static void run(Derived1 &dst, const Derived2 &src) { @@ -285,8 +285,8 @@ void Part::lazyAssign(const Other& other) /** \returns a lvalue pseudo-expression allowing to perform special operations on \c *this. * - * The \a Mode parameter can have the following values: \c Upper, \c StrictlyUpper, \c Lower, - * \c StrictlyLower, \c SelfAdjoint. + * The \a Mode parameter can have the following values: \c UpperTriangular, \c StrictlyUpperTriangular, \c LowerTriangular, + * \c StrictlyLowerTriangular, \c SelfAdjoint. * * \addexample PartExample \label How to write to a triangular part of a matrix * @@ -305,44 +305,44 @@ inline Part MatrixBase::part() /** \returns true if *this is approximately equal to an upper triangular matrix, * within the precision given by \a prec. * - * \sa isLower(), extract(), part(), marked() + * \sa isLowerTriangular(), extract(), part(), marked() */ template -bool MatrixBase::isUpper(RealScalar prec) const +bool MatrixBase::isUpperTriangular(RealScalar prec) const { if(cols() != rows()) return false; - RealScalar maxAbsOnUpperPart = static_cast(-1); + RealScalar maxAbsOnUpperTriangularPart = static_cast(-1); for(int j = 0; j < cols(); ++j) for(int i = 0; i <= j; ++i) { RealScalar absValue = ei_abs(coeff(i,j)); - if(absValue > maxAbsOnUpperPart) maxAbsOnUpperPart = absValue; + if(absValue > maxAbsOnUpperTriangularPart) maxAbsOnUpperTriangularPart = absValue; } for(int j = 0; j < cols()-1; ++j) for(int i = j+1; i < rows(); ++i) - if(!ei_isMuchSmallerThan(coeff(i, j), maxAbsOnUpperPart, prec)) return false; + if(!ei_isMuchSmallerThan(coeff(i, j), maxAbsOnUpperTriangularPart, prec)) return false; return true; } /** \returns true if *this is approximately equal to a lower triangular matrix, * within the precision given by \a prec. * - * \sa isUpper(), extract(), part(), marked() + * \sa isUpperTriangular(), extract(), part(), marked() */ template -bool MatrixBase::isLower(RealScalar prec) const +bool MatrixBase::isLowerTriangular(RealScalar prec) const { if(cols() != rows()) return false; - RealScalar maxAbsOnLowerPart = static_cast(-1); + RealScalar maxAbsOnLowerTriangularPart = static_cast(-1); for(int j = 0; j < cols(); ++j) for(int i = j; i < rows(); ++i) { RealScalar absValue = ei_abs(coeff(i,j)); - if(absValue > maxAbsOnLowerPart) maxAbsOnLowerPart = absValue; + if(absValue > maxAbsOnLowerTriangularPart) maxAbsOnLowerTriangularPart = absValue; } for(int j = 1; j < cols(); ++j) for(int i = 0; i < j; ++i) - if(!ei_isMuchSmallerThan(coeff(i, j), maxAbsOnLowerPart, prec)) return false; + if(!ei_isMuchSmallerThan(coeff(i, j), maxAbsOnLowerTriangularPart, prec)) return false; return true; } diff --git a/Eigen/src/Core/SolveTriangular.h b/Eigen/src/Core/SolveTriangular.h index b58dab01d..1c586b865 100644 --- a/Eigen/src/Core/SolveTriangular.h +++ b/Eigen/src/Core/SolveTriangular.h @@ -30,9 +30,9 @@ template struct ei_is_part::value ? -1 // this is to solve ambiguous specializations : int(Lhs::Flags) & (RowMajorBit|SparseBit) @@ -56,14 +56,14 @@ struct ei_solve_triangular_selector typedef typename Rhs::Scalar Scalar; static void run(const Lhs& lhs, Rhs& other) { - const bool IsLower = (UpLo==Lower); + const bool IsLowerTriangular = (UpLo==LowerTriangular); const int size = lhs.cols(); /* We perform the inverse product per block of 4 rows such that we perfectly match * our optimized matrix * vector product. blockyStart represents the number of rows * we have process first using the non-block version. */ int blockyStart = (std::max(size-5,0)/4)*4; - if (IsLower) + if (IsLowerTriangular) blockyStart = size - blockyStart; else blockyStart -= 1; @@ -72,15 +72,15 @@ struct ei_solve_triangular_selector // process first rows using the non block version if(!(Lhs::Flags & UnitDiagBit)) { - if (IsLower) + if (IsLowerTriangular) other.coeffRef(0,c) = other.coeff(0,c)/lhs.coeff(0, 0); else other.coeffRef(size-1,c) = other.coeff(size-1, c)/lhs.coeff(size-1, size-1); } - for(int i=(IsLower ? 1 : size-2); IsLower ? iblockyStart; i += (IsLower ? 1 : -1) ) + for(int i=(IsLowerTriangular ? 1 : size-2); IsLowerTriangular ? iblockyStart; i += (IsLowerTriangular ? 1 : -1) ) { Scalar tmp = other.coeff(i,c) - - (IsLower ? ((lhs.row(i).start(i)) * other.col(c).start(i)).coeff(0,0) + - (IsLowerTriangular ? ((lhs.row(i).start(i)) * other.col(c).start(i)).coeff(0,0) : ((lhs.row(i).end(size-i-1)) * other.col(c).end(size-i-1)).coeff(0,0)); if (Lhs::Flags & UnitDiagBit) other.coeffRef(i,c) = tmp; @@ -89,15 +89,15 @@ struct ei_solve_triangular_selector } // now let's process the remaining rows 4 at once - for(int i=blockyStart; IsLower ? i0; ) + for(int i=blockyStart; IsLowerTriangular ? i0; ) { int startBlock = i; - int endBlock = startBlock + (IsLower ? 4 : -4); + int endBlock = startBlock + (IsLowerTriangular ? 4 : -4); /* Process the i cols times 4 rows block, and keep the result in a temporary vector */ // FIXME use fixed size block but take care to small fixed size matrices... Matrix btmp(4); - if (IsLower) + if (IsLowerTriangular) btmp = lhs.block(startBlock,0,4,i) * other.col(c).start(i); else btmp = lhs.block(i-3,i+1,4,size-1-i) * other.col(c).end(size-1-i); @@ -106,21 +106,21 @@ struct ei_solve_triangular_selector * btmp stores the diagonal coefficients used to update the remaining part of the result. */ { - Scalar tmp = other.coeff(startBlock,c)-btmp.coeff(IsLower?0:3); + Scalar tmp = other.coeff(startBlock,c)-btmp.coeff(IsLowerTriangular?0:3); if (Lhs::Flags & UnitDiagBit) other.coeffRef(i,c) = tmp; else other.coeffRef(i,c) = tmp/lhs.coeff(i,i); } - i += IsLower ? 1 : -1; - for (;IsLower ? iendBlock; i += IsLower ? 1 : -1) + i += IsLowerTriangular ? 1 : -1; + for (;IsLowerTriangular ? iendBlock; i += IsLowerTriangular ? 1 : -1) { - int remainingSize = IsLower ? i-startBlock : startBlock-i; + int remainingSize = IsLowerTriangular ? i-startBlock : startBlock-i; Scalar tmp = other.coeff(i,c) - - btmp.coeff(IsLower ? remainingSize : 3-remainingSize) - - ( lhs.row(i).segment(IsLower ? startBlock : i+1, remainingSize) - * other.col(c).segment(IsLower ? startBlock : i+1, remainingSize)).coeff(0,0); + - btmp.coeff(IsLowerTriangular ? remainingSize : 3-remainingSize) + - ( lhs.row(i).segment(IsLowerTriangular ? startBlock : i+1, remainingSize) + * other.col(c).segment(IsLowerTriangular ? startBlock : i+1, remainingSize)).coeff(0,0); if (Lhs::Flags & UnitDiagBit) other.coeffRef(i,c) = tmp; @@ -133,10 +133,10 @@ struct ei_solve_triangular_selector }; // Implements the following configurations: -// - inv(Lower, ColMajor) * Column vector -// - inv(Lower,UnitDiag,ColMajor) * Column vector -// - inv(Upper, ColMajor) * Column vector -// - inv(Upper,UnitDiag,ColMajor) * Column vector +// - inv(LowerTriangular, ColMajor) * Column vector +// - inv(LowerTriangular,UnitDiag,ColMajor) * Column vector +// - inv(UpperTriangular, ColMajor) * Column vector +// - inv(UpperTriangular,UnitDiag,ColMajor) * Column vector template struct ei_solve_triangular_selector { @@ -146,7 +146,7 @@ struct ei_solve_triangular_selector static void run(const Lhs& lhs, Rhs& other) { - static const bool IsLower = (UpLo==Lower); + static const bool IsLowerTriangular = (UpLo==LowerTriangular); const int size = lhs.cols(); for(int c=0 ; c * we can process using the block version. */ int blockyEnd = (std::max(size-5,0)/4)*4; - if (!IsLower) + if (!IsLowerTriangular) blockyEnd = size-1 - blockyEnd; - for(int i=IsLower ? 0 : size-1; IsLower ? iblockyEnd;) + for(int i=IsLowerTriangular ? 0 : size-1; IsLowerTriangular ? iblockyEnd;) { /* Let's process the 4x4 sub-matrix as usual. * btmp stores the diagonal coefficients used to update the remaining part of the result. */ int startBlock = i; - int endBlock = startBlock + (IsLower ? 4 : -4); + int endBlock = startBlock + (IsLowerTriangular ? 4 : -4); Matrix btmp; - for (;IsLower ? iendBlock; - i += IsLower ? 1 : -1) + for (;IsLowerTriangular ? iendBlock; + i += IsLowerTriangular ? 1 : -1) { if(!(Lhs::Flags & UnitDiagBit)) other.coeffRef(i,c) /= lhs.coeff(i,i); - int remainingSize = IsLower ? endBlock-i-1 : i-endBlock-1; + int remainingSize = IsLowerTriangular ? endBlock-i-1 : i-endBlock-1; if (remainingSize>0) - other.col(c).segment((IsLower ? i : endBlock) + 1, remainingSize) -= + other.col(c).segment((IsLowerTriangular ? i : endBlock) + 1, remainingSize) -= other.coeffRef(i,c) - * Block(lhs, (IsLower ? i : endBlock) + 1, i, remainingSize, 1); - btmp.coeffRef(IsLower ? i-startBlock : remainingSize) = -other.coeffRef(i,c); + * Block(lhs, (IsLowerTriangular ? i : endBlock) + 1, i, remainingSize, 1); + btmp.coeffRef(IsLowerTriangular ? i-startBlock : remainingSize) = -other.coeffRef(i,c); } /* Now we can efficiently update the remaining part of the result as a matrix * vector product. @@ -187,11 +187,11 @@ struct ei_solve_triangular_selector // FIXME this is cool but what about conjugate/adjoint expressions ? do we want to evaluate them ? // this is a more general problem though. ei_cache_friendly_product_colmajor_times_vector( - IsLower ? size-endBlock : endBlock+1, - &(lhs.const_cast_derived().coeffRef(IsLower ? endBlock : 0, IsLower ? startBlock : endBlock+1)), + IsLowerTriangular ? size-endBlock : endBlock+1, + &(lhs.const_cast_derived().coeffRef(IsLowerTriangular ? endBlock : 0, IsLowerTriangular ? startBlock : endBlock+1)), lhs.stride(), - btmp, &(other.coeffRef(IsLower ? endBlock : 0, c))); -// if (IsLower) + btmp, &(other.coeffRef(IsLowerTriangular ? endBlock : 0, c))); +// if (IsLowerTriangular) // other.col(c).end(size-endBlock) += (lhs.block(endBlock, startBlock, size-endBlock, endBlock-startBlock) // * other.col(c).block(startBlock,endBlock-startBlock)).lazy(); // else @@ -201,7 +201,7 @@ struct ei_solve_triangular_selector /* Now we have to process the remaining part as usual */ int i; - for(i=blockyEnd; IsLower ? i0; i += (IsLower ? 1 : -1) ) + for(i=blockyEnd; IsLowerTriangular ? i0; i += (IsLowerTriangular ? 1 : -1) ) { if(!(Lhs::Flags & UnitDiagBit)) other.coeffRef(i,c) /= lhs.coeff(i,i); @@ -209,7 +209,7 @@ struct ei_solve_triangular_selector /* NOTE we cannot use lhs.col(i).end(size-i-1) because Part::coeffRef gets called by .col() to * get the address of the start of the row */ - if(IsLower) + if(IsLowerTriangular) other.col(c).end(size-i-1) -= other.coeffRef(i,c) * Block(lhs, i+1,i, size-i-1,1); else other.col(c).start(i) -= other.coeffRef(i,c) * Block(lhs, 0,i, i, 1); diff --git a/Eigen/src/Core/Transpose.h b/Eigen/src/Core/Transpose.h index 2a9530124..1137e42dd 100644 --- a/Eigen/src/Core/Transpose.h +++ b/Eigen/src/Core/Transpose.h @@ -167,7 +167,7 @@ struct ei_inplace_transpose_selector; template struct ei_inplace_transpose_selector { // square matrix static void run(MatrixType& m) { - m.template part().swap(m.transpose()); + m.template part().swap(m.transpose()); } }; @@ -175,7 +175,7 @@ template struct ei_inplace_transpose_selector { // non square matrix static void run(MatrixType& m) { if (m.rows()==m.cols()) - m.template part().swap(m.transpose()); + m.template part().swap(m.transpose()); else m.set(m.transpose().eval()); } diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h index 20b427df1..203f3b294 100644 --- a/Eigen/src/Core/util/Constants.h +++ b/Eigen/src/Core/util/Constants.h @@ -185,16 +185,16 @@ const unsigned int HereditaryBits = RowMajorBit | SparseBit; // Possible values for the Mode parameter of part() and of extract() -const unsigned int Upper = UpperTriangularBit; -const unsigned int StrictlyUpper = UpperTriangularBit | ZeroDiagBit; -const unsigned int Lower = LowerTriangularBit; -const unsigned int StrictlyLower = LowerTriangularBit | ZeroDiagBit; +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; // additional possible values for the Mode parameter of extract() -const unsigned int UnitUpper = UpperTriangularBit | UnitDiagBit; -const unsigned int UnitLower = LowerTriangularBit | UnitDiagBit; -const unsigned int Diagonal = Upper | Lower; +const unsigned int UnitUpperTriangular = UpperTriangularBit | UnitDiagBit; +const unsigned int UnitLowerTriangular = LowerTriangularBit | UnitDiagBit; +const unsigned int Diagonal = UpperTriangular | LowerTriangular; enum { Aligned, Unaligned }; enum { ForceAligned, AsRequested }; diff --git a/Eigen/src/LU/LU.h b/Eigen/src/LU/LU.h index 41e33b48f..f3c93a4e8 100644 --- a/Eigen/src/LU/LU.h +++ b/Eigen/src/LU/LU.h @@ -114,7 +114,7 @@ template class LU * * \sa matrixLU(), matrixU() */ - inline const Part matrixL() const + inline const Part matrixL() const { return m_lu; } @@ -123,7 +123,7 @@ template class LU * * \sa matrixLU(), matrixL() */ - inline const Part matrixU() const + inline const Part matrixU() const { return m_lu; } @@ -441,7 +441,7 @@ void LU::computeKernel(KernelMatrixType *result) const y(-m_lu.corner(TopRight, m_rank, dimker)); m_lu.corner(TopLeft, m_rank, m_rank) - .template marked() + .template marked() .solveTriangularInPlace(y); for(int i = 0; i < m_rank; ++i) @@ -510,7 +510,7 @@ bool LU::solve( l.setZero(); l.corner(Eigen::TopLeft,rows,smalldim) = m_lu.corner(Eigen::TopLeft,rows,smalldim); - l.template marked().solveTriangularInPlace(c); + l.template marked().solveTriangularInPlace(c); // Step 3 if(!isSurjective()) @@ -527,7 +527,7 @@ bool LU::solve( MatrixType::MaxRowsAtCompileTime, OtherDerived::MaxColsAtCompileTime> d(c.corner(TopLeft, m_rank, c.cols())); m_lu.corner(TopLeft, m_rank, m_rank) - .template marked() + .template marked() .solveTriangularInPlace(d); // Step 4 diff --git a/Eigen/src/QR/HessenbergDecomposition.h b/Eigen/src/QR/HessenbergDecomposition.h index 21597bb02..6d0ff794e 100644 --- a/Eigen/src/QR/HessenbergDecomposition.h +++ b/Eigen/src/QR/HessenbergDecomposition.h @@ -243,7 +243,7 @@ HessenbergDecomposition::matrixH(void) const int n = m_matrix.rows(); MatrixType matH = m_matrix; if (n>2) - matH.corner(BottomLeft,n-2, n-2).template part().setZero(); + matH.corner(BottomLeft,n-2, n-2).template part().setZero(); return matH; } diff --git a/Eigen/src/QR/QR.h b/Eigen/src/QR/QR.h index c19205816..88914af3a 100644 --- a/Eigen/src/QR/QR.h +++ b/Eigen/src/QR/QR.h @@ -60,11 +60,11 @@ template class QR bool isFullRank() const { return ei_isMuchSmallerThan(m_hCoeffs.cwise().abs().minCoeff(), Scalar(1)); } /** \returns a read-only expression of the matrix R of the actual the QR decomposition */ - const Part, Upper> + const Part, UpperTriangular> matrixR(void) const { int cols = m_qr.cols(); - return MatrixRBlockType(m_qr, 0, 0, cols, cols).nestByValue().template part(); + return MatrixRBlockType(m_qr, 0, 0, cols, cols).nestByValue().template part(); } MatrixType matrixQ(void) const; diff --git a/Eigen/src/QR/SelfAdjointEigenSolver.h b/Eigen/src/QR/SelfAdjointEigenSolver.h index 95842c160..36188bcf9 100644 --- a/Eigen/src/QR/SelfAdjointEigenSolver.h +++ b/Eigen/src/QR/SelfAdjointEigenSolver.h @@ -254,22 +254,22 @@ compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors cholB.matrixL().solveTriangularInPlace(matC); // FIXME since we currently do not support A * inv(L'), let's do (inv(L) A')' : matC = matC.adjoint().eval(); - cholB.matrixL().template marked().solveTriangularInPlace(matC); + cholB.matrixL().template marked().solveTriangularInPlace(matC); matC = matC.adjoint().eval(); // this version works too: // matC = matC.transpose(); -// cholB.matrixL().conjugate().template marked().solveTriangularInPlace(matC); +// cholB.matrixL().conjugate().template marked().solveTriangularInPlace(matC); // matC = matC.transpose(); // FIXME: this should work: (currently it only does for small matrices) // Transpose trMatC(matC); -// cholB.matrixL().conjugate().eval().template marked().solveTriangularInPlace(trMatC); +// cholB.matrixL().conjugate().eval().template marked().solveTriangularInPlace(trMatC); compute(matC, computeEigenvectors); if (computeEigenvectors) { // transform back the eigen vectors: evecs = inv(U) * evecs - cholB.matrixL().adjoint().template marked().solveTriangularInPlace(m_eivec); + cholB.matrixL().adjoint().template marked().solveTriangularInPlace(m_eivec); for (int i=0; i::matrixT(void) const matT.corner(TopRight,n-1, n-1).diagonal() = subDiagonal().conjugate(); if (n>2) { - matT.corner(TopRight,n-2, n-2).template part().setZero(); - matT.corner(BottomLeft,n-2, n-2).template part().setZero(); + matT.corner(TopRight,n-2, n-2).template part().setZero(); + matT.corner(BottomLeft,n-2, n-2).template part().setZero(); } return matT; } @@ -223,17 +223,17 @@ void Tridiagonalization::_compute(MatrixType& matA, CoeffVectorType& /* This is the initial algorithm which minimize operation counts and maximize * the use of Eigen's expression. Unfortunately, the first matrix-vector product - * using Part is very very slow */ + * using Part is very very slow */ #ifdef EIGEN_NEVER_DEFINED // matrix - vector product - hCoeffs.end(n-i-1) = (matA.corner(BottomRight,n-i-1,n-i-1).template part() + hCoeffs.end(n-i-1) = (matA.corner(BottomRight,n-i-1,n-i-1).template part() * (h * matA.col(i).end(n-i-1))).lazy(); // simple axpy hCoeffs.end(n-i-1) += (h * Scalar(-0.5) * matA.col(i).end(n-i-1).dot(hCoeffs.end(n-i-1))) * matA.col(i).end(n-i-1); // rank-2 update //Block B(matA,i+1,i,n-i-1,1); - matA.corner(BottomRight,n-i-1,n-i-1).template part() -= + matA.corner(BottomRight,n-i-1,n-i-1).template part() -= (matA.col(i).end(n-i-1) * hCoeffs.end(n-i-1).adjoint()).lazy() + (hCoeffs.end(n-i-1) * matA.col(i).end(n-i-1).adjoint()).lazy(); #endif @@ -256,7 +256,7 @@ void Tridiagonalization::_compute(MatrixType& matA, CoeffVectorType& Block(matA,b+4,b,n-b-4,4).adjoint() * Block(matA,b+4,i,n-b-4,1); // the 4x4 block diagonal: Block(hCoeffs, b, 0, 4,1) += - (Block(matA,b,b,4,4).template part() + (Block(matA,b,b,4,4).template part() * (h * Block(matA,b,i,4,1))).lazy(); } #endif diff --git a/Eigen/src/Sparse/CholmodSupport.h b/Eigen/src/Sparse/CholmodSupport.h index c98ca2125..6689a4669 100644 --- a/Eigen/src/Sparse/CholmodSupport.h +++ b/Eigen/src/Sparse/CholmodSupport.h @@ -69,9 +69,9 @@ cholmod_sparse SparseMatrix::asCholmodMatrix() if (Flags & SelfAdjoint) { - if (Flags & Upper) + if (Flags & UpperTriangular) res.stype = 1; - else if (Flags & Lower) + else if (Flags & LowerTriangular) res.stype = -1; else res.stype = 0; diff --git a/Eigen/src/Sparse/SparseLDLT.h b/Eigen/src/Sparse/SparseLDLT.h index b29148864..b1f58b4b1 100644 --- a/Eigen/src/Sparse/SparseLDLT.h +++ b/Eigen/src/Sparse/SparseLDLT.h @@ -79,7 +79,7 @@ class SparseLDLT protected: typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::Real RealScalar; - typedef SparseMatrix CholMatrixType; + typedef SparseMatrix CholMatrixType; typedef Matrix VectorType; enum { diff --git a/Eigen/src/Sparse/SparseLLT.h b/Eigen/src/Sparse/SparseLLT.h index f9eb459dc..0b63f80ab 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::Real RealScalar; - typedef SparseMatrix CholMatrixType; + typedef SparseMatrix CholMatrixType; enum { SupernodalFactorIsDirty = 0x10000, diff --git a/Eigen/src/Sparse/SparseLU.h b/Eigen/src/Sparse/SparseLU.h index 4eb591fe7..142592050 100644 --- a/Eigen/src/Sparse/SparseLU.h +++ b/Eigen/src/Sparse/SparseLU.h @@ -41,7 +41,7 @@ class SparseLU protected: typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::Real RealScalar; - typedef SparseMatrix LUMatrixType; + typedef SparseMatrix LUMatrixType; enum { MatrixLUIsDirty = 0x10000 diff --git a/Eigen/src/Sparse/SuperLUSupport.h b/Eigen/src/Sparse/SuperLUSupport.h index aeab835aa..6f554fa6e 100644 --- a/Eigen/src/Sparse/SuperLUSupport.h +++ b/Eigen/src/Sparse/SuperLUSupport.h @@ -162,9 +162,9 @@ struct SluMatrixMapHelper > res.setScalarType(); // FIXME the following is not very accurate - if (Flags & Upper) + if (Flags & UpperTriangular) res.Mtype = SLU_TRU; - if (Flags & Lower) + if (Flags & LowerTriangular) res.Mtype = SLU_TRL; if (Flags & SelfAdjoint) ei_assert(false && "SelfAdjoint matrix shape not supported by SuperLU"); @@ -213,8 +213,8 @@ class SparseLU : public SparseLU typedef Matrix Vector; typedef Matrix IntRowVectorType; typedef Matrix IntColVectorType; - typedef SparseMatrix LMatrixType; - typedef SparseMatrix UMatrixType; + typedef SparseMatrix LMatrixType; + typedef SparseMatrix UMatrixType; using Base::m_flags; using Base::m_status; diff --git a/Eigen/src/Sparse/TaucsSupport.h b/Eigen/src/Sparse/TaucsSupport.h index 22fef1ceb..15fbb83ac 100644 --- a/Eigen/src/Sparse/TaucsSupport.h +++ b/Eigen/src/Sparse/TaucsSupport.h @@ -50,13 +50,13 @@ taucs_ccs_matrix SparseMatrix::asTaucsMatrix() ei_assert(false && "Scalar type not supported by TAUCS"); } - if (Flags & Upper) + if (Flags & UpperTriangular) res.flags |= TAUCS_UPPER; - if (Flags & Lower) + if (Flags & LowerTriangular) res.flags |= TAUCS_LOWER; if (Flags & SelfAdjoint) res.flags |= (NumTraits::IsComplex ? TAUCS_HERMITIAN : TAUCS_SYMMETRIC); - else if ((Flags & Upper) || (Flags & Lower)) + else if ((Flags & UpperTriangular) || (Flags & LowerTriangular)) res.flags |= TAUCS_TRIANGULAR; return res; diff --git a/Eigen/src/Sparse/TriangularSolver.h b/Eigen/src/Sparse/TriangularSolver.h index c32af328d..d44938273 100644 --- a/Eigen/src/Sparse/TriangularSolver.h +++ b/Eigen/src/Sparse/TriangularSolver.h @@ -27,7 +27,7 @@ // forward substitution, row-major template -struct ei_solve_triangular_selector +struct ei_solve_triangular_selector { typedef typename Rhs::Scalar Scalar; static void run(const Lhs& lhs, Rhs& other) @@ -59,7 +59,7 @@ struct ei_solve_triangular_selector // backward substitution, row-major template -struct ei_solve_triangular_selector +struct ei_solve_triangular_selector { typedef typename Rhs::Scalar Scalar; static void run(const Lhs& lhs, Rhs& other) @@ -92,7 +92,7 @@ struct ei_solve_triangular_selector // forward substitution, col-major template -struct ei_solve_triangular_selector +struct ei_solve_triangular_selector { typedef typename Rhs::Scalar Scalar; static void run(const Lhs& lhs, Rhs& other) @@ -120,7 +120,7 @@ struct ei_solve_triangular_selector // backward substitution, col-major template -struct ei_solve_triangular_selector +struct ei_solve_triangular_selector { typedef typename Rhs::Scalar Scalar; static void run(const Lhs& lhs, Rhs& other) diff --git a/Eigen/src/Sparse/UmfPackSupport.h b/Eigen/src/Sparse/UmfPackSupport.h index 00fa7e57d..1bb40a7c3 100644 --- a/Eigen/src/Sparse/UmfPackSupport.h +++ b/Eigen/src/Sparse/UmfPackSupport.h @@ -127,8 +127,8 @@ class SparseLU : public SparseLU typedef Matrix Vector; typedef Matrix IntRowVectorType; typedef Matrix IntColVectorType; - typedef SparseMatrix LMatrixType; - typedef SparseMatrix UMatrixType; + typedef SparseMatrix LMatrixType; + typedef SparseMatrix UMatrixType; using Base::m_flags; using Base::m_status; diff --git a/bench/btl/libs/eigen2/eigen2_interface.hh b/bench/btl/libs/eigen2/eigen2_interface.hh index 0f5b6f640..40d8f9c70 100644 --- a/bench/btl/libs/eigen2/eigen2_interface.hh +++ b/bench/btl/libs/eigen2/eigen2_interface.hh @@ -134,11 +134,11 @@ public : } static inline void trisolve_lower(const gene_matrix & L, const gene_vector& B, gene_vector& X, int N){ - X = L.template marked().solveTriangular(B); + X = L.template marked().solveTriangular(B); } static inline void trisolve_lower_matrix(const gene_matrix & L, const gene_matrix& B, gene_matrix& X, int N){ - X = L.template marked().solveTriangular(B); + X = L.template marked().solveTriangular(B); } static inline void cholesky(const gene_matrix & X, gene_matrix & C, int N){ diff --git a/bench/sparse_cholesky.cpp b/bench/sparse_cholesky.cpp index 5e07fac9b..ec8078e4c 100644 --- a/bench/sparse_cholesky.cpp +++ b/bench/sparse_cholesky.cpp @@ -37,8 +37,8 @@ X \ } timer.stop(); } -// typedef SparseMatrix EigenSparseTriMatrix; -typedef SparseMatrix EigenSparseSelfAdjointMatrix; +// typedef SparseMatrix EigenSparseTriMatrix; +typedef SparseMatrix EigenSparseSelfAdjointMatrix; void fillSpdMatrix(float density, int rows, int cols, EigenSparseSelfAdjointMatrix& dst) { diff --git a/bench/sparse_trisolver.cpp b/bench/sparse_trisolver.cpp index 021433043..158a381d8 100644 --- a/bench/sparse_trisolver.cpp +++ b/bench/sparse_trisolver.cpp @@ -34,8 +34,8 @@ X \ } timer.stop(); } -typedef SparseMatrix EigenSparseTriMatrix; -typedef SparseMatrix EigenSparseTriMatrixRow; +typedef SparseMatrix EigenSparseTriMatrix; +typedef SparseMatrix EigenSparseTriMatrixRow; void fillMatrix(float density, int rows, int cols, EigenSparseTriMatrix& dst) { @@ -83,11 +83,11 @@ int main(int argc, char *argv[]) eiToDense(sm1, m1); m2 = m1; - BENCH(x = m1.marked().solveTriangular(b);) + BENCH(x = m1.marked().solveTriangular(b);) std::cout << " colmajor^-1 * b:\t" << timer.value() << endl; // std::cerr << x.transpose() << "\n"; - BENCH(x = m2.marked().solveTriangular(b);) + BENCH(x = m2.marked().solveTriangular(b);) std::cout << " rowmajor^-1 * b:\t" << timer.value() << endl; // std::cerr << x.transpose() << "\n"; } diff --git a/doc/QuickStartGuide.dox b/doc/QuickStartGuide.dox index aa40606ed..81df4aa3d 100644 --- a/doc/QuickStartGuide.dox +++ b/doc/QuickStartGuide.dox @@ -513,20 +513,20 @@ Read/write access to special parts of a matrix can be achieved. See \link Matrix Extract triangular matrices \n from a given matrix m: \code -m.part() -m.part() -m.part() -m.part() -m.part() -m.part()\endcode +m.part() +m.part() +m.part() +m.part() +m.part() +m.part()\endcode Write to triangular parts \n of a matrix m: \code -m1.part() = m2; -m1.part() = m2; -m1.part() = m2; -m1.part() = m2;\endcode +m1.part() = m2; +m1.part() = m2; +m1.part() = m2; +m1.part() = m2;\endcode Special: take advantage of symmetry \n (selfadjointness) when copying \n an expression into a matrix diff --git a/doc/snippets/MatrixBase_extract.cpp b/doc/snippets/MatrixBase_extract.cpp index bedf4dfbd..02f44c83f 100644 --- a/doc/snippets/MatrixBase_extract.cpp +++ b/doc/snippets/MatrixBase_extract.cpp @@ -1,8 +1,8 @@ Matrix3i m = Matrix3i::Random(); cout << "Here is the matrix m:" << endl << m << endl; cout << "Here is the upper-triangular matrix extracted from m:" << endl - << m.part() << endl; + << m.part() << endl; cout << "Here is the strictly-upper-triangular matrix extracted from m:" << endl - << m.part() << endl; + << m.part() << endl; cout << "Here is the unit-lower-triangular matrix extracted from m:" << endl - << m.part() << endl; + << m.part() << endl; diff --git a/doc/snippets/MatrixBase_marked.cpp b/doc/snippets/MatrixBase_marked.cpp index d4c11c6fc..5c05ce81e 100644 --- a/doc/snippets/MatrixBase_marked.cpp +++ b/doc/snippets/MatrixBase_marked.cpp @@ -1,9 +1,9 @@ Matrix3d m = Matrix3d::Zero(); -m.part().setOnes(); +m.part().setOnes(); cout << "Here is the matrix m:" << endl << m << endl; Matrix3d n = Matrix3d::Ones(); -n.part() *= 2; +n.part() *= 2; cout << "Here is the matrix n:" << endl << n << endl; cout << "And now here is m.inverse()*n, taking advantage of the fact that" " m is upper-triangular:" << endl - << m.marked().solveTriangular(n); + << m.marked().solveTriangular(n); diff --git a/doc/snippets/MatrixBase_part.cpp b/doc/snippets/MatrixBase_part.cpp index fbc55b1ca..1abbd6888 100644 --- a/doc/snippets/MatrixBase_part.cpp +++ b/doc/snippets/MatrixBase_part.cpp @@ -1,5 +1,5 @@ Matrix3d m = Matrix3i::Zero(); -m.part().setOnes(); +m.part().setOnes(); cout << "Here is the matrix m:" << endl << m << endl; cout << "And let us now compute m*m.adjoint() in a very optimized way" << endl << "taking advantage of the symmetry." << endl; diff --git a/doc/snippets/class_LU_2.cpp b/doc/snippets/class_LU_2.cpp index ccf89a995..0c72d34a4 100644 --- a/doc/snippets/class_LU_2.cpp +++ b/doc/snippets/class_LU_2.cpp @@ -7,7 +7,7 @@ cout << "Here is, up to permutations, its LU decomposition matrix:" << endl << lu.matrixLU() << endl; cout << "Here is the actual L matrix in this decomposition:" << endl; Matrix5x5 l = Matrix5x5::Identity(); -l.block<5,3>(0,0).part() = lu.matrixLU(); +l.block<5,3>(0,0).part() = lu.matrixLU(); cout << l << endl; cout << "Let us now reconstruct the original matrix m:" << endl; Matrix5x3 x = l * lu.matrixU(); diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index fd39dc82e..4255128ad 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -131,9 +131,11 @@ macro(ei_add_test testname) target_link_libraries(${targetname} ${EXTERNAL_LIBS}) if(${ARGC} GREATER 2) - if(ARGV2) + string(STRIP "${ARGV2}" ARGV2_stripped) + string(LENGTH "${ARGV2_stripped}" ARGV2_stripped_length) + if(${ARGV2_stripped_length} GREATER 0) target_link_libraries(${targetname} ${ARGV2}) - endif(ARGV2) + endif(${ARGV2_stripped_length} GREATER 0) endif(${ARGC} GREATER 2) if(WIN32) diff --git a/test/sparse_solvers.cpp b/test/sparse_solvers.cpp index c3050ce7d..da9a245a2 100644 --- a/test/sparse_solvers.cpp +++ b/test/sparse_solvers.cpp @@ -44,13 +44,13 @@ template void sparse_solvers(int rows, int cols) // lower initSparse(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords); - VERIFY_IS_APPROX(refMat2.template marked().solveTriangular(vec2), - m2.template marked().solveTriangular(vec3)); + VERIFY_IS_APPROX(refMat2.template marked().solveTriangular(vec2), + m2.template marked().solveTriangular(vec3)); // upper initSparse(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, &zeroCoords, &nonzeroCoords); - VERIFY_IS_APPROX(refMat2.template marked().solveTriangular(vec2), - m2.template marked().solveTriangular(vec3)); + VERIFY_IS_APPROX(refMat2.template marked().solveTriangular(vec2), + m2.template marked().solveTriangular(vec3)); // TODO test row major } @@ -70,7 +70,7 @@ template void sparse_solvers(int rows, int cols) refMat2.diagonal() *= 0.5; refMat2.llt().solve(b, &refX); - typedef SparseMatrix SparseSelfAdjointMatrix; + typedef SparseMatrix SparseSelfAdjointMatrix; x = b; SparseLLT (m2).solveInPlace(x); //VERIFY(refX.isApprox(x,test_precision()) && "LLT: default"); @@ -107,7 +107,7 @@ template void sparse_solvers(int rows, int cols) refMat2.diagonal() *= 0.5; refMat2.ldlt().solve(b, &refX); - typedef SparseMatrix SparseSelfAdjointMatrix; + typedef SparseMatrix SparseSelfAdjointMatrix; x = b; SparseLDLT ldlt(m2); if (ldlt.succeeded()) diff --git a/test/triangular.cpp b/test/triangular.cpp index 34afa7b3c..7021e0776 100644 --- a/test/triangular.cpp +++ b/test/triangular.cpp @@ -51,14 +51,14 @@ template void triangular(const MatrixType& m) v2 = VectorType::Random(rows), vzero = VectorType::Zero(rows); - MatrixType m1up = m1.template part(); - MatrixType m2up = m2.template part(); + MatrixType m1up = m1.template part(); + MatrixType m2up = m2.template part(); if (rows*cols>1) { - VERIFY(m1up.isUpper()); - VERIFY(m2up.transpose().isLower()); - VERIFY(!m2.isLower()); + VERIFY(m1up.isUpperTriangular()); + VERIFY(m2up.transpose().isLowerTriangular()); + VERIFY(!m2.isLowerTriangular()); } // VERIFY_IS_APPROX(m1up.transpose() * m2, m1.upper().transpose().lower() * m2); @@ -66,20 +66,20 @@ template void triangular(const MatrixType& m) // test overloaded operator+= r1.setZero(); r2.setZero(); - r1.template part() += m1; + r1.template part() += m1; r2 += m1up; VERIFY_IS_APPROX(r1,r2); // test overloaded operator= m1.setZero(); - m1.template part() = (m2.transpose() * m2).lazy(); + m1.template part() = (m2.transpose() * m2).lazy(); m3 = m2.transpose() * m2; - VERIFY_IS_APPROX(m3.template part().transpose(), m1); + VERIFY_IS_APPROX(m3.template part().transpose(), m1); // test overloaded operator= m1.setZero(); - m1.template part() = (m2.transpose() * m2).lazy(); - VERIFY_IS_APPROX(m3.template part(), m1); + m1.template part() = (m2.transpose() * m2).lazy(); + VERIFY_IS_APPROX(m3.template part(), m1); VERIFY_IS_APPROX(m3.template part(), m3.diagonal().asDiagonal()); @@ -89,30 +89,30 @@ template void triangular(const MatrixType& m) Transpose trm4(m4); // test back and forward subsitution - m3 = m1.template part(); - VERIFY(m3.template marked().solveTriangular(m3).cwise().abs().isIdentity(test_precision())); - VERIFY(m3.transpose().template marked() + m3 = m1.template part(); + VERIFY(m3.template marked().solveTriangular(m3).cwise().abs().isIdentity(test_precision())); + VERIFY(m3.transpose().template marked() .solveTriangular(m3.transpose()).cwise().abs().isIdentity(test_precision())); // check M * inv(L) using in place API m4 = m3; - m3.transpose().template marked().solveTriangularInPlace(trm4); + m3.transpose().template marked().solveTriangularInPlace(trm4); VERIFY(m4.cwise().abs().isIdentity(test_precision())); - m3 = m1.template part(); - VERIFY(m3.template marked().solveTriangular(m3).cwise().abs().isIdentity(test_precision())); - VERIFY(m3.transpose().template marked() + m3 = m1.template part(); + VERIFY(m3.template marked().solveTriangular(m3).cwise().abs().isIdentity(test_precision())); + VERIFY(m3.transpose().template marked() .solveTriangular(m3.transpose()).cwise().abs().isIdentity(test_precision())); // check M * inv(U) using in place API m4 = m3; - m3.transpose().template marked().solveTriangularInPlace(trm4); + m3.transpose().template marked().solveTriangularInPlace(trm4); VERIFY(m4.cwise().abs().isIdentity(test_precision())); - m3 = m1.template part(); - VERIFY(m2.isApprox(m3 * (m3.template marked().solveTriangular(m2)), largerEps)); - m3 = m1.template part(); - VERIFY(m2.isApprox(m3 * (m3.template marked().solveTriangular(m2)), largerEps)); + m3 = m1.template part(); + VERIFY(m2.isApprox(m3 * (m3.template marked().solveTriangular(m2)), largerEps)); + m3 = m1.template part(); + VERIFY(m2.isApprox(m3 * (m3.template marked().solveTriangular(m2)), largerEps)); - VERIFY((m1.template part() * m2.template part()).isUpper()); + VERIFY((m1.template part() * m2.template part()).isUpperTriangular()); } -- cgit v1.2.3