GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-06-02 17:59:46 +08:00
Code Issues Packages Projects Releases Wiki Activity

remove the Triangular suffix to Upper, Lower, UnitLower, etc,

Browse Source 
and remove the respective bit flags
This commit is contained in:
Gael Guennebaud 2010-01-07 21:15:32 +01:00
parent 82ec250a0f
commit c5d7c9f0de
58 changed files with 591 additions and 563 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
Eigen/Core
View File

@ -23,10 +23,6 @@
// License and a copy of the GNU General Public License along with // License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>. // Eigen. If not, see <http://www.gnu.org/licenses/>.
#ifdef EIGEN2_SUPPORT
#include "Eigen2Support"
#endif
#ifndef EIGEN_CORE_H #ifndef EIGEN_CORE_H
#define EIGEN_CORE_H #define EIGEN_CORE_H
@ -225,5 +221,9 @@ struct Dense {};
#include "src/Core/util/EnableMSVCWarnings.h" #include "src/Core/util/EnableMSVCWarnings.h"
#ifdef EIGEN2_SUPPORT
#include "Eigen2Support"
#endif
#endif // EIGEN_CORE_H #endif // EIGEN_CORE_H
/* vim: set filetype=cpp et sw=2 ts=2 ai: */ /* vim: set filetype=cpp et sw=2 ts=2 ai: */

12
Eigen/Eigen2Support
View File

@ -25,14 +25,13 @@
#ifndef EIGEN2SUPPORT_H #ifndef EIGEN2SUPPORT_H
#define EIGEN2SUPPORT_H #define EIGEN2SUPPORT_H
#include "src/Core/util/DisableMSVCWarnings.h" #include "Array"
#ifndef EIGEN2_SUPPORT #if (!defined(EIGEN2_SUPPORT)) || (!defined(EIGEN_CORE_H))
#define EIGEN2_SUPPORT #error Eigen2 support must be enabled by defining EIGEN2_SUPPORT before any other Eigen header
#endif #endif
#include "Core" #include "src/Core/util/DisableMSVCWarnings.h"
#include "Array"
namespace Eigen { namespace Eigen {
@ -40,8 +39,9 @@ namespace Eigen {
* This module provides a couple of deprecated functions improving the compatibility with Eigen2. * This module provides a couple of deprecated functions improving the compatibility with Eigen2.
* *
* \code * \code
* #include <Eigen/Eigen2Support> * #define EIGEN2_SUPPORT
* \endcode * \endcode
*
*/ */
#include "src/Eigen2Support/Flagged.h" #include "src/Eigen2Support/Flagged.h"

4
Eigen/src/Array/Reverse.h
View File

@ -59,9 +59,7 @@ struct ei_traits<Reverse<MatrixType, Direction> >
LinearAccess = ( (Direction==BothDirections) && (int(_MatrixTypeNested::Flags)&PacketAccessBit) ) LinearAccess = ( (Direction==BothDirections) && (int(_MatrixTypeNested::Flags)&PacketAccessBit) )
? LinearAccessBit : 0, ? LinearAccessBit : 0,
Flags = (int(_MatrixTypeNested::Flags) & (HereditaryBits | PacketAccessBit | LinearAccess)) Flags = (int(_MatrixTypeNested::Flags) & (HereditaryBits | PacketAccessBit | LinearAccess)),
| (int(_MatrixTypeNested::Flags)&UpperTriangularBit ? LowerTriangularBit : 0)
| (int(_MatrixTypeNested::Flags)&LowerTriangularBit ? UpperTriangularBit : 0),
CoeffReadCost = _MatrixTypeNested::CoeffReadCost CoeffReadCost = _MatrixTypeNested::CoeffReadCost
}; };

6
Eigen/src/Cholesky/LDLT.h
View File

@ -80,7 +80,7 @@ template<typename _MatrixType> class LDLT
} }
/** \returns the lower triangular matrix L */ /** \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."); ei_assert(m_isInitialized && "LDLT is not initialized.");
return m_matrix; return m_matrix;
@ -301,13 +301,13 @@ bool LDLT<MatrixType>::solveInPlace(MatrixBase<Derived> &bAndX) const
// y = L^-1 z // y = L^-1 z
//matrixL().solveInPlace(bAndX); //matrixL().solveInPlace(bAndX);
m_matrix.template triangularView<UnitLowerTriangular>().solveInPlace(bAndX); m_matrix.template triangularView<UnitLower>().solveInPlace(bAndX);
// w = D^-1 y // w = D^-1 y
bAndX = m_matrix.diagonal().asDiagonal().inverse() * bAndX; bAndX = m_matrix.diagonal().asDiagonal().inverse() * bAndX;
// u = L^-T w // u = L^-T w
m_matrix.adjoint().template triangularView<UnitUpperTriangular>().solveInPlace(bAndX); m_matrix.adjoint().template triangularView<UnitUpper>().solveInPlace(bAndX);
// x = P^T u // x = P^T u
for (int i = size-1; i >= 0; --i) bAndX.row(m_transpositions.coeff(i)).swap(bAndX.row(i)); for (int i = size-1; i >= 0; --i) bAndX.row(m_transpositions.coeff(i)).swap(bAndX.row(i));

28
Eigen/src/Cholesky/LLT.h
View File

@ -148,7 +148,7 @@ template<typename _MatrixType, int _UpLo> class LLT
// forward declaration (defined at the end of this file) // forward declaration (defined at the end of this file)
template<int UpLo> struct ei_llt_inplace; template<int UpLo> struct ei_llt_inplace;
template<> struct ei_llt_inplace<LowerTriangular> template<> struct ei_llt_inplace<Lower>
{ {
template<typename MatrixType> template<typename MatrixType>
static bool unblocked(MatrixType& mat) 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); Block<MatrixType,Dynamic,Dynamic> A22(m,k+bs,k+bs,rs,rs);
if(!unblocked(A11)) return false; if(!unblocked(A11)) return false;
if(rs>0) A11.adjoint().template triangularView<UpperTriangular>().template solveInPlace<OnTheRight>(A21); if(rs>0) A11.adjoint().template triangularView<Upper>().template solveInPlace<OnTheRight>(A21);
if(rs>0) A22.template selfadjointView<LowerTriangular>().rankUpdate(A21,-1); // bottleneck if(rs>0) A22.template selfadjointView<Lower>().rankUpdate(A21,-1); // bottleneck
} }
return true; return true;
} }
}; };
template<> struct ei_llt_inplace<UpperTriangular> template<> struct ei_llt_inplace<Upper>
{ {
template<typename MatrixType> template<typename MatrixType>
static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat) static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat)
{ {
Transpose<MatrixType> matt(mat); Transpose<MatrixType> matt(mat);
return ei_llt_inplace<LowerTriangular>::unblocked(matt); return ei_llt_inplace<Lower>::unblocked(matt);
} }
template<typename MatrixType> template<typename MatrixType>
static EIGEN_STRONG_INLINE bool blocked(MatrixType& mat) static EIGEN_STRONG_INLINE bool blocked(MatrixType& mat)
{ {
Transpose<MatrixType> matt(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<MatrixType, Lower> MatrixL;
typedef TriangularView<typename MatrixType::AdjointReturnType, UpperTriangular> MatrixU; typedef TriangularView<typename MatrixType::AdjointReturnType, Upper> MatrixU;
inline static MatrixL getL(const MatrixType& m) { return m; } inline static MatrixL getL(const MatrixType& m) { return m; }
inline static MatrixU getU(const MatrixType& m) { return m.adjoint(); } inline static MatrixU getU(const MatrixType& m) { return m.adjoint(); }
static bool inplace_decomposition(MatrixType& m) 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<typename MatrixType::AdjointReturnType, Lower> MatrixL;
typedef TriangularView<MatrixType, UpperTriangular> MatrixU; typedef TriangularView<MatrixType, Upper> MatrixU;
inline static MatrixL getL(const MatrixType& m) { return m.adjoint(); } inline static MatrixL getL(const MatrixType& m) { return m.adjoint(); }
inline static MatrixU getU(const MatrixType& m) { return m; } inline static MatrixU getU(const MatrixType& m) { return m; }
static bool inplace_decomposition(MatrixType& 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 /** Computes / recomputes the Cholesky decomposition A = LL^* = U^*U of \a matrix

7
Eigen/src/Core/MatrixBase.h
View File

@ -405,13 +405,6 @@ template<typename Derived> class MatrixBase
template<int Size> const VectorBlock<Derived,Size> start() const; template<int Size> const VectorBlock<Derived,Size> start() const;
template<int Size> VectorBlock<Derived,Size> end(); template<int Size> VectorBlock<Derived,Size> end();
template<int Size> const VectorBlock<Derived,Size> end() const; 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 #endif
protected: protected:

4
Eigen/src/Core/SelfAdjointView.h
View File

@ -31,7 +31,7 @@
* \brief Expression of a selfadjoint matrix from a triangular part of a dense matrix * \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 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 * 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() * 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 typename ei_unref<MatrixTypeNested>::type _MatrixTypeNested;
typedef MatrixType ExpressionType; typedef MatrixType ExpressionType;
enum { enum {
Mode = TriangularPart | SelfAdjointBit, Mode = TriangularPart | SelfAdjoint,
Flags = _MatrixTypeNested::Flags & (HereditaryBits) Flags = _MatrixTypeNested::Flags & (HereditaryBits)
& (~(PacketAccessBit | DirectAccessBit | LinearAccessBit)), // FIXME these flags should be preserved & (~(PacketAccessBit | DirectAccessBit | LinearAccessBit)), // FIXME these flags should be preserved
CoeffReadCost = _MatrixTypeNested::CoeffReadCost CoeffReadCost = _MatrixTypeNested::CoeffReadCost

56
Eigen/src/Core/SolveTriangular.h
View File

@ -57,7 +57,7 @@ struct ei_triangular_solver_selector<Lhs,Rhs,OnTheLeft,Mode,NoUnrolling,RowMajor
typedef ei_blas_traits<Lhs> LhsProductTraits; typedef ei_blas_traits<Lhs> LhsProductTraits;
typedef typename LhsProductTraits::ExtractType ActualLhsType; typedef typename LhsProductTraits::ExtractType ActualLhsType;
enum { enum {
IsLowerTriangular = ((Mode&LowerTriangularBit)==LowerTriangularBit) IsLower = ((Mode&Lower)==Lower)
}; };
static void run(const Lhs& lhs, Rhs& other) 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); ActualLhsType actualLhs = LhsProductTraits::extract(lhs);
const int size = lhs.cols(); const int size = lhs.cols();
for(int pi=IsLowerTriangular ? 0 : size; for(int pi=IsLower ? 0 : size;
IsLowerTriangular ? pi<size : pi>0; IsLower ? pi<size : pi>0;
IsLowerTriangular ? pi+=PanelWidth : pi-=PanelWidth) 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) if (r > 0)
{ {
// let's directly call the low level product function because: // let's directly call the low level product function because:
// 1 - it is faster to compile // 1 - it is faster to compile
// 2 - it is slighlty faster at runtime // 2 - it is slighlty faster at runtime
int startRow = IsLowerTriangular ? pi : pi-actualPanelWidth; int startRow = IsLower ? pi : pi-actualPanelWidth;
int startCol = IsLowerTriangular ? 0 : pi; int startCol = IsLower ? 0 : pi;
VectorBlock<Rhs,Dynamic> target(other,startRow,actualPanelWidth); VectorBlock<Rhs,Dynamic> target(other,startRow,actualPanelWidth);
ei_cache_friendly_product_rowmajor_times_vector<LhsProductTraits::NeedToConjugate,false>( 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) for(int k=0; k<actualPanelWidth; ++k)
{ {
int i = IsLowerTriangular ? pi+k : pi-k-1; int i = IsLower ? pi+k : pi-k-1;
int s = IsLowerTriangular ? pi : i+1; int s = IsLower ? pi : i+1;
if (k>0) if (k>0)
other.coeffRef(i) -= (lhs.row(i).segment(s,k).transpose().cwiseProduct(other.segment(s,k))).sum(); 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); 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; typedef typename LhsProductTraits::ExtractType ActualLhsType;
enum { enum {
PacketSize = ei_packet_traits<Scalar>::size, PacketSize = ei_packet_traits<Scalar>::size,
IsLowerTriangular = ((Mode&LowerTriangularBit)==LowerTriangularBit) IsLower = ((Mode&Lower)==Lower)
}; };
static void run(const Lhs& lhs, Rhs& other) 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); ActualLhsType actualLhs = LhsProductTraits::extract(lhs);
const int size = lhs.cols(); const int size = lhs.cols();
for(int pi=IsLowerTriangular ? 0 : size; for(int pi=IsLower ? 0 : size;
IsLowerTriangular ? pi<size : pi>0; IsLower ? pi<size : pi>0;
IsLowerTriangular ? pi+=PanelWidth : pi-=PanelWidth) IsLower ? pi+=PanelWidth : pi-=PanelWidth)
{ {
int actualPanelWidth = std::min(IsLowerTriangular ? size - pi : pi, PanelWidth); int actualPanelWidth = std::min(IsLower ? size - pi : pi, PanelWidth);
int startBlock = IsLowerTriangular ? pi : pi-actualPanelWidth; int startBlock = IsLower ? pi : pi-actualPanelWidth;
int endBlock = IsLowerTriangular ? pi + actualPanelWidth : 0; int endBlock = IsLower ? pi + actualPanelWidth : 0;
for(int k=0; k<actualPanelWidth; ++k) for(int k=0; k<actualPanelWidth; ++k)
{ {
int i = IsLowerTriangular ? pi+k : pi-k-1; int i = IsLower ? pi+k : pi-k-1;
if(!(Mode & UnitDiagBit)) if(!(Mode & UnitDiag))
other.coeffRef(i) /= lhs.coeff(i,i); other.coeffRef(i) /= lhs.coeff(i,i);
int r = actualPanelWidth - k - 1; // remaining size int r = actualPanelWidth - k - 1; // remaining size
int s = IsLowerTriangular ? i+1 : i-r; int s = IsLower ? i+1 : i-r;
if (r>0) if (r>0)
other.segment(s,r) -= other.coeffRef(i) * Block<Lhs,Dynamic,1>(lhs, s, i, r, 1); 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) if (r > 0)
{ {
// let's directly call the low level product function because: // 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> template<typename Lhs, typename Rhs, int Mode, int Index, int Size>
struct ei_triangular_solver_unroller<Lhs,Rhs,Mode,Index,Size,false> { struct ei_triangular_solver_unroller<Lhs,Rhs,Mode,Index,Size,false> {
enum { enum {
IsLowerTriangular = ((Mode&LowerTriangularBit)==LowerTriangularBit), IsLower = ((Mode&Lower)==Lower),
I = IsLowerTriangular ? Index : Size - Index - 1, I = IsLower ? Index : Size - Index - 1,
S = IsLowerTriangular ? 0 : I+1 S = IsLower ? 0 : I+1
}; };
static void run(const Lhs& lhs, Rhs& rhs) static void run(const Lhs& lhs, Rhs& rhs)
{ {
if (Index>0) if (Index>0)
rhs.coeffRef(I) -= ((lhs.row(I).template segment<Index>(S).transpose()).cwiseProduct(rhs.template segment<Index>(S))).sum(); 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); rhs.coeffRef(I) /= lhs.coeff(I,I);
ei_triangular_solver_unroller<Lhs,Rhs,Mode,Index+1,Size>::run(lhs,rhs); 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(); OtherDerived& other = _other.const_cast_derived();
ei_assert(cols() == rows()); ei_assert(cols() == rows());
ei_assert( (Side==OnTheLeft && cols() == other.rows()) || (Side==OnTheRight && cols() == other.cols()) ); ei_assert( (Side==OnTheLeft && cols() == other.rows()) || (Side==OnTheRight && cols() == other.cols()) );
ei_assert(!(Mode & ZeroDiagBit)); ei_assert(!(Mode & ZeroDiag));
ei_assert(Mode & (UpperTriangularBit|LowerTriangularBit)); ei_assert(Mode & (Upper|Lower));
enum { copy = ei_traits<OtherDerived>::Flags & RowMajorBit && OtherDerived::IsVectorAtCompileTime }; enum { copy = ei_traits<OtherDerived>::Flags & RowMajorBit && OtherDerived::IsVectorAtCompileTime };
typedef typename ei_meta_if<copy, typedef typename ei_meta_if<copy,

9
Eigen/src/Core/Transpose.h
View File

@ -49,10 +49,7 @@ struct ei_traits<Transpose<MatrixType> > : ei_traits<MatrixType>
ColsAtCompileTime = MatrixType::RowsAtCompileTime, ColsAtCompileTime = MatrixType::RowsAtCompileTime,
MaxRowsAtCompileTime = MatrixType::MaxColsAtCompileTime, MaxRowsAtCompileTime = MatrixType::MaxColsAtCompileTime,
MaxColsAtCompileTime = MatrixType::MaxRowsAtCompileTime, MaxColsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
Flags = ((int(_MatrixTypeNested::Flags) ^ RowMajorBit) Flags = (int(_MatrixTypeNested::Flags) ^ RowMajorBit),
& ~(LowerTriangularBit | UpperTriangularBit))
| (int(_MatrixTypeNested::Flags)&UpperTriangularBit ? LowerTriangularBit : 0)
| (int(_MatrixTypeNested::Flags)&LowerTriangularBit ? UpperTriangularBit : 0),
CoeffReadCost = _MatrixTypeNested::CoeffReadCost CoeffReadCost = _MatrixTypeNested::CoeffReadCost
}; };
}; };
@ -229,7 +226,7 @@ struct ei_inplace_transpose_selector;
template<typename MatrixType> template<typename MatrixType>
struct ei_inplace_transpose_selector<MatrixType,true> { // square matrix struct ei_inplace_transpose_selector<MatrixType,true> { // square matrix
static void run(MatrixType& m) { 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 struct ei_inplace_transpose_selector<MatrixType,false> { // non square matrix
static void run(MatrixType& m) { static void run(MatrixType& m) {
if (m.rows()==m.cols()) if (m.rows()==m.cols())
m.template triangularView<StrictlyUpperTriangular>().swap(m.transpose()); m.template triangularView<StrictlyUpper>().swap(m.transpose());
else else
m = m.transpose().eval(); m = m.transpose().eval();
} }

91
Eigen/src/Core/TriangularMatrix.h
View File

@ -46,7 +46,7 @@ template<typename Derived> class TriangularBase : public AnyMatrixBase<Derived>
}; };
typedef typename ei_traits<Derived>::Scalar Scalar; 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 rows() const { return derived().rows(); }
inline int cols() const { return derived().cols(); } 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) void check_coordinates(int row, int col)
{ {
ei_assert(col>=0 && col<cols() && row>=0 && row<rows()); ei_assert(col>=0 && col<cols() && row>=0 && row<rows());
ei_assert( (Mode==UpperTriangular && col>=row) ei_assert( (Mode==Upper && col>=row)
|| (Mode==LowerTriangular && col<=row) || (Mode==Lower && col<=row)
|| ((Mode==StrictlyUpperTriangular || Mode==UnitUpperTriangular) && col>row) || ((Mode==StrictlyUpper || Mode==UnitUpper) && col>row)
|| ((Mode==StrictlyLowerTriangular || Mode==UnitLowerTriangular) && col<row)); || ((Mode==StrictlyLower || Mode==UnitLower) && col<row));
} }
#ifdef EIGEN_INTERNAL_DEBUGGING #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 * \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 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, * \param Mode the kind of triangular matrix expression to construct. Can be Upper,
* LowerTriangular, UpperSelfadjoint, or LowerSelfadjoint. This is in fact a bit field; * Lower, UpperSelfadjoint, or LowerSelfadjoint. This is in fact a bit field;
* it must have either UpperBit or LowerBit, and additionnaly it may have either * it must have either Upper or Lower, and additionnaly it may have either
* TraingularBit or SelfadjointBit. * UnitDiag or Selfadjoint.
* *
* This class represents a triangular part of a matrix, not necessarily square. Strictly speaking, for rectangular * 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 * 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 { enum {
Mode = _Mode, Mode = _Mode,
TransposeMode = (Mode & UpperTriangularBit ? LowerTriangularBit : 0) TransposeMode = (Mode & Upper ? Lower : 0)
| (Mode & LowerTriangularBit ? UpperTriangularBit : 0) | (Mode & Lower ? Upper : 0)
| (Mode & (ZeroDiagBit | UnitDiagBit)) | (Mode & (UnitDiag))
| (Mode & (ZeroDiag))
}; };
inline TriangularView(const MatrixType& matrix) : m_matrix(matrix) 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 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); return SelfAdjointView<_MatrixTypeNested,Mode>(m_matrix);
} }
SelfAdjointView<_MatrixTypeNested,Mode> selfadjointView() 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); 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()); 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: protected:
const MatrixTypeNested m_matrix; 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_triangular_assignment_selector<Derived1, Derived2, Mode, UnrollCount-1, ClearOpposite>::run(dst, src);
ei_assert( Mode == UpperTriangular || Mode == LowerTriangular ei_assert( Mode == Upper || Mode == Lower
|| Mode == StrictlyUpperTriangular || Mode == StrictlyLowerTriangular || Mode == StrictlyUpper || Mode == StrictlyLower
|| Mode == UnitUpperTriangular || Mode == UnitLowerTriangular); || Mode == UnitUpper || Mode == UnitLower);
if((Mode == UpperTriangular && row <= col) if((Mode == Upper && row <= col)
|| (Mode == LowerTriangular && row >= col) || (Mode == Lower && row >= col)
|| (Mode == StrictlyUpperTriangular && row < col) || (Mode == StrictlyUpper && row < col)
|| (Mode == StrictlyLowerTriangular && row > col) || (Mode == StrictlyLower && row > col)
|| (Mode == UnitUpperTriangular && row < col) || (Mode == UnitUpper && row < col)
|| (Mode == UnitLowerTriangular && row > col)) || (Mode == UnitLower && row > col))
dst.copyCoeff(row, col, src); dst.copyCoeff(row, col, src);
else if(ClearOpposite) else if(ClearOpposite)
{ {
if (Mode&UnitDiagBit && row==col) if (Mode&UnitDiag && row==col)
dst.coeffRef(row, col) = 1; dst.coeffRef(row, col) = 1;
else else
dst.coeffRef(row, col) = 0; 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> 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) 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> 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) 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> 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) 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> 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) 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> 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) 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> 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) inline static void run(Derived1 &dst, const Derived2 &src)
{ {
@ -590,8 +601,8 @@ EIGEN_DEPRECATED TriangularView<Derived, Mode> MatrixBase<Derived>::part()
/** \nonstableyet /** \nonstableyet
* \returns an expression of a triangular view extracted from the current matrix * \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, * The parameter \a Mode can have the following values: \c Upper, \c StrictlyUpper, \c UnitUpper,
* \c LowerTriangular, \c StrictlyLowerTriangular, \c UnitLowerTriangular. * \c Lower, \c StrictlyLower, \c UnitLower.
* *
* Example: \include MatrixBase_extract.cpp * Example: \include MatrixBase_extract.cpp
* Output: \verbinclude MatrixBase_extract.out * 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, /** \returns true if *this is approximately equal to an upper triangular matrix,
* within the precision given by \a prec. * within the precision given by \a prec.
* *
* \sa isLowerTriangular(), extract(), part(), marked() * \sa isLower(), extract(), part(), marked()
*/ */
template<typename Derived> template<typename Derived>
bool MatrixBase<Derived>::isUpperTriangular(RealScalar prec) const 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) for(int j = 0; j < cols(); ++j)
{ {
int maxi = std::min(j, rows()-1); int maxi = std::min(j, rows()-1);
for(int i = 0; i <= maxi; ++i) for(int i = 0; i <= maxi; ++i)
{ {
RealScalar absValue = ei_abs(coeff(i,j)); 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 j = 0; j < cols(); ++j)
for(int i = j+1; i < rows(); ++i) for(int i = j+1; i < rows(); ++i)
if(ei_abs(coeff(i, j)) > threshold) return false; 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, /** \returns true if *this is approximately equal to a lower triangular matrix,
* within the precision given by \a prec. * within the precision given by \a prec.
* *
* \sa isUpperTriangular(), extract(), part(), marked() * \sa isUpper(), extract(), part(), marked()
*/ */
template<typename Derived> template<typename Derived>
bool MatrixBase<Derived>::isLowerTriangular(RealScalar prec) const 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 j = 0; j < cols(); ++j)
for(int i = j; i < rows(); ++i) for(int i = j; i < rows(); ++i)
{ {
RealScalar absValue = ei_abs(coeff(i,j)); 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) for(int j = 1; j < cols(); ++j)
{ {
int maxi = std::min(j, rows()-1); int maxi = std::min(j, rows()-1);

16
Eigen/src/Core/products/SelfadjointMatrixMatrix.h
View File

@ -368,10 +368,10 @@ struct SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,RhsMode,false>
SelfadjointProductMatrix(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) {} SelfadjointProductMatrix(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) {}
enum { enum {
LhsUpLo = LhsMode&(UpperTriangularBit|LowerTriangularBit), LhsUpLo = LhsMode&(Upper|Lower),
LhsIsSelfAdjoint = (LhsMode&SelfAdjointBit)==SelfAdjointBit, LhsIsSelfAdjoint = (LhsMode&SelfAdjoint)==SelfAdjoint,
RhsUpLo = RhsMode&(UpperTriangularBit|LowerTriangularBit), RhsUpLo = RhsMode&(Upper|Lower),
RhsIsSelfAdjoint = (RhsMode&SelfAdjointBit)==SelfAdjointBit RhsIsSelfAdjoint = (RhsMode&SelfAdjoint)==SelfAdjoint
}; };
template<typename Dest> void scaleAndAddTo(Dest& dst, Scalar alpha) const 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); * RhsBlasTraits::extractScalarFactor(m_rhs);
ei_product_selfadjoint_matrix<Scalar, ei_product_selfadjoint_matrix<Scalar,
EIGEN_LOGICAL_XOR(LhsUpLo==UpperTriangular, EIGEN_LOGICAL_XOR(LhsUpLo==Upper,
ei_traits<Lhs>::Flags &RowMajorBit) ? RowMajor : ColMajor, LhsIsSelfAdjoint, ei_traits<Lhs>::Flags &RowMajorBit) ? RowMajor : ColMajor, LhsIsSelfAdjoint,
NumTraits<Scalar>::IsComplex && EIGEN_LOGICAL_XOR(LhsUpLo==UpperTriangular,bool(LhsBlasTraits::NeedToConjugate)), NumTraits<Scalar>::IsComplex && EIGEN_LOGICAL_XOR(LhsUpLo==Upper,bool(LhsBlasTraits::NeedToConjugate)),
EIGEN_LOGICAL_XOR(RhsUpLo==UpperTriangular, EIGEN_LOGICAL_XOR(RhsUpLo==Upper,
ei_traits<Rhs>::Flags &RowMajorBit) ? RowMajor : ColMajor, RhsIsSelfAdjoint, 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> ei_traits<Dest>::Flags&RowMajorBit ? RowMajor : ColMajor>
::run( ::run(
lhs.rows(), rhs.cols(), // sizes lhs.rows(), rhs.cols(), // sizes

4
Eigen/src/Core/products/SelfadjointMatrixVector.h
View File

@ -42,7 +42,7 @@ static EIGEN_DONT_INLINE void ei_product_selfadjoint_vector(
enum { enum {
IsRowMajor = StorageOrder==RowMajor ? 1 : 0, IsRowMajor = StorageOrder==RowMajor ? 1 : 0,
IsLower = UpLo == LowerTriangularBit ? 1 : 0, IsLower = UpLo == Lower ? 1 : 0,
FirstTriangular = IsRowMajor == IsLower FirstTriangular = IsRowMajor == IsLower
}; };
@ -170,7 +170,7 @@ struct SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,0,true>
EIGEN_PRODUCT_PUBLIC_INTERFACE(SelfadjointProductMatrix) EIGEN_PRODUCT_PUBLIC_INTERFACE(SelfadjointProductMatrix)
enum { enum {
LhsUpLo = LhsMode&(UpperTriangularBit|LowerTriangularBit) LhsUpLo = LhsMode&(Upper|Lower)
}; };
SelfadjointProductMatrix(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) {} SelfadjointProductMatrix(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) {}

14
Eigen/src/Core/products/SelfadjointProduct.h
View File

@ -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) 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); ::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 // 1 - before the diagonal => processed with gebp or skipped
// 2 - the actual_mc x actual_mc symmetric block => processed with a special kernel // 2 - the actual_mc x actual_mc symmetric block => processed with a special kernel
// 3 - after the diagonal => processed with gebp or skipped // 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)); 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>() 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); (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; 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)); 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); int actualBlockSize = std::min<int>(BlockSize,size - j);
const Scalar* actual_b = blockB+j*depth*PacketSize; 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); gebp_kernel(res+j*resStride, resStride, blockA, actual_b, j, depth, actualBlockSize);
// selfadjoint micro block // selfadjoint micro block
@ -186,13 +186,13 @@ struct ei_sybb_kernel
for(int j1=0; j1<actualBlockSize; ++j1) for(int j1=0; j1<actualBlockSize; ++j1)
{ {
Scalar* r = res + (j+j1)*resStride + i; Scalar* r = res + (j+j1)*resStride + i;
for(int i1=UpLo==LowerTriangular ? j1 : 0; for(int i1=UpLo==Lower ? j1 : 0;
UpLo==LowerTriangular ? i1<actualBlockSize : i1<=j1; ++i1) UpLo==Lower ? i1<actualBlockSize : i1<=j1; ++i1)
r[i1] += buffer(i1,j1); r[i1] += buffer(i1,j1);
} }
} }
if(UpLo==LowerTriangular) if(UpLo==Lower)
{ {
int i = j+actualBlockSize; int i = j+actualBlockSize;
gebp_kernel(res+j*resStride+i, resStride, blockA+depth*i, actual_b, size-i, depth, actualBlockSize); gebp_kernel(res+j*resStride+i, resStride, blockA+depth*i, actual_b, size-i, depth, actualBlockSize);

6
Eigen/src/Core/products/SelfadjointRank2Update.h
View File

@ -33,7 +33,7 @@ template<typename Scalar, typename UType, typename VType, int UpLo>
struct ei_selfadjoint_rank2_update_selector; struct ei_selfadjoint_rank2_update_selector;
template<typename Scalar, typename UType, typename VType> 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) 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> 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) 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, 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 ^ UBlasTraits::NeedToConjugate,_ActualUType>::ret>::type,
typename ei_cleantype<typename ei_conj_expr_if<IsRowMajor ^ VBlasTraits::NeedToConjugate,_ActualVType>::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); ::run(const_cast<Scalar*>(_expression().data()),_expression().stride(),actualU,actualV,actualAlpha);
return *this; return *this;

56
Eigen/src/Core/products/TriangularMatrixMatrix.h
View File

@ -75,7 +75,7 @@ struct ei_product_triangular_matrix_matrix<Scalar,Mode,LhsIsTriangular,
Scalar alpha) Scalar alpha)
{ {
ei_product_triangular_matrix_matrix<Scalar, ei_product_triangular_matrix_matrix<Scalar,
(Mode&UnitDiagBit) | (Mode&UpperTriangular) ? LowerTriangular : UpperTriangular, (Mode&UnitDiag) | (Mode&Upper) ? Lower : Upper,
(!LhsIsTriangular), (!LhsIsTriangular),
RhsStorageOrder==RowMajor ? ColMajor : RowMajor, RhsStorageOrder==RowMajor ? ColMajor : RowMajor,
ConjugateRhs, ConjugateRhs,
@ -113,7 +113,7 @@ struct ei_product_triangular_matrix_matrix<Scalar,Mode,true,
typedef ei_product_blocking_traits<Scalar> Blocking; typedef ei_product_blocking_traits<Scalar> Blocking;
enum { enum {
SmallPanelWidth = EIGEN_ENUM_MAX(Blocking::mr,Blocking::nr), 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 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_lhs<Scalar,Blocking::mr,LhsStorageOrder> pack_lhs;
ei_gemm_pack_rhs<Scalar,Blocking::nr,RhsStorageOrder> pack_rhs; ei_gemm_pack_rhs<Scalar,Blocking::nr,RhsStorageOrder> pack_rhs;
for(int k2=IsLowerTriangular ? size : 0; for(int k2=IsLower ? size : 0;
IsLowerTriangular ? k2>0 : k2<size; IsLower ? k2>0 : k2<size;
IsLowerTriangular ? k2-=kc : k2+=kc) IsLower ? k2-=kc : k2+=kc)
{ {
const int actual_kc = std::min(IsLowerTriangular ? k2 : size-k2, kc); const int actual_kc = std::min(IsLower ? k2 : size-k2, kc);
int actual_k2 = IsLowerTriangular ? k2-actual_kc : k2; int actual_k2 = IsLower ? k2-actual_kc : k2;
pack_rhs(blockB, &rhs(actual_k2,0), rhsStride, alpha, actual_kc, cols); 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) for (int k1=0; k1<actual_kc; k1+=SmallPanelWidth)
{ {
int actualPanelWidth = std::min<int>(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 startBlock = actual_k2+k1;
int blockBOffset = 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 // To this end we do an extra triangular copy to small temporary buffer
for (int k=0;k<actualPanelWidth;++k) for (int k=0;k<actualPanelWidth;++k)
{ {
if (!(Mode&UnitDiagBit)) if (!(Mode&UnitDiag))
triangularBuffer.coeffRef(k,k) = lhs(startBlock+k,startBlock+k); 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); triangularBuffer.coeffRef(i,k) = lhs(startBlock+i,startBlock+k);
} }
pack_lhs(blockA, triangularBuffer.data(), triangularBuffer.stride(), actualPanelWidth, actualPanelWidth); 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 // GEBP with remaining micro panel
if (lengthTarget>0) 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); 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 // the part below the diagonal => GEPP
{ {
int start = IsLowerTriangular ? k2 : 0; int start = IsLower ? k2 : 0;
int end = IsLowerTriangular ? size : actual_k2; int end = IsLower ? size : actual_k2;
for(int i2=start; i2<end; i2+=mc) for(int i2=start; i2<end; i2+=mc)
{ {
const int actual_mc = std::min(i2+mc,end)-i2; 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; typedef ei_product_blocking_traits<Scalar> Blocking;
enum { enum {
SmallPanelWidth = EIGEN_ENUM_MAX(Blocking::mr,Blocking::nr), 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 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> pack_rhs;
ei_gemm_pack_rhs<Scalar,Blocking::nr,RhsStorageOrder,true> pack_rhs_panel; ei_gemm_pack_rhs<Scalar,Blocking::nr,RhsStorageOrder,true> pack_rhs_panel;
for(int k2=IsLowerTriangular ? 0 : size; for(int k2=IsLower ? 0 : size;
IsLowerTriangular ? k2<size : k2>0; IsLower ? k2<size : k2>0;
IsLowerTriangular ? k2+=kc : k2-=kc) 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);
int actual_k2 = IsLowerTriangular ? k2 : k2-actual_kc; int actual_k2 = IsLower ? k2 : k2-actual_kc;
int rs = IsLowerTriangular ? actual_k2 : size - k2; int rs = IsLower ? actual_k2 : size - k2;
Scalar* geb = blockB+actual_kc*actual_kc*Blocking::PacketSize; 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 // 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 actualPanelWidth = std::min<int>(actual_kc-j2, SmallPanelWidth);
int actual_j2 = actual_k2 + j2; int actual_j2 = actual_k2 + j2;
int panelOffset = IsLowerTriangular ? j2+actualPanelWidth : 0; int panelOffset = IsLower ? j2+actualPanelWidth : 0;
int panelLength = IsLowerTriangular ? actual_kc-j2-actualPanelWidth : j2; int panelLength = IsLower ? actual_kc-j2-actualPanelWidth : j2;
// general part // general part
pack_rhs_panel(blockB+j2*actual_kc*Blocking::PacketSize, pack_rhs_panel(blockB+j2*actual_kc*Blocking::PacketSize,
&rhs(actual_k2+panelOffset, actual_j2), rhsStride, alpha, &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 // append the triangular part via a temporary buffer
for (int j=0;j<actualPanelWidth;++j) for (int j=0;j<actualPanelWidth;++j)
{ {
if (!(Mode&UnitDiagBit)) if (!(Mode&UnitDiag))
triangularBuffer.coeffRef(j,j) = rhs(actual_j2+j,actual_j2+j); 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); 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) for (int j2=0; j2<actual_kc; j2+=SmallPanelWidth)
{ {
int actualPanelWidth = std::min<int>(actual_kc-j2, SmallPanelWidth); int actualPanelWidth = std::min<int>(actual_kc-j2, SmallPanelWidth);
int panelLength = IsLowerTriangular ? actual_kc-j2 : j2+actualPanelWidth; int panelLength = IsLower ? actual_kc-j2 : j2+actualPanelWidth;
int blockOffset = IsLowerTriangular ? j2 : 0; int blockOffset = IsLower ? j2 : 0;
gebp_kernel(res+i2+(actual_k2+j2)*resStride, resStride, gebp_kernel(res+i2+(actual_k2+j2)*resStride, resStride,
blockA, blockB+j2*actual_kc*Blocking::PacketSize, 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 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); blockA, geb, actual_mc, actual_kc, rs);
} }
} }

24
Eigen/src/Core/products/TriangularMatrixVector.h
View File

@ -34,8 +34,8 @@ struct ei_product_triangular_vector_selector<Lhs,Rhs,Result,Mode,ConjLhs,ConjRhs
{ {
typedef typename Rhs::Scalar Scalar; typedef typename Rhs::Scalar Scalar;
enum { enum {
IsLowerTriangular = ((Mode&LowerTriangularBit)==LowerTriangularBit), IsLower = ((Mode&Lower)==Lower),
HasUnitDiag = (Mode & UnitDiagBit)==UnitDiagBit HasUnitDiag = (Mode & UnitDiag)==UnitDiag
}; };
static EIGEN_DONT_INLINE void run(const Lhs& lhs, const Rhs& rhs, Result& res, typename ei_traits<Lhs>::Scalar alpha) 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) for (int k=0; k<actualPanelWidth; ++k)
{ {
int i = pi + k; int i = pi + k;
int s = IsLowerTriangular ? (HasUnitDiag ? i+1 : i ) : pi; int s = IsLower ? (HasUnitDiag ? i+1 : i ) : pi;
int r = IsLowerTriangular ? actualPanelWidth-k : k+1; int r = IsLower ? actualPanelWidth-k : k+1;
if ((!HasUnitDiag) || (--r)>0) if ((!HasUnitDiag) || (--r)>0)
res.segment(s,r) += (alpha * cjRhs.coeff(i)) * cjLhs.col(i).segment(s,r); res.segment(s,r) += (alpha * cjRhs.coeff(i)) * cjLhs.col(i).segment(s,r);
if (HasUnitDiag) if (HasUnitDiag)
res.coeffRef(i) += alpha * cjRhs.coeff(i); res.coeffRef(i) += alpha * cjRhs.coeff(i);
} }
int r = IsLowerTriangular ? size - pi - actualPanelWidth : pi; int r = IsLower ? size - pi - actualPanelWidth : pi;
if (r>0) if (r>0)
{ {
int s = IsLowerTriangular ? pi+actualPanelWidth : 0; int s = IsLower ? pi+actualPanelWidth : 0;
ei_cache_friendly_product_colmajor_times_vector<ConjLhs,ConjRhs>( ei_cache_friendly_product_colmajor_times_vector<ConjLhs,ConjRhs>(
r, r,
&(lhs.const_cast_derived().coeffRef(s,pi)), lhs.stride(), &(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; typedef typename Rhs::Scalar Scalar;
enum { enum {
IsLowerTriangular = ((Mode&LowerTriangularBit)==LowerTriangularBit), IsLower = ((Mode&Lower)==Lower),
HasUnitDiag = (Mode & UnitDiagBit)==UnitDiagBit HasUnitDiag = (Mode & UnitDiag)==UnitDiag
}; };
static void run(const Lhs& lhs, const Rhs& rhs, Result& res, typename ei_traits<Lhs>::Scalar alpha) 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) for (int k=0; k<actualPanelWidth; ++k)
{ {
int i = pi + k; int i = pi + k;
int s = IsLowerTriangular ? pi : (HasUnitDiag ? i+1 : i); int s = IsLower ? pi : (HasUnitDiag ? i+1 : i);
int r = IsLowerTriangular ? k+1 : actualPanelWidth-k; int r = IsLower ? k+1 : actualPanelWidth-k;
if ((!HasUnitDiag) || (--r)>0) if ((!HasUnitDiag) || (--r)>0)
res.coeffRef(i) += alpha * (cjLhs.row(i).segment(s,r).cwiseProduct(cjRhs.segment(s,r).transpose())).sum(); res.coeffRef(i) += alpha * (cjLhs.row(i).segment(s,r).cwiseProduct(cjRhs.segment(s,r).transpose())).sum();
if (HasUnitDiag) if (HasUnitDiag)
res.coeffRef(i) += alpha * cjRhs.coeff(i); res.coeffRef(i) += alpha * cjRhs.coeff(i);
} }
int r = IsLowerTriangular ? pi : size - pi - actualPanelWidth; int r = IsLower ? pi : size - pi - actualPanelWidth;
if (r>0) 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); Block<Result,Dynamic,1> target(res,pi,0,actualPanelWidth,1);
ei_cache_friendly_product_rowmajor_times_vector<ConjLhs,ConjRhs>( ei_cache_friendly_product_rowmajor_times_vector<ConjLhs,ConjRhs>(
&(lhs.const_cast_derived().coeffRef(pi,s)), lhs.stride(), &(lhs.const_cast_derived().coeffRef(pi,s)), lhs.stride(),

70
Eigen/src/Core/products/TriangularSolverMatrix.h
View File

@ -36,7 +36,7 @@ struct ei_triangular_solve_matrix<Scalar,Side,Mode,Conjugate,TriStorageOrder,Row
{ {
ei_triangular_solve_matrix< ei_triangular_solve_matrix<
Scalar, Side==OnTheLeft?OnTheRight:OnTheLeft, Scalar, Side==OnTheLeft?OnTheRight:OnTheLeft,
(Mode&UnitDiagBit) | ((Mode&UpperTriangular) ? LowerTriangular : UpperTriangular), (Mode&UnitDiag) | ((Mode&Upper) ? Lower : Upper),
NumTraits<Scalar>::IsComplex && Conjugate, NumTraits<Scalar>::IsComplex && Conjugate,
TriStorageOrder==RowMajor ? ColMajor : RowMajor, ColMajor> TriStorageOrder==RowMajor ? ColMajor : RowMajor, ColMajor>
::run(size, cols, tri, triStride, _other, otherStride); ::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; typedef ei_product_blocking_traits<Scalar> Blocking;
enum { enum {
SmallPanelWidth = EIGEN_ENUM_MAX(Blocking::mr,Blocking::nr), 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 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_gebp_kernel<Scalar, Blocking::mr, Blocking::nr, ei_conj_helper<Conjugate,false> > gebp_kernel;
ei_gemm_pack_lhs<Scalar,Blocking::mr,TriStorageOrder> pack_lhs; ei_gemm_pack_lhs<Scalar,Blocking::mr,TriStorageOrder> pack_lhs;
for(int k2=IsLowerTriangular ? 0 : size; for(int k2=IsLower ? 0 : size;
IsLowerTriangular ? k2<size : k2>0; IsLower ? k2<size : k2>0;
IsLowerTriangular ? k2+=kc : k2-=kc) 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, // 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 // 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) for (int k=0; k<actualPanelWidth; ++k)
{ {
// TODO write a small kernel handling this (can be shared with trsv) // TODO write a small kernel handling this (can be shared with trsv)
int i = IsLowerTriangular ? k2+k1+k : k2-k1-k-1; int i = IsLower ? k2+k1+k : k2-k1-k-1;
int s = IsLowerTriangular ? k2+k1 : i+1; int s = IsLower ? k2+k1 : i+1;
int rs = actualPanelWidth - k - 1; // remaining size 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) for (int j=0; j<cols; ++j)
{ {
if (TriStorageOrder==RowMajor) if (TriStorageOrder==RowMajor)
@ -120,7 +120,7 @@ struct ei_triangular_solve_matrix<Scalar,OnTheLeft,Mode,Conjugate,TriStorageOrde
} }
else else
{ {
int s = IsLowerTriangular ? i+1 : i-rs; int s = IsLower ? i+1 : i-rs;
Scalar b = (other(i,j) *= a); Scalar b = (other(i,j) *= a);
Scalar* r = &other(s,j); Scalar* r = &other(s,j);
const Scalar* l = &tri(s,i); 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 lengthTarget = actual_kc-k1-actualPanelWidth;
int startBlock = IsLowerTriangular ? k2+k1 : k2-k1-actualPanelWidth; int startBlock = IsLower ? k2+k1 : k2-k1-actualPanelWidth;
int blockBOffset = IsLowerTriangular ? k1 : lengthTarget; int blockBOffset = IsLower ? k1 : lengthTarget;
// update the respective rows of B from other // update the respective rows of B from other
ei_gemm_pack_rhs<Scalar, Blocking::nr, ColMajor, true>() ei_gemm_pack_rhs<Scalar, Blocking::nr, ColMajor, true>()
@ -141,7 +141,7 @@ struct ei_triangular_solve_matrix<Scalar,OnTheLeft,Mode,Conjugate,TriStorageOrde
// GEBP // GEBP
if (lengthTarget>0) 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); 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 // R2 = A2 * B => GEPP
{ {
int start = IsLowerTriangular ? k2+kc : 0; int start = IsLower ? k2+kc : 0;
int end = IsLowerTriangular ? size : k2-kc; int end = IsLower ? size : k2-kc;
for(int i2=start; i2<end; i2+=mc) for(int i2=start; i2<end; i2+=mc)
{ {
const int actual_mc = std::min(mc,end-i2); const int actual_mc = std::min(mc,end-i2);
if (actual_mc>0) 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); 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 { enum {
RhsStorageOrder = TriStorageOrder, RhsStorageOrder = TriStorageOrder,
SmallPanelWidth = EIGEN_ENUM_MAX(Blocking::mr,Blocking::nr), 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 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_rhs<Scalar,Blocking::nr,RhsStorageOrder,true> pack_rhs_panel;
ei_gemm_pack_lhs<Scalar, Blocking::mr, ColMajor, false, true> pack_lhs_panel; ei_gemm_pack_lhs<Scalar, Blocking::mr, ColMajor, false, true> pack_lhs_panel;
for(int k2=IsLowerTriangular ? size : 0; for(int k2=IsLower ? size : 0;
IsLowerTriangular ? k2>0 : k2<size; IsLower ? k2>0 : k2<size;
IsLowerTriangular ? k2-=kc : k2+=kc) IsLower ? k2-=kc : k2+=kc)
{ {
const int actual_kc = std::min(IsLowerTriangular ? k2 : size-k2, kc); const int actual_kc = std::min(IsLower ? k2 : size-k2, kc);
int actual_k2 = IsLowerTriangular ? k2-actual_kc : k2 ; int actual_k2 = IsLower ? k2-actual_kc : k2 ;
int startPanel = IsLowerTriangular ? 0 : k2+actual_kc; int startPanel = IsLower ? 0 : k2+actual_kc;
int rs = IsLowerTriangular ? actual_k2 : size - actual_k2 - actual_kc; int rs = IsLower ? actual_k2 : size - actual_k2 - actual_kc;
Scalar* geb = blockB+actual_kc*actual_kc*Blocking::PacketSize; Scalar* geb = blockB+actual_kc*actual_kc*Blocking::PacketSize;
if (rs>0) pack_rhs(geb, &rhs(actual_k2,startPanel), triStride, -1, actual_kc, rs); 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 actualPanelWidth = std::min<int>(actual_kc-j2, SmallPanelWidth);
int actual_j2 = actual_k2 + j2; int actual_j2 = actual_k2 + j2;
int panelOffset = IsLowerTriangular ? j2+actualPanelWidth : 0; int panelOffset = IsLower ? j2+actualPanelWidth : 0;
int panelLength = IsLowerTriangular ? actual_kc-j2-actualPanelWidth : j2; int panelLength = IsLower ? actual_kc-j2-actualPanelWidth : j2;
if (panelLength>0) if (panelLength>0)
pack_rhs_panel(blockB+j2*actual_kc*Blocking::PacketSize, 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 // triangular solver kernel
{ {
// for each small block of the diagonal (=> vertical panels of rhs) // 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) ? (actual_kc - ((actual_kc%SmallPanelWidth) ? (actual_kc%SmallPanelWidth)
: SmallPanelWidth)) : SmallPanelWidth))
: 0; : 0;
IsLowerTriangular ? j2>=0 : j2<actual_kc; IsLower ? j2>=0 : j2<actual_kc;
IsLowerTriangular ? j2-=SmallPanelWidth : j2+=SmallPanelWidth) IsLower ? j2-=SmallPanelWidth : j2+=SmallPanelWidth)
{ {
int actualPanelWidth = std::min<int>(actual_kc-j2, SmallPanelWidth); int actualPanelWidth = std::min<int>(actual_kc-j2, SmallPanelWidth);
int absolute_j2 = actual_k2 + j2; int absolute_j2 = actual_k2 + j2;
int panelOffset = IsLowerTriangular ? j2+actualPanelWidth : 0; int panelOffset = IsLower ? j2+actualPanelWidth : 0;
int panelLength = IsLowerTriangular ? actual_kc - j2 - actualPanelWidth : j2; int panelLength = IsLower ? actual_kc - j2 - actualPanelWidth : j2;
// GEBP // GEBP
if(panelLength>0) if(panelLength>0)
@ -269,17 +269,17 @@ struct ei_triangular_solve_matrix<Scalar,OnTheRight,Mode,Conjugate,TriStorageOrd
// unblocked triangular solve // unblocked triangular solve
for (int k=0; k<actualPanelWidth; ++k) 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); Scalar* r = &lhs(i2,j);
for (int k3=0; k3<k; ++k3) for (int k3=0; k3<k; ++k3)
{ {
Scalar b = conj(rhs(IsLowerTriangular ? j+1+k3 : absolute_j2+k3,j)); Scalar b = conj(rhs(IsLower ? j+1+k3 : absolute_j2+k3,j));
Scalar* a = &lhs(i2,IsLowerTriangular ? j+1+k3 : absolute_j2+k3); Scalar* a = &lhs(i2,IsLower ? j+1+k3 : absolute_j2+k3);
for (int i=0; i<actual_mc; ++i) for (int i=0; i<actual_mc; ++i)
r[i] -= a[i] * b; 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) for (int i=0; i<actual_mc; ++i)
r[i] *= b; r[i] *= b;
} }

38
Eigen/src/Core/util/Constants.h
View File

@ -149,44 +149,18 @@ const unsigned int DirectAccessBit = 0x20;
* means the first coefficient packet is guaranteed to be aligned */ * means the first coefficient packet is guaranteed to be aligned */
const unsigned int AlignedBit = 0x40; 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 // list of flags that are inherited by default
const unsigned int HereditaryBits = RowMajorBit const unsigned int HereditaryBits = RowMajorBit
| EvalBeforeNestingBit | EvalBeforeNestingBit
| EvalBeforeAssigningBit; | EvalBeforeAssigningBit;
// Possible values for the Mode parameter of part() // Possible values for the Mode parameter of triangularView()
const unsigned int UpperTriangular = UpperTriangularBit; enum {
const unsigned int StrictlyUpperTriangular = UpperTriangularBit | ZeroDiagBit; Lower=0x1, Upper=0x2, UnitDiag=0x4, ZeroDiag=0x8,
const unsigned int LowerTriangular = LowerTriangularBit; UnitLower=UnitDiag|Lower, UnitUpper=UnitDiag|Upper,
const unsigned int StrictlyLowerTriangular = LowerTriangularBit | ZeroDiagBit; StrictlyLower=ZeroDiag|Lower, StrictlyUpper=ZeroDiag|Upper,
const unsigned int SelfAdjoint = SelfAdjointBit; SelfAdjoint=0x10};
const unsigned int UnitUpperTriangular = UpperTriangularBit | UnitDiagBit;
const unsigned int UnitLowerTriangular = LowerTriangularBit | UnitDiagBit;
enum { Unaligned=0, Aligned=1 }; enum { Unaligned=0, Aligned=1 };
enum { ConditionalJumpCost = 5 }; enum { ConditionalJumpCost = 5 };

2
Eigen/src/Core/util/ForwardDeclarations.h
View File

@ -132,7 +132,7 @@ template<typename MatrixType> class ColPivHouseholderQR;
template<typename MatrixType> class FullPivHouseholderQR; template<typename MatrixType> class FullPivHouseholderQR;
template<typename MatrixType> class SVD; template<typename MatrixType> class SVD;
template<typename MatrixType, unsigned int Options = 0> class JacobiSVD; 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 MatrixType> class LDLT;
template<typename VectorsType, typename CoeffsType> class HouseholderSequence; template<typename VectorsType, typename CoeffsType> class HouseholderSequence;
template<typename Scalar> class PlanarRotation; template<typename Scalar> class PlanarRotation;

3
Eigen/src/Core/util/XprHelper.h
View File

@ -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 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 /** \internal Helper base class to add a scalar multiple operator

6
Eigen/src/Eigen2Support/Flagged.h
View File

@ -109,6 +109,12 @@ template<typename ExpressionType, unsigned int Added, unsigned int Removed> clas
const ExpressionType& _expression() const { return m_matrix; } 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: protected:
ExpressionTypeNested m_matrix; ExpressionTypeNested m_matrix;
}; };

53
Eigen/src/Eigen2Support/TriangularSolver.h
View File

@ -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

2
Eigen/src/Eigenvalues/HessenbergDecomposition.h
View File

@ -199,7 +199,7 @@ HessenbergDecomposition<MatrixType>::matrixH() const
int n = m_matrix.rows(); int n = m_matrix.rows();
MatrixType matH = m_matrix; MatrixType matH = m_matrix;
if (n>2) 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; return matH;
} }

8
Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
View File

@ -247,11 +247,11 @@ compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors
matC.adjointInPlace(); matC.adjointInPlace();
// this version works too: // this version works too:
// matC = matC.transpose(); // matC = matC.transpose();
// cholB.matrixL().conjugate().template marked<LowerTriangular>().solveTriangularInPlace(matC); // cholB.matrixL().conjugate().template marked<Lower>().solveTriangularInPlace(matC);
// matC = matC.transpose(); // matC = matC.transpose();
// FIXME: this should work: (currently it only does for small matrices) // FIXME: this should work: (currently it only does for small matrices)
// Transpose<MatrixType> trMatC(matC); // 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); 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> inline Matrix<typename NumTraits<typename ei_traits<Derived>::Scalar>::Real, ei_traits<Derived>::ColsAtCompileTime, 1>
MatrixBase<Derived>::eigenvalues() const MatrixBase<Derived>::eigenvalues() const
{ {
ei_assert(Flags&SelfAdjointBit); ei_assert(Flags&SelfAdjoint);
return SelfAdjointEigenSolver<typename Derived::PlainMatrixType>(eval(),false).eigenvalues(); return SelfAdjointEigenSolver<typename Derived::PlainMatrixType>(eval(),false).eigenvalues();
} }
@ -316,7 +316,7 @@ template<typename Derived>
inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real
MatrixBase<Derived>::operatorNorm() const MatrixBase<Derived>::operatorNorm() const
{ {
return ei_operatorNorm_selector<Derived, Flags&SelfAdjointBit> return ei_operatorNorm_selector<Derived, Flags&SelfAdjoint>
::operatorNorm(derived()); ::operatorNorm(derived());
} }

8
Eigen/src/Eigenvalues/Tridiagonalization.h
View File

@ -172,8 +172,8 @@ Tridiagonalization<MatrixType>::matrixT(void) const
matT.corner(TopRight,n-1, n-1).diagonal() = subDiagonal().template cast<Scalar>().conjugate(); matT.corner(TopRight,n-1, n-1).diagonal() = subDiagonal().template cast<Scalar>().conjugate();
if (n>2) if (n>2)
{ {
matT.corner(TopRight,n-2, n-2).template triangularView<UpperTriangular>().setZero(); matT.corner(TopRight,n-2, n-2).template triangularView<Upper>().setZero();
matT.corner(BottomLeft,n-2, n-2).template triangularView<LowerTriangular>().setZero(); matT.corner(BottomLeft,n-2, n-2).template triangularView<Lower>().setZero();
} }
return matT; 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) // 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; 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))); * (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); 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); .rankUpdate(matA.col(i).tail(remainingSize), hCoeffs.tail(remainingSize), -1);
matA.col(i).coeffRef(i+1) = beta; matA.col(i).coeffRef(i+1) = beta;

21
Eigen/src/LU/Determinant.h
View File

@ -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)); * (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, template<typename Derived,
int DeterminantType = int DeterminantType = Derived::RowsAtCompileTime
(Derived::Flags & (UpperTriangularBit | LowerTriangularBit))
? TriangularDeterminant : Derived::RowsAtCompileTime
> struct ei_determinant_impl > struct ei_determinant_impl
{ {
static inline typename ei_traits<Derived>::Scalar run(const Derived& m) 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> template<typename Derived> struct ei_determinant_impl<Derived, 1>
{ {
static inline typename ei_traits<Derived>::Scalar run(const Derived& m) static inline typename ei_traits<Derived>::Scalar run(const Derived& m)

8
Eigen/src/LU/FullPivLU.h
View File

@ -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.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());
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) for(int i = 0; i < rank(); ++i)
m.col(i).swap(m.col(pivots.coeff(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 // 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. // negative of the RHS, but for performance we just put the negative sign elsewhere, see below.
m.corner(TopLeft, rank(), rank()) m.corner(TopLeft, rank(), rank())
.template triangularView<UpperTriangular>().solveInPlace( .template triangularView<Upper>().solveInPlace(
m.corner(TopRight, rank(), dimker) m.corner(TopRight, rank(), dimker)
); );
@ -638,7 +638,7 @@ struct ei_solve_retval<FullPivLU<_MatrixType>, Rhs>
// Step 2 // Step 2
dec().matrixLU() dec().matrixLU()
.corner(Eigen::TopLeft,smalldim,smalldim) .corner(Eigen::TopLeft,smalldim,smalldim)
.template triangularView<UnitLowerTriangular>() .template triangularView<UnitLower>()
.solveInPlace(c.corner(Eigen::TopLeft, smalldim, c.cols())); .solveInPlace(c.corner(Eigen::TopLeft, smalldim, c.cols()));
if(rows>cols) if(rows>cols)
{ {
@ -650,7 +650,7 @@ struct ei_solve_retval<FullPivLU<_MatrixType>, Rhs>
// Step 3 // Step 3
dec().matrixLU() dec().matrixLU()
.corner(TopLeft, nonzero_pivots, nonzero_pivots) .corner(TopLeft, nonzero_pivots, nonzero_pivots)
.template triangularView<UpperTriangular>() .template triangularView<Upper>()
.solveInPlace(c.corner(TopLeft, nonzero_pivots, c.cols())); .solveInPlace(c.corner(TopLeft, nonzero_pivots, c.cols()));
// Step 4 // Step 4

6
Eigen/src/LU/PartialPivLU.h
View File

@ -349,7 +349,7 @@ struct ei_partial_lu_impl
A_2.row(i).swap(A_2.row(row_transpositions[i])); A_2.row(i).swap(A_2.row(row_transpositions[i]));
// A12 = A11^-1 A12 // A12 = A11^-1 A12
A11.template triangularView<UnitLowerTriangular>().solveInPlace(A12); A11.template triangularView<UnitLower>().solveInPlace(A12);
A22.noalias() -= A21 * A12; A22.noalias() -= A21 * A12;
} }
@ -426,10 +426,10 @@ struct ei_solve_retval<PartialPivLU<_MatrixType>, Rhs>
dst = dec().permutationP() * rhs(); dst = dec().permutationP() * rhs();
// Step 2 // Step 2
dec().matrixLU().template triangularView<UnitLowerTriangular>().solveInPlace(dst); dec().matrixLU().template triangularView<UnitLower>().solveInPlace(dst);
// Step 3 // Step 3
dec().matrixLU().template triangularView<UpperTriangular>().solveInPlace(dst); dec().matrixLU().template triangularView<Upper>().solveInPlace(dst);
} }
}; };

6
Eigen/src/QR/ColPivHouseholderQR.h
View File

@ -381,7 +381,7 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const
m_nonzero_pivots = k; m_nonzero_pivots = k;
m_hCoeffs.tail(size-k).setZero(); m_hCoeffs.tail(size-k).setZero();
m_qr.corner(BottomRight,rows-k,cols-k) m_qr.corner(BottomRight,rows-k,cols-k)
.template triangularView<StrictlyLowerTriangular>() .template triangularView<StrictlyLower>()
.setZero(); .setZero();
break; break;
} }
@ -453,7 +453,7 @@ struct ei_solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs>
dec().matrixQR() dec().matrixQR()
.corner(TopLeft, nonzero_pivots, nonzero_pivots) .corner(TopLeft, nonzero_pivots, nonzero_pivots)
.template triangularView<UpperTriangular>() .template triangularView<Upper>()
.solveInPlace(c.corner(TopLeft, nonzero_pivots, c.cols())); .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()) d.corner(TopLeft, nonzero_pivots, c.cols())
= dec().matrixQR() = dec().matrixQR()
.corner(TopLeft, nonzero_pivots, nonzero_pivots) .corner(TopLeft, nonzero_pivots, nonzero_pivots)
.template triangularView<UpperTriangular>() .template triangularView<Upper>()
* c.corner(TopLeft, nonzero_pivots, c.cols()); * 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); for(int i = 0; i < nonzero_pivots; ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i);

2
Eigen/src/QR/FullPivHouseholderQR.h
View File

@ -376,7 +376,7 @@ struct ei_solve_retval<FullPivHouseholderQR<_MatrixType>, Rhs>
} }
dec().matrixQR() dec().matrixQR()
.corner(TopLeft, dec().rank(), dec().rank()) .corner(TopLeft, dec().rank(), dec().rank())
.template triangularView<UpperTriangular>() .template triangularView<Upper>()
.solveInPlace(c.corner(TopLeft, dec().rank(), c.cols())); .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); for(int i = 0; i < dec().rank(); ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i);

2
Eigen/src/QR/HouseholderQR.h
View File

@ -231,7 +231,7 @@ struct ei_solve_retval<HouseholderQR<_MatrixType>, Rhs>
dec().matrixQR() dec().matrixQR()
.corner(TopLeft, rank, rank) .corner(TopLeft, rank, rank)
.template triangularView<UpperTriangular>() .template triangularView<Upper>()
.solveInPlace(c.corner(TopLeft, rank, c.cols())); .solveInPlace(c.corner(TopLeft, rank, c.cols()));
dst.corner(TopLeft, rank, c.cols()) = c.corner(TopLeft, rank, c.cols()); dst.corner(TopLeft, rank, c.cols()) = c.corner(TopLeft, rank, c.cols());

4
Eigen/src/SVD/JacobiSVD.h
View File

@ -234,7 +234,7 @@ struct ei_svd_precondition_if_more_rows_than_cols<MatrixType, Options, true>
if(rows > cols) if(rows > cols)
{ {
FullPivHouseholderQR<MatrixType> qr(matrix); 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(ComputeU) svd.m_matrixU = qr.matrixQ();
if(ComputeV) svd.m_matrixV = qr.colsPermutation(); 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> MatrixOptions,MaxColsAtCompileTime,MaxRowsAtCompileTime>
TransposeTypeWithSameStorageOrder; TransposeTypeWithSameStorageOrder;
FullPivHouseholderQR<TransposeTypeWithSameStorageOrder> qr(matrix.adjoint()); 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(ComputeV) svd.m_matrixV = qr.matrixQ();
if(ComputeU) svd.m_matrixU = qr.colsPermutation(); if(ComputeU) svd.m_matrixU = qr.colsPermutation();
return true; return true;

4
Eigen/src/Sparse/CholmodSupport.h
View File

@ -76,9 +76,9 @@ cholmod_sparse SparseMatrixBase<Derived>::asCholmodMatrix()
if (Derived::Flags & SelfAdjoint) if (Derived::Flags & SelfAdjoint)
{ {
if (Derived::Flags & UpperTriangular) if (Derived::Flags & Upper)
res.stype = 1; res.stype = 1;
else if (Derived::Flags & LowerTriangular) else if (Derived::Flags & Lower)
res.stype = -1; res.stype = -1;
else else
res.stype = 0; res.stype = 0;

8
Eigen/src/Sparse/RandomSetter.h
View File

@ -179,8 +179,8 @@ class RandomSetter
SwapStorage = 1 - MapTraits<ScalarWrapper>::IsSorted, SwapStorage = 1 - MapTraits<ScalarWrapper>::IsSorted,
TargetRowMajor = (SparseMatrixType::Flags & RowMajorBit) ? 1 : 0, TargetRowMajor = (SparseMatrixType::Flags & RowMajorBit) ? 1 : 0,
SetterRowMajor = SwapStorage ? 1-TargetRowMajor : TargetRowMajor, SetterRowMajor = SwapStorage ? 1-TargetRowMajor : TargetRowMajor,
IsUpperTriangular = SparseMatrixType::Flags & UpperTriangularBit, IsUpper = SparseMatrixType::Flags & Upper,
IsLowerTriangular = SparseMatrixType::Flags & LowerTriangularBit IsLower = SparseMatrixType::Flags & Lower
}; };
public: public:
@ -303,8 +303,8 @@ class RandomSetter
/** \returns a reference to the coefficient at given coordinates \a row, \a col */ /** \returns a reference to the coefficient at given coordinates \a row, \a col */
Scalar& operator() (int row, int col) Scalar& operator() (int row, int col)
{ {
ei_assert(((!IsUpperTriangular) || (row<=col)) && "Invalid access to an upper triangular matrix"); ei_assert(((!IsUpper) || (row<=col)) && "Invalid access to an upper triangular matrix");
ei_assert(((!IsLowerTriangular) || (col<=row)) && "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 outer = SetterRowMajor ? row : col;
const int inner = SetterRowMajor ? col : row; const int inner = SetterRowMajor ? col : row;
const int outerMajor = outer >> OuterPacketBits; // index of the packet/map const int outerMajor = outer >> OuterPacketBits; // index of the packet/map

4
Eigen/src/Sparse/SparseLDLT.h
View File

@ -333,12 +333,12 @@ bool SparseLDLT<MatrixType, Backend>::solveInPlace(MatrixBase<Derived> &b) const
return false; return false;
if (m_matrix.nonZeros()>0) // otherwise L==I 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); b = b.cwiseQuotient(m_diag);
// FIXME should be .adjoint() but it fails to compile... // FIXME should be .adjoint() but it fails to compile...
if (m_matrix.nonZeros()>0) // otherwise L==I 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; return true;
} }

8
Eigen/src/Sparse/SparseLLT.h
View File

@ -41,7 +41,7 @@ class SparseLLT
protected: protected:
typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
typedef SparseMatrix<Scalar,LowerTriangular> CholMatrixType; typedef SparseMatrix<Scalar> CholMatrixType;
enum { enum {
SupernodalFactorIsDirty = 0x10000, SupernodalFactorIsDirty = 0x10000,
@ -193,15 +193,15 @@ bool SparseLLT<MatrixType, Backend>::solveInPlace(MatrixBase<Derived> &b) const
const int size = m_matrix.rows(); const int size = m_matrix.rows();
ei_assert(size==b.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... // FIXME should be simply .adjoint() but it fails to compile...
if (NumTraits<Scalar>::IsComplex) if (NumTraits<Scalar>::IsComplex)
{ {
CholMatrixType aux = m_matrix.conjugate(); CholMatrixType aux = m_matrix.conjugate();
aux.transpose().template triangularView<UpperTriangular>().solveInPlace(b); aux.transpose().template triangularView<Upper>().solveInPlace(b);
} }
else else
m_matrix.transpose().template triangularView<UpperTriangular>().solveInPlace(b); m_matrix.transpose().template triangularView<Upper>().solveInPlace(b);
return true; return true;
} }

2
Eigen/src/Sparse/SparseLU.h
View File

@ -47,7 +47,7 @@ class SparseLU
protected: protected:
typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
typedef SparseMatrix<Scalar,LowerTriangular> LUMatrixType; typedef SparseMatrix<Scalar> LUMatrixType;
enum { enum {
MatrixLUIsDirty = 0x10000 MatrixLUIsDirty = 0x10000

4
Eigen/src/Sparse/SparseMatrixBase.h
View File

@ -532,8 +532,8 @@ template<typename Derived> class SparseMatrixBase : public AnyMatrixBase<Derived
// bool isIdentity(RealScalar prec = dummy_precision<Scalar>()) const; // bool isIdentity(RealScalar prec = dummy_precision<Scalar>()) const;
// bool isDiagonal(RealScalar prec = dummy_precision<Scalar>()) const; // bool isDiagonal(RealScalar prec = dummy_precision<Scalar>()) const;
// bool isUpperTriangular(RealScalar prec = dummy_precision<Scalar>()) const; // bool isUpper(RealScalar prec = dummy_precision<Scalar>()) const;
// bool isLowerTriangular(RealScalar prec = dummy_precision<Scalar>()) const; // bool isLower(RealScalar prec = dummy_precision<Scalar>()) const;
// template<typename OtherDerived> // template<typename OtherDerived>
// bool isOrthogonal(const MatrixBase<OtherDerived>& other, // bool isOrthogonal(const MatrixBase<OtherDerived>& other,

8
Eigen/src/Sparse/SparseSelfAdjointView.h
View File

@ -31,7 +31,7 @@
* \brief Pseudo expression to manipulate a triangular sparse matrix as a selfadjoint matrix. * \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 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 * 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() * with given dense storage of the coefficients. It is the return type of MatrixBase::selfadjointView()
@ -168,9 +168,9 @@ class SparseSelfAdjointTimeDenseProduct
enum { enum {
LhsIsRowMajor = (_Lhs::Flags&RowMajorBit)==RowMajorBit, LhsIsRowMajor = (_Lhs::Flags&RowMajorBit)==RowMajorBit,
ProcessFirstHalf = ProcessFirstHalf =
((UpLo&(UpperTriangularBit|LowerTriangularBit))==(UpperTriangularBit|LowerTriangularBit)) ((UpLo&(Upper|Lower))==(Upper|Lower))
|| ( (UpLo&UpperTriangularBit) && !LhsIsRowMajor) || ( (UpLo&Upper) && !LhsIsRowMajor)
|| ( (UpLo&LowerTriangularBit) && LhsIsRowMajor), || ( (UpLo&Lower) && LhsIsRowMajor),
ProcessSecondHalf = !ProcessFirstHalf ProcessSecondHalf = !ProcessFirstHalf
}; };
for (int j=0; j<m_lhs.outerSize(); ++j) for (int j=0; j<m_lhs.outerSize(); ++j)

4
Eigen/src/Sparse/SparseTriangularView.h
View File

@ -33,8 +33,8 @@ struct ei_traits<SparseTriangularView<MatrixType,Mode> >
template<typename MatrixType, int Mode> class SparseTriangularView template<typename MatrixType, int Mode> class SparseTriangularView
: public SparseMatrixBase<SparseTriangularView<MatrixType,Mode> > : public SparseMatrixBase<SparseTriangularView<MatrixType,Mode> >
{ {
enum { SkipFirst = (Mode==LowerTriangular && !(MatrixType::Flags&RowMajorBit)) enum { SkipFirst = (Mode==Lower && !(MatrixType::Flags&RowMajorBit))
|| (Mode==UpperTriangular && (MatrixType::Flags&RowMajorBit)) }; || (Mode==Upper && (MatrixType::Flags&RowMajorBit)) };
public: public:
class InnerIterator; class InnerIterator;

12
Eigen/src/Sparse/SuperLUSupport.h
View File

@ -193,9 +193,9 @@ struct SluMatrix : SuperMatrix
res.setScalarType<typename MatrixType::Scalar>(); res.setScalarType<typename MatrixType::Scalar>();
// FIXME the following is not very accurate // FIXME the following is not very accurate
if (MatrixType::Flags & UpperTriangular) if (MatrixType::Flags & Upper)
res.Mtype = SLU_TRU; res.Mtype = SLU_TRU;
if (MatrixType::Flags & LowerTriangular) if (MatrixType::Flags & Lower)
res.Mtype = SLU_TRL; res.Mtype = SLU_TRL;
if (MatrixType::Flags & SelfAdjoint) if (MatrixType::Flags & SelfAdjoint)
ei_assert(false && "SelfAdjoint matrix shape not supported by SuperLU"); ei_assert(false && "SelfAdjoint matrix shape not supported by SuperLU");
@ -251,9 +251,9 @@ struct SluMatrixMapHelper<SparseMatrixBase<Derived> >
res.setScalarType<typename MatrixType::Scalar>(); res.setScalarType<typename MatrixType::Scalar>();
// FIXME the following is not very accurate // FIXME the following is not very accurate
if (MatrixType::Flags & UpperTriangular) if (MatrixType::Flags & Upper)
res.Mtype = SLU_TRU; res.Mtype = SLU_TRU;
if (MatrixType::Flags & LowerTriangular) if (MatrixType::Flags & Lower)
res.Mtype = SLU_TRL; res.Mtype = SLU_TRL;
if (MatrixType::Flags & SelfAdjoint) if (MatrixType::Flags & SelfAdjoint)
ei_assert(false && "SelfAdjoint matrix shape not supported by SuperLU"); 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<Scalar,Dynamic,1> Vector;
typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType; typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType; typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
typedef SparseMatrix<Scalar,LowerTriangular|UnitDiagBit> LMatrixType; typedef SparseMatrix<Scalar,Lower|UnitDiag> LMatrixType;
typedef SparseMatrix<Scalar,UpperTriangular> UMatrixType; typedef SparseMatrix<Scalar,Upper> UMatrixType;
using Base::m_flags; using Base::m_flags;
using Base::m_status; using Base::m_status;

6
Eigen/src/Sparse/TaucsSupport.h
View File

@ -50,13 +50,13 @@ taucs_ccs_matrix SparseMatrixBase<Derived>::asTaucsMatrix()
ei_assert(false && "Scalar type not supported by TAUCS"); ei_assert(false && "Scalar type not supported by TAUCS");
} }
if (Flags & UpperTriangular) if (Flags & Upper)
res.flags |= TAUCS_UPPER; res.flags |= TAUCS_UPPER;
if (Flags & LowerTriangular) if (Flags & Lower)
res.flags |= TAUCS_LOWER; res.flags |= TAUCS_LOWER;
if (Flags & SelfAdjoint) if (Flags & SelfAdjoint)
res.flags |= (NumTraits<Scalar>::IsComplex ? TAUCS_HERMITIAN : TAUCS_SYMMETRIC); 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; res.flags |= TAUCS_TRIANGULAR;
return res; return res;

56
Eigen/src/Sparse/TriangularSolver.h
View File

@ -26,17 +26,17 @@
#define EIGEN_SPARSETRIANGULARSOLVER_H #define EIGEN_SPARSETRIANGULARSOLVER_H
template<typename Lhs, typename Rhs, int Mode, template<typename Lhs, typename Rhs, int Mode,
int UpLo = (Mode & LowerTriangularBit) int UpLo = (Mode & Lower)
? LowerTriangular ? Lower
: (Mode & UpperTriangularBit) : (Mode & Upper)
? UpperTriangular ? Upper
: -1, : -1,
int StorageOrder = int(ei_traits<Lhs>::Flags) & RowMajorBit> int StorageOrder = int(ei_traits<Lhs>::Flags) & RowMajorBit>
struct ei_sparse_solve_triangular_selector; struct ei_sparse_solve_triangular_selector;
// forward substitution, row-major // forward substitution, row-major
template<typename Lhs, typename Rhs, int Mode> 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; typedef typename Rhs::Scalar Scalar;
static void run(const Lhs& lhs, Rhs& other) static void run(const Lhs& lhs, Rhs& other)
@ -56,7 +56,7 @@ struct ei_sparse_solve_triangular_selector<Lhs,Rhs,Mode,LowerTriangular,RowMajor
break; break;
tmp -= lastVal * other.coeff(lastIndex,col); tmp -= lastVal * other.coeff(lastIndex,col);
} }
if (Mode & UnitDiagBit) if (Mode & UnitDiag)
other.coeffRef(i,col) = tmp; other.coeffRef(i,col) = tmp;
else else
{ {
@ -70,7 +70,7 @@ struct ei_sparse_solve_triangular_selector<Lhs,Rhs,Mode,LowerTriangular,RowMajor
// backward substitution, row-major // backward substitution, row-major
template<typename Lhs, typename Rhs, int Mode> 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; typedef typename Rhs::Scalar Scalar;
static void run(const Lhs& lhs, Rhs& other) 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); tmp -= it.value() * other.coeff(it.index(),col);
} }
if (Mode & UnitDiagBit) if (Mode & UnitDiag)
other.coeffRef(i,col) = tmp; other.coeffRef(i,col) = tmp;
else else
{ {
@ -103,7 +103,7 @@ struct ei_sparse_solve_triangular_selector<Lhs,Rhs,Mode,UpperTriangular,RowMajor
// forward substitution, col-major // forward substitution, col-major
template<typename Lhs, typename Rhs, int Mode> 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; typedef typename Rhs::Scalar Scalar;
static void run(const Lhs& lhs, Rhs& other) 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 if (tmp!=Scalar(0)) // optimization when other is actually sparse
{ {
typename Lhs::InnerIterator it(lhs, i); typename Lhs::InnerIterator it(lhs, i);
if(!(Mode & UnitDiagBit)) if(!(Mode & UnitDiag))
{ {
ei_assert(it.index()==i); ei_assert(it.index()==i);
tmp /= it.value(); tmp /= it.value();
@ -133,7 +133,7 @@ struct ei_sparse_solve_triangular_selector<Lhs,Rhs,Mode,LowerTriangular,ColMajor
// backward substitution, col-major // backward substitution, col-major
template<typename Lhs, typename Rhs, int Mode> 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; typedef typename Rhs::Scalar Scalar;
static void run(const Lhs& lhs, Rhs& other) 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); Scalar& tmp = other.coeffRef(i,col);
if (tmp!=Scalar(0)) // optimization when other is actually sparse 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 // FIXME lhs.coeff(i,i) might not be always efficient while it must simply be the
// last element of the column ! // 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() == m_matrix.rows());
ei_assert(m_matrix.cols() == other.rows()); ei_assert(m_matrix.cols() == other.rows());
ei_assert(!(Mode & ZeroDiagBit)); ei_assert(!(Mode & ZeroDiag));
ei_assert(Mode & (UpperTriangularBit|LowerTriangularBit)); ei_assert(Mode & (Upper|Lower));
enum { copy = ei_traits<OtherDerived>::Flags & RowMajorBit }; enum { copy = ei_traits<OtherDerived>::Flags & RowMajorBit };
@ -194,10 +194,10 @@ SparseTriangularView<ExpressionType,Mode>::solve(const MatrixBase<OtherDerived>&
// pure sparse path // pure sparse path
template<typename Lhs, typename Rhs, int Mode, template<typename Lhs, typename Rhs, int Mode,
int UpLo = (Mode & LowerTriangularBit) int UpLo = (Mode & Lower)
? LowerTriangular ? Lower
: (Mode & UpperTriangularBit) : (Mode & Upper)
? UpperTriangular ? Upper
: -1, : -1,
int StorageOrder = int(Lhs::Flags) & (RowMajorBit)> int StorageOrder = int(Lhs::Flags) & (RowMajorBit)>
struct ei_sparse_solve_triangular_sparse_selector; 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; typedef typename Rhs::Scalar Scalar;
static void run(const Lhs& lhs, Rhs& other) static void run(const Lhs& lhs, Rhs& other)
{ {
const bool IsLowerTriangular = (UpLo==LowerTriangular); const bool IsLower = (UpLo==Lower);
AmbiVector<Scalar> tempVector(other.rows()*2); AmbiVector<Scalar> tempVector(other.rows()*2);
tempVector.setBounds(0,other.rows()); 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(); tempVector.coeffRef(rhsIt.index()) = rhsIt.value();
} }
for(int i=IsLowerTriangular?0:lhs.cols()-1; for(int i=IsLower?0:lhs.cols()-1;
IsLowerTriangular?i<lhs.cols():i>=0; IsLower?i<lhs.cols():i>=0;
i+=IsLowerTriangular?1:-1) i+=IsLower?1:-1)
{ {
tempVector.restart(); tempVector.restart();
Scalar& ci = tempVector.coeffRef(i); Scalar& ci = tempVector.coeffRef(i);
@ -237,9 +237,9 @@ struct ei_sparse_solve_triangular_sparse_selector<Lhs,Rhs,Mode,UpLo,ColMajor>
{ {
// find // find
typename Lhs::InnerIterator it(lhs, i); typename Lhs::InnerIterator it(lhs, i);
if(!(Mode & UnitDiagBit)) if(!(Mode & UnitDiag))
{ {
if (IsLowerTriangular) if (IsLower)
{ {
ei_assert(it.index()==i); ei_assert(it.index()==i);
ci /= it.value(); ci /= it.value();
@ -248,7 +248,7 @@ struct ei_sparse_solve_triangular_sparse_selector<Lhs,Rhs,Mode,UpLo,ColMajor>
ci /= lhs.coeff(i,i); ci /= lhs.coeff(i,i);
} }
tempVector.restart(); tempVector.restart();
if (IsLowerTriangular) if (IsLower)
{ {
if (it.index()==i) if (it.index()==i)
++it; ++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() == m_matrix.rows());
ei_assert(m_matrix.cols() == other.rows()); ei_assert(m_matrix.cols() == other.rows());
ei_assert(!(Mode & ZeroDiagBit)); ei_assert(!(Mode & ZeroDiag));
ei_assert(Mode & (UpperTriangularBit|LowerTriangularBit)); ei_assert(Mode & (Upper|Lower));
// enum { copy = ei_traits<OtherDerived>::Flags & RowMajorBit }; // enum { copy = ei_traits<OtherDerived>::Flags & RowMajorBit };
@ -311,7 +311,7 @@ template<typename Derived>
template<typename OtherDerived> template<typename OtherDerived>
void SparseMatrixBase<Derived>::solveTriangularInPlace(MatrixBase<OtherDerived>& other) const 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 */ /** \deprecated */

4
Eigen/src/Sparse/UmfPackSupport.h
View File

@ -127,8 +127,8 @@ class SparseLU<MatrixType,UmfPack> : public SparseLU<MatrixType>
typedef Matrix<Scalar,Dynamic,1> Vector; typedef Matrix<Scalar,Dynamic,1> Vector;
typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType; typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType; typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
typedef SparseMatrix<Scalar,LowerTriangular|UnitDiagBit> LMatrixType; typedef SparseMatrix<Scalar,Lower|UnitDiagBit> LMatrixType;
typedef SparseMatrix<Scalar,UpperTriangular> UMatrixType; typedef SparseMatrix<Scalar,Upper> UMatrixType;
using Base::m_flags; using Base::m_flags;
using Base::m_status; using Base::m_status;

19
doc/A05_PortingFrom2To3.dox
View File

@ -29,13 +29,30 @@ vec.head(length)
vec.head<length>() vec.head<length>()
vec.tail(length) vec.tail(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 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>\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 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>\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 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>\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> </table>
\section CoefficientWiseOperations Coefficient wise operations \section CoefficientWiseOperations Coefficient wise operations

4
test/bandmatrix.cpp
View File

@ -64,8 +64,8 @@ template<typename MatrixType> void bandmatrix(const MatrixType& _m)
int a = std::max(0,cols-d-supers); int a = std::max(0,cols-d-supers);
int b = std::max(0,rows-d-subs); int b = std::max(0,rows-d-subs);
if(a>0) dm1.block(0,d+supers,rows,a).setZero(); 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(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<LowerTriangular>().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(); 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"; //std::cerr << m.m_data << "\n\n" << m.toDense() << "\n\n" << dm1 << "\n\n";
VERIFY_IS_APPROX(dm1,m.toDenseMatrix()); VERIFY_IS_APPROX(dm1,m.toDenseMatrix());

10
test/cholesky.cpp
View File

@ -58,12 +58,12 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
symm += a1 * a1.adjoint(); symm += a1 * a1.adjoint();
} }
SquareMatrixType symmUp = symm.template triangularView<UpperTriangular>(); SquareMatrixType symmUp = symm.template triangularView<Upper>();
SquareMatrixType symmLo = symm.template triangularView<LowerTriangular>(); SquareMatrixType symmLo = symm.template triangularView<Lower>();
// to test if really Cholesky only uses the upper triangular part, uncomment the following // to test if really Cholesky only uses the upper triangular part, uncomment the following
// FIXME: currently that fails !! // FIXME: currently that fails !!
//symm.template part<StrictlyLowerTriangular>().setZero(); //symm.template part<StrictlyLower>().setZero();
#ifdef HAS_GSL #ifdef HAS_GSL
// if (ei_is_same_type<RealScalar,double>::ret) // if (ei_is_same_type<RealScalar,double>::ret)
@ -94,7 +94,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
#endif #endif
{ {
LLT<SquareMatrixType,LowerTriangular> chollo(symmLo); LLT<SquareMatrixType,Lower> chollo(symmLo);
VERIFY_IS_APPROX(symm, chollo.matrixL().toDenseMatrix() * chollo.matrixL().adjoint().toDenseMatrix()); VERIFY_IS_APPROX(symm, chollo.matrixL().toDenseMatrix() * chollo.matrixL().adjoint().toDenseMatrix());
vecX = chollo.solve(vecB); vecX = chollo.solve(vecB);
VERIFY_IS_APPROX(symm * vecX, 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); VERIFY_IS_APPROX(symm * matX, matB);
// test the upper mode // 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()); VERIFY_IS_APPROX(symm, cholup.matrixL().toDenseMatrix() * cholup.matrixL().adjoint().toDenseMatrix());
vecX = cholup.solve(vecB); vecX = cholup.solve(vecB);
VERIFY_IS_APPROX(symm * vecX, vecB); VERIFY_IS_APPROX(symm * vecX, vecB);

32
test/product_notemporary.cpp
View File

@ -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 // 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.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( m3.noalias() -= (s1 * m1).template triangularView<Lower>() * 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<Upper>() * (m2+m2), 1);
VERIFY_EVALUATION_COUNT( rm3.noalias() = (s1 * m1.adjoint()).template triangularView<UnitUpperTriangular>() * m2.adjoint(), 0); 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.template triangularView<Lower>().solveInPlace(m3), 0);
VERIFY_EVALUATION_COUNT( m1.adjoint().template triangularView<LowerTriangular>().solveInPlace(m3.transpose()), 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() -= (s1 * m1).adjoint().template selfadjointView<Lower>() * (-m2*s3).adjoint(), 0);
VERIFY_EVALUATION_COUNT( m3.noalias() = s2 * m2.adjoint() * (s1 * m1.adjoint()).template selfadjointView<UpperTriangular>(), 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<LowerTriangular>() * m2.adjoint(), 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<Lower>() * (-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<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<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<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>() * 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); 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); 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() void test_product_notemporary()

24
test/product_selfadjoint.cpp
View File

@ -53,25 +53,25 @@ template<typename MatrixType> void product_selfadjoint(const MatrixType& m)
m1 = (m1.adjoint() + m1).eval(); m1 = (m1.adjoint() + m1).eval();
// rank2 update // rank2 update
m2 = m1.template triangularView<LowerTriangular>(); m2 = m1.template triangularView<Lower>();
m2.template selfadjointView<LowerTriangular>().rankUpdate(v1,v2); m2.template selfadjointView<Lower>().rankUpdate(v1,v2);
VERIFY_IS_APPROX(m2, (m1 + v1 * v2.adjoint()+ v2 * v1.adjoint()).template triangularView<LowerTriangular>().toDenseMatrix()); VERIFY_IS_APPROX(m2, (m1 + v1 * v2.adjoint()+ v2 * v1.adjoint()).template triangularView<Lower>().toDenseMatrix());
m2 = m1.template triangularView<UpperTriangular>(); m2 = m1.template triangularView<Upper>();
m2.template selfadjointView<UpperTriangular>().rankUpdate(-v1,s2*v2,s3); m2.template selfadjointView<Upper>().rankUpdate(-v1,s2*v2,s3);
VERIFY_IS_APPROX(m2, (m1 + (-s2*s3) * (v1 * v2.adjoint()+ v2 * v1.adjoint())).template triangularView<UpperTriangular>().toDenseMatrix()); VERIFY_IS_APPROX(m2, (m1 + (-s2*s3) * (v1 * v2.adjoint()+ v2 * v1.adjoint())).template triangularView<Upper>().toDenseMatrix());
m2 = m1.template triangularView<UpperTriangular>(); m2 = m1.template triangularView<Upper>();
m2.template selfadjointView<UpperTriangular>().rankUpdate(-r1.adjoint(),r2.adjoint()*s3,s1); 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<UpperTriangular>().toDenseMatrix()); VERIFY_IS_APPROX(m2, (m1 + (-s3*s1) * (r1.adjoint() * r2 + r2.adjoint() * r1)).template triangularView<Upper>().toDenseMatrix());
if (rows>1) if (rows>1)
{ {
m2 = m1.template triangularView<LowerTriangular>(); m2 = m1.template triangularView<Lower>();
m2.block(1,1,rows-1,cols-1).template selfadjointView<LowerTriangular>().rankUpdate(v1.tail(rows-1),v2.head(cols-1)); m2.block(1,1,rows-1,cols-1).template selfadjointView<Lower>().rankUpdate(v1.tail(rows-1),v2.head(cols-1));
m3 = m1; 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(); 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());
} }
} }

38
test/product_symm.cpp
View File

@ -28,10 +28,10 @@ template<int OtherSize> struct symm_extra {
template<typename M1, typename M2, typename Scalar> template<typename M1, typename M2, typename Scalar>
static void run(M1& m1, M1& m2, M2& rhs2, M2& rhs22, M2& rhs23, Scalar s1, Scalar s2) static void run(M1& m1, M1& m2, M2& rhs2, M2& rhs22, M2& rhs23, Scalar s1, Scalar s2)
{ {
m2 = m1.template triangularView<LowerTriangular>(); m2 = m1.template triangularView<Lower>();
VERIFY_IS_APPROX(rhs22 = (rhs2) * (m2).template selfadjointView<LowerTriangular>(), VERIFY_IS_APPROX(rhs22 = (rhs2) * (m2).template selfadjointView<Lower>(),
rhs23 = (rhs2) * (m1)); 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)); 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>(), Scalar s1 = ei_random<Scalar>(),
s2 = ei_random<Scalar>(); s2 = ei_random<Scalar>();
m2 = m1.template triangularView<LowerTriangular>(); m2 = m1.template triangularView<Lower>();
VERIFY_IS_APPROX(rhs12 = (s1*m2).template selfadjointView<LowerTriangular>() * (s2*rhs1), VERIFY_IS_APPROX(rhs12 = (s1*m2).template selfadjointView<Lower>() * (s2*rhs1),
rhs13 = (s1*m1) * (s2*rhs1)); rhs13 = (s1*m1) * (s2*rhs1));
m2 = m1.template triangularView<UpperTriangular>(); rhs12.setRandom(); rhs13 = rhs12; m2 = m1.template triangularView<Upper>(); rhs12.setRandom(); rhs13 = rhs12;
VERIFY_IS_APPROX(rhs12 += (s1*m2).template selfadjointView<UpperTriangular>() * (s2*rhs1), VERIFY_IS_APPROX(rhs12 += (s1*m2).template selfadjointView<Upper>() * (s2*rhs1),
rhs13 += (s1*m1) * (s2*rhs1)); rhs13 += (s1*m1) * (s2*rhs1));
m2 = m1.template triangularView<LowerTriangular>(); m2 = m1.template triangularView<Lower>();
VERIFY_IS_APPROX(rhs12 = (s1*m2).template selfadjointView<LowerTriangular>() * (s2*rhs2.adjoint()), VERIFY_IS_APPROX(rhs12 = (s1*m2).template selfadjointView<Lower>() * (s2*rhs2.adjoint()),
rhs13 = (s1*m1) * (s2*rhs2.adjoint())); rhs13 = (s1*m1) * (s2*rhs2.adjoint()));
m2 = m1.template triangularView<UpperTriangular>(); m2 = m1.template triangularView<Upper>();
VERIFY_IS_APPROX(rhs12 = (s1*m2).template selfadjointView<UpperTriangular>() * (s2*rhs2.adjoint()), VERIFY_IS_APPROX(rhs12 = (s1*m2).template selfadjointView<Upper>() * (s2*rhs2.adjoint()),
rhs13 = (s1*m1) * (s2*rhs2.adjoint())); rhs13 = (s1*m1) * (s2*rhs2.adjoint()));
m2 = m1.template triangularView<UpperTriangular>(); m2 = m1.template triangularView<Upper>();
VERIFY_IS_APPROX(rhs12 = (s1*m2.adjoint()).template selfadjointView<LowerTriangular>() * (s2*rhs2.adjoint()), VERIFY_IS_APPROX(rhs12 = (s1*m2.adjoint()).template selfadjointView<Lower>() * (s2*rhs2.adjoint()),
rhs13 = (s1*m1.adjoint()) * (s2*rhs2.adjoint())); rhs13 = (s1*m1.adjoint()) * (s2*rhs2.adjoint()));
// test row major = <...> // test row major = <...>
m2 = m1.template triangularView<LowerTriangular>(); rhs12.setRandom(); rhs13 = rhs12; m2 = m1.template triangularView<Lower>(); rhs12.setRandom(); rhs13 = rhs12;
VERIFY_IS_APPROX(rhs12 -= (s1*m2).template selfadjointView<LowerTriangular>() * (s2*rhs3), VERIFY_IS_APPROX(rhs12 -= (s1*m2).template selfadjointView<Lower>() * (s2*rhs3),
rhs13 -= (s1*m1) * (s2 * rhs3)); rhs13 -= (s1*m1) * (s2 * rhs3));
m2 = m1.template triangularView<UpperTriangular>(); m2 = m1.template triangularView<Upper>();
VERIFY_IS_APPROX(rhs12 = (s1*m2.adjoint()).template selfadjointView<LowerTriangular>() * (s2*rhs3).conjugate(), VERIFY_IS_APPROX(rhs12 = (s1*m2.adjoint()).template selfadjointView<Lower>() * (s2*rhs3).conjugate(),
rhs13 = (s1*m1.adjoint()) * (s2*rhs3).conjugate()); rhs13 = (s1*m1.adjoint()) * (s2*rhs3).conjugate());
m2 = m1.template triangularView<UpperTriangular>(); rhs13 = rhs12; m2 = m1.template triangularView<Upper>(); rhs13 = rhs12;
VERIFY_IS_APPROX(rhs12.noalias() += s1 * ((m2.adjoint()).template selfadjointView<LowerTriangular>() * (s2*rhs3).conjugate()), VERIFY_IS_APPROX(rhs12.noalias() += s1 * ((m2.adjoint()).template selfadjointView<Lower>() * (s2*rhs3).conjugate()),
rhs13 += (s1*m1.adjoint()) * (s2*rhs3).conjugate()); rhs13 += (s1*m1.adjoint()) * (s2*rhs3).conjugate());
// test matrix * selfadjoint // test matrix * selfadjoint

24
test/product_syrk.cpp
View File

@ -45,28 +45,28 @@ template<typename MatrixType> void syrk(const MatrixType& m)
Scalar s1 = ei_random<Scalar>(); Scalar s1 = ei_random<Scalar>();
m2.setZero(); m2.setZero();
VERIFY_IS_APPROX((m2.template selfadjointView<LowerTriangular>().rankUpdate(rhs2,s1)._expression()), VERIFY_IS_APPROX((m2.template selfadjointView<Lower>().rankUpdate(rhs2,s1)._expression()),
((s1 * rhs2 * rhs2.adjoint()).eval().template triangularView<LowerTriangular>().toDenseMatrix())); ((s1 * rhs2 * rhs2.adjoint()).eval().template triangularView<Lower>().toDenseMatrix()));
m2.setZero(); m2.setZero();
VERIFY_IS_APPROX(m2.template selfadjointView<UpperTriangular>().rankUpdate(rhs2,s1)._expression(), VERIFY_IS_APPROX(m2.template selfadjointView<Upper>().rankUpdate(rhs2,s1)._expression(),
(s1 * rhs2 * rhs2.adjoint()).eval().template triangularView<UpperTriangular>().toDenseMatrix()); (s1 * rhs2 * rhs2.adjoint()).eval().template triangularView<Upper>().toDenseMatrix());
m2.setZero(); m2.setZero();
VERIFY_IS_APPROX(m2.template selfadjointView<LowerTriangular>().rankUpdate(rhs1.adjoint(),s1)._expression(), VERIFY_IS_APPROX(m2.template selfadjointView<Lower>().rankUpdate(rhs1.adjoint(),s1)._expression(),
(s1 * rhs1.adjoint() * rhs1).eval().template triangularView<LowerTriangular>().toDenseMatrix()); (s1 * rhs1.adjoint() * rhs1).eval().template triangularView<Lower>().toDenseMatrix());
m2.setZero(); m2.setZero();
VERIFY_IS_APPROX(m2.template selfadjointView<UpperTriangular>().rankUpdate(rhs1.adjoint(),s1)._expression(), VERIFY_IS_APPROX(m2.template selfadjointView<Upper>().rankUpdate(rhs1.adjoint(),s1)._expression(),
(s1 * rhs1.adjoint() * rhs1).eval().template triangularView<UpperTriangular>().toDenseMatrix()); (s1 * rhs1.adjoint() * rhs1).eval().template triangularView<Upper>().toDenseMatrix());
m2.setZero(); m2.setZero();
VERIFY_IS_APPROX(m2.template selfadjointView<LowerTriangular>().rankUpdate(rhs3.adjoint(),s1)._expression(), VERIFY_IS_APPROX(m2.template selfadjointView<Lower>().rankUpdate(rhs3.adjoint(),s1)._expression(),
(s1 * rhs3.adjoint() * rhs3).eval().template triangularView<LowerTriangular>().toDenseMatrix()); (s1 * rhs3.adjoint() * rhs3).eval().template triangularView<Lower>().toDenseMatrix());
m2.setZero(); m2.setZero();
VERIFY_IS_APPROX(m2.template selfadjointView<UpperTriangular>().rankUpdate(rhs3.adjoint(),s1)._expression(), VERIFY_IS_APPROX(m2.template selfadjointView<Upper>().rankUpdate(rhs3.adjoint(),s1)._expression(),
(s1 * rhs3.adjoint() * rhs3).eval().template triangularView<UpperTriangular>().toDenseMatrix()); (s1 * rhs3.adjoint() * rhs3).eval().template triangularView<Upper>().toDenseMatrix());
} }
void test_product_syrk() void test_product_syrk()

36
test/product_trmm.cpp
View File

@ -36,27 +36,27 @@ template<typename Scalar> void trmm(int size,int othersize)
s2 = ei_random<Scalar>(); s2 = ei_random<Scalar>();
tri.setRandom(); tri.setRandom();
loTri = tri.template triangularView<LowerTriangular>(); loTri = tri.template triangularView<Lower>();
upTri = tri.template triangularView<UpperTriangular>(); upTri = tri.template triangularView<Upper>();
ge1.setRandom(); ge1.setRandom();
ge2.setRandom(); ge2.setRandom();
VERIFY_IS_APPROX( ge3 = tri.template triangularView<LowerTriangular>() * ge1, loTri * ge1); VERIFY_IS_APPROX( ge3 = tri.template triangularView<Lower>() * ge1, loTri * ge1);
VERIFY_IS_APPROX(rge3 = tri.template triangularView<LowerTriangular>() * ge1, loTri * ge1); VERIFY_IS_APPROX(rge3 = tri.template triangularView<Lower>() * ge1, loTri * ge1);
VERIFY_IS_APPROX( ge3 = tri.template triangularView<UpperTriangular>() * ge1, upTri * ge1); VERIFY_IS_APPROX( ge3 = tri.template triangularView<Upper>() * ge1, upTri * ge1);
VERIFY_IS_APPROX(rge3 = tri.template triangularView<UpperTriangular>() * ge1, upTri * ge1); VERIFY_IS_APPROX(rge3 = tri.template triangularView<Upper>() * ge1, upTri * ge1);
VERIFY_IS_APPROX( ge3 = (s1*tri.adjoint()).template triangularView<UpperTriangular>() * (s2*ge1), s1*loTri.adjoint() * (s2*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<UpperTriangular>() * ge1, loTri.adjoint() * ge1); VERIFY_IS_APPROX(rge3 = tri.adjoint().template triangularView<Upper>() * ge1, loTri.adjoint() * ge1);
VERIFY_IS_APPROX( ge3 = tri.adjoint().template triangularView<LowerTriangular>() * ge1, upTri.adjoint() * ge1); VERIFY_IS_APPROX( ge3 = tri.adjoint().template triangularView<Lower>() * ge1, upTri.adjoint() * ge1);
VERIFY_IS_APPROX(rge3 = tri.adjoint().template triangularView<LowerTriangular>() * ge1, upTri.adjoint() * ge1); VERIFY_IS_APPROX(rge3 = tri.adjoint().template triangularView<Lower>() * ge1, upTri.adjoint() * ge1);
VERIFY_IS_APPROX( ge3 = tri.template triangularView<LowerTriangular>() * ge2.adjoint(), loTri * ge2.adjoint()); VERIFY_IS_APPROX( ge3 = tri.template triangularView<Lower>() * ge2.adjoint(), loTri * ge2.adjoint());
VERIFY_IS_APPROX(rge3 = tri.template triangularView<LowerTriangular>() * 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<UpperTriangular>() * ge2.adjoint(), upTri * ge2.adjoint()); VERIFY_IS_APPROX( ge3 = tri.template triangularView<Upper>() * ge2.adjoint(), upTri * ge2.adjoint());
VERIFY_IS_APPROX(rge3 = tri.template triangularView<UpperTriangular>() * 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<UpperTriangular>() * ge2.adjoint(), ei_conj(s1) * loTri.adjoint() * 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<UpperTriangular>() * ge2.adjoint(), 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<LowerTriangular>() * ge2.adjoint(), upTri.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<LowerTriangular>() * 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() void test_product_trmm()

50
test/product_trmv.cpp
View File

@ -44,37 +44,37 @@ template<typename MatrixType> void trmv(const MatrixType& m)
m1 = MatrixType::Random(rows, cols); m1 = MatrixType::Random(rows, cols);
// check with a column-major matrix // check with a column-major matrix
m3 = m1.template triangularView<Eigen::LowerTriangular>(); m3 = m1.template triangularView<Eigen::Lower>();
VERIFY((m3 * v1).isApprox(m1.template triangularView<Eigen::LowerTriangular>() * v1, largerEps)); VERIFY((m3 * v1).isApprox(m1.template triangularView<Eigen::Lower>() * v1, largerEps));
m3 = m1.template triangularView<Eigen::UpperTriangular>(); m3 = m1.template triangularView<Eigen::Upper>();
VERIFY((m3 * v1).isApprox(m1.template triangularView<Eigen::UpperTriangular>() * v1, largerEps)); VERIFY((m3 * v1).isApprox(m1.template triangularView<Eigen::Upper>() * v1, largerEps));
m3 = m1.template triangularView<Eigen::UnitLowerTriangular>(); m3 = m1.template triangularView<Eigen::UnitLower>();
VERIFY((m3 * v1).isApprox(m1.template triangularView<Eigen::UnitLowerTriangular>() * v1, largerEps)); VERIFY((m3 * v1).isApprox(m1.template triangularView<Eigen::UnitLower>() * v1, largerEps));
m3 = m1.template triangularView<Eigen::UnitUpperTriangular>(); m3 = m1.template triangularView<Eigen::UnitUpper>();
VERIFY((m3 * v1).isApprox(m1.template triangularView<Eigen::UnitUpperTriangular>() * v1, largerEps)); VERIFY((m3 * v1).isApprox(m1.template triangularView<Eigen::UnitUpper>() * v1, largerEps));
// check conjugated and scalar multiple expressions (col-major) // check conjugated and scalar multiple expressions (col-major)
m3 = m1.template triangularView<Eigen::LowerTriangular>(); m3 = m1.template triangularView<Eigen::Lower>();
VERIFY(((s1*m3).conjugate() * v1).isApprox((s1*m1).conjugate().template triangularView<Eigen::LowerTriangular>() * v1, largerEps)); VERIFY(((s1*m3).conjugate() * v1).isApprox((s1*m1).conjugate().template triangularView<Eigen::Lower>() * v1, largerEps));
m3 = m1.template triangularView<Eigen::UpperTriangular>(); m3 = m1.template triangularView<Eigen::Upper>();
VERIFY((m3.conjugate() * v1.conjugate()).isApprox(m1.conjugate().template triangularView<Eigen::UpperTriangular>() * v1.conjugate(), largerEps)); VERIFY((m3.conjugate() * v1.conjugate()).isApprox(m1.conjugate().template triangularView<Eigen::Upper>() * v1.conjugate(), largerEps));
// check with a row-major matrix // check with a row-major matrix
m3 = m1.template triangularView<Eigen::UpperTriangular>(); m3 = m1.template triangularView<Eigen::Upper>();
VERIFY((m3.transpose() * v1).isApprox(m1.transpose().template triangularView<Eigen::LowerTriangular>() * v1, largerEps)); VERIFY((m3.transpose() * v1).isApprox(m1.transpose().template triangularView<Eigen::Lower>() * v1, largerEps));
m3 = m1.template triangularView<Eigen::LowerTriangular>(); m3 = m1.template triangularView<Eigen::Lower>();
VERIFY((m3.transpose() * v1).isApprox(m1.transpose().template triangularView<Eigen::UpperTriangular>() * v1, largerEps)); VERIFY((m3.transpose() * v1).isApprox(m1.transpose().template triangularView<Eigen::Upper>() * v1, largerEps));
m3 = m1.template triangularView<Eigen::UnitUpperTriangular>(); m3 = m1.template triangularView<Eigen::UnitUpper>();
VERIFY((m3.transpose() * v1).isApprox(m1.transpose().template triangularView<Eigen::UnitLowerTriangular>() * v1, largerEps)); VERIFY((m3.transpose() * v1).isApprox(m1.transpose().template triangularView<Eigen::UnitLower>() * v1, largerEps));
m3 = m1.template triangularView<Eigen::UnitLowerTriangular>(); m3 = m1.template triangularView<Eigen::UnitLower>();
VERIFY((m3.transpose() * v1).isApprox(m1.transpose().template triangularView<Eigen::UnitUpperTriangular>() * v1, largerEps)); VERIFY((m3.transpose() * v1).isApprox(m1.transpose().template triangularView<Eigen::UnitUpper>() * v1, largerEps));
// check conjugated and scalar multiple expressions (row-major) // check conjugated and scalar multiple expressions (row-major)
m3 = m1.template triangularView<Eigen::UpperTriangular>(); m3 = m1.template triangularView<Eigen::Upper>();
VERIFY((m3.adjoint() * v1).isApprox(m1.adjoint().template triangularView<Eigen::LowerTriangular>() * v1, largerEps)); VERIFY((m3.adjoint() * v1).isApprox(m1.adjoint().template triangularView<Eigen::Lower>() * v1, largerEps));
m3 = m1.template triangularView<Eigen::LowerTriangular>(); m3 = m1.template triangularView<Eigen::Lower>();
VERIFY((m3.adjoint() * (s1*v1.conjugate())).isApprox(m1.adjoint().template triangularView<Eigen::UpperTriangular>() * (s1*v1.conjugate()), largerEps)); VERIFY((m3.adjoint() * (s1*v1.conjugate())).isApprox(m1.adjoint().template triangularView<Eigen::Upper>() * (s1*v1.conjugate()), largerEps));
m3 = m1.template triangularView<Eigen::UnitUpperTriangular>(); m3 = m1.template triangularView<Eigen::UnitUpper>();
// TODO check with sub-matrices // TODO check with sub-matrices
} }

32
test/product_trsolve.cpp
View File

@ -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); 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); 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.conjugate().template triangularView<Lower>(), cmRhs);
VERIFY_TRSM(cmLhs .template triangularView<UpperTriangular>(), cmRhs); VERIFY_TRSM(cmLhs .template triangularView<Upper>(), cmRhs);
VERIFY_TRSM(cmLhs .template triangularView<LowerTriangular>(), rmRhs); VERIFY_TRSM(cmLhs .template triangularView<Lower>(), rmRhs);
VERIFY_TRSM(cmLhs.conjugate().template triangularView<UpperTriangular>(), rmRhs); VERIFY_TRSM(cmLhs.conjugate().template triangularView<Upper>(), rmRhs);
VERIFY_TRSM(cmLhs.conjugate().template triangularView<UnitLowerTriangular>(), cmRhs); VERIFY_TRSM(cmLhs.conjugate().template triangularView<UnitLower>(), cmRhs);
VERIFY_TRSM(cmLhs .template triangularView<UnitUpperTriangular>(), rmRhs); VERIFY_TRSM(cmLhs .template triangularView<UnitUpper>(), rmRhs);
VERIFY_TRSM(rmLhs .template triangularView<LowerTriangular>(), cmRhs); VERIFY_TRSM(rmLhs .template triangularView<Lower>(), cmRhs);
VERIFY_TRSM(rmLhs.conjugate().template triangularView<UnitUpperTriangular>(), rmRhs); VERIFY_TRSM(rmLhs.conjugate().template triangularView<UnitUpper>(), rmRhs);
VERIFY_TRSM_ONTHERIGHT(cmLhs.conjugate().template triangularView<LowerTriangular>(), cmRhs); VERIFY_TRSM_ONTHERIGHT(cmLhs.conjugate().template triangularView<Lower>(), cmRhs);
VERIFY_TRSM_ONTHERIGHT(cmLhs .template triangularView<UpperTriangular>(), cmRhs); VERIFY_TRSM_ONTHERIGHT(cmLhs .template triangularView<Upper>(), cmRhs);
VERIFY_TRSM_ONTHERIGHT(cmLhs .template triangularView<LowerTriangular>(), rmRhs); VERIFY_TRSM_ONTHERIGHT(cmLhs .template triangularView<Lower>(), rmRhs);
VERIFY_TRSM_ONTHERIGHT(cmLhs.conjugate().template triangularView<UpperTriangular>(), rmRhs); VERIFY_TRSM_ONTHERIGHT(cmLhs.conjugate().template triangularView<Upper>(), rmRhs);
VERIFY_TRSM_ONTHERIGHT(cmLhs.conjugate().template triangularView<UnitLowerTriangular>(), cmRhs); VERIFY_TRSM_ONTHERIGHT(cmLhs.conjugate().template triangularView<UnitLower>(), cmRhs);
VERIFY_TRSM_ONTHERIGHT(cmLhs .template triangularView<UnitUpperTriangular>(), rmRhs); VERIFY_TRSM_ONTHERIGHT(cmLhs .template triangularView<UnitUpper>(), rmRhs);
VERIFY_TRSM_ONTHERIGHT(rmLhs .template triangularView<LowerTriangular>(), cmRhs); VERIFY_TRSM_ONTHERIGHT(rmLhs .template triangularView<Lower>(), cmRhs);
VERIFY_TRSM_ONTHERIGHT(rmLhs.conjugate().template triangularView<UnitUpperTriangular>(), rmRhs); VERIFY_TRSM_ONTHERIGHT(rmLhs.conjugate().template triangularView<UnitUpper>(), rmRhs);
} }
void test_product_trsolve() void test_product_trsolve()

4
test/qr_colpivoting.cpp
View File

@ -47,7 +47,7 @@ template<typename MatrixType> void qr()
MatrixQType q = qr.householderQ(); MatrixQType q = qr.householderQ();
VERIFY_IS_UNITARY(q); 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(); MatrixType c = q * r * qr.colsPermutation().inverse();
VERIFY_IS_APPROX(m1, c); VERIFY_IS_APPROX(m1, c);
@ -72,7 +72,7 @@ template<typename MatrixType, int Cols2> void qr_fixedsize()
VERIFY(!qr.isInvertible()); VERIFY(!qr.isInvertible());
VERIFY(!qr.isSurjective()); 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(); Matrix<Scalar,Rows,Cols> c = qr.householderQ() * r * qr.colsPermutation().inverse();
VERIFY_IS_APPROX(m1, c); VERIFY_IS_APPROX(m1, c);

6
test/sparse_product.cpp
View File

@ -115,9 +115,9 @@ template<typename SparseMatrixType> void sparse_product(const SparseMatrixType&
VERIFY_IS_APPROX(mS, refS); VERIFY_IS_APPROX(mS, refS);
VERIFY_IS_APPROX(x=mS*b, refX=refS*b); 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=mUp.template selfadjointView<Upper>()*b, refX=refS*b);
VERIFY_IS_APPROX(x=mLo.template selfadjointView<LowerTriangular>()*b, refX=refS*b); VERIFY_IS_APPROX(x=mLo.template selfadjointView<Lower>()*b, refX=refS*b);
VERIFY_IS_APPROX(x=mS.template selfadjointView<UpperTriangular|LowerTriangular>()*b, refX=refS*b); VERIFY_IS_APPROX(x=mS.template selfadjointView<Upper|Lower>()*b, refX=refS*b);
} }
} }

41
test/sparse_solvers.cpp
View File

@ -65,23 +65,23 @@ template<typename Scalar> void sparse_solvers(int rows, int cols)
// lower - dense // lower - dense
initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords); initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords);
VERIFY_IS_APPROX(refMat2.template triangularView<LowerTriangular>().solve(vec2), VERIFY_IS_APPROX(refMat2.template triangularView<Lower>().solve(vec2),
m2.template triangularView<LowerTriangular>().solve(vec3)); m2.template triangularView<Lower>().solve(vec3));
// upper - dense // upper - dense
initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, &zeroCoords, &nonzeroCoords); initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, &zeroCoords, &nonzeroCoords);
VERIFY_IS_APPROX(refMat2.template triangularView<UpperTriangular>().solve(vec2), VERIFY_IS_APPROX(refMat2.template triangularView<Upper>().solve(vec2),
m2.template triangularView<UpperTriangular>().solve(vec3)); m2.template triangularView<Upper>().solve(vec3));
// lower - transpose // lower - transpose
initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords); initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords);
VERIFY_IS_APPROX(refMat2.transpose().template triangularView<UpperTriangular>().solve(vec2), VERIFY_IS_APPROX(refMat2.transpose().template triangularView<Upper>().solve(vec2),
m2.transpose().template triangularView<UpperTriangular>().solve(vec3)); m2.transpose().template triangularView<Upper>().solve(vec3));
// upper - transpose // upper - transpose
initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, &zeroCoords, &nonzeroCoords); initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, &zeroCoords, &nonzeroCoords);
VERIFY_IS_APPROX(refMat2.transpose().template triangularView<LowerTriangular>().solve(vec2), VERIFY_IS_APPROX(refMat2.transpose().template triangularView<Lower>().solve(vec2),
m2.transpose().template triangularView<LowerTriangular>().solve(vec3)); m2.transpose().template triangularView<Lower>().solve(vec3));
SparseMatrix<Scalar> matB(rows, rows); SparseMatrix<Scalar> matB(rows, rows);
DenseMatrix refMatB = DenseMatrix::Zero(rows, rows); DenseMatrix refMatB = DenseMatrix::Zero(rows, rows);
@ -89,21 +89,21 @@ template<typename Scalar> void sparse_solvers(int rows, int cols)
// lower - sparse // lower - sparse
initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular); initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular);
initSparse<Scalar>(density, refMatB, matB); initSparse<Scalar>(density, refMatB, matB);
refMat2.template triangularView<LowerTriangular>().solveInPlace(refMatB); refMat2.template triangularView<Lower>().solveInPlace(refMatB);
m2.template triangularView<LowerTriangular>().solveInPlace(matB); m2.template triangularView<Lower>().solveInPlace(matB);
VERIFY_IS_APPROX(matB.toDense(), refMatB); VERIFY_IS_APPROX(matB.toDense(), refMatB);
// upper - sparse // upper - sparse
initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular); initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular);
initSparse<Scalar>(density, refMatB, matB); initSparse<Scalar>(density, refMatB, matB);
refMat2.template triangularView<UpperTriangular>().solveInPlace(refMatB); refMat2.template triangularView<Upper>().solveInPlace(refMatB);
m2.template triangularView<UpperTriangular>().solveInPlace(matB); m2.template triangularView<Upper>().solveInPlace(matB);
VERIFY_IS_APPROX(matB, refMatB); VERIFY_IS_APPROX(matB, refMatB);
// test deprecated API // test deprecated API
initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords); initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords);
VERIFY_IS_APPROX(refMat2.template triangularView<LowerTriangular>().solve(vec2), VERIFY_IS_APPROX(refMat2.template triangularView<Lower>().solve(vec2),
m2.template triangularView<LowerTriangular>().solve(vec3)); m2.template triangularView<Lower>().solve(vec3));
} }
// test LLT // test LLT
@ -118,29 +118,28 @@ template<typename Scalar> void sparse_solvers(int rows, int cols)
initSPD(density, refMat2, m2); initSPD(density, refMat2, m2);
refX = refMat2.llt().solve(b); refX = refMat2.llt().solve(b);
typedef SparseMatrix<Scalar,LowerTriangular|SelfAdjoint> SparseSelfAdjointMatrix;
if (!NumTraits<Scalar>::IsComplex) if (!NumTraits<Scalar>::IsComplex)
{ {
x = b; x = b;
SparseLLT<SparseSelfAdjointMatrix> (m2).solveInPlace(x); SparseLLT<SparseMatrix<Scalar> > (m2).solveInPlace(x);
VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: default"); VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: default");
} }
#ifdef EIGEN_CHOLMOD_SUPPORT #ifdef EIGEN_CHOLMOD_SUPPORT
x = b; x = b;
SparseLLT<SparseSelfAdjointMatrix,Cholmod>(m2).solveInPlace(x); SparseLLT<SparseMatrix<Scalar> ,Cholmod>(m2).solveInPlace(x);
VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: cholmod"); VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: cholmod");
#endif #endif
#ifdef EIGEN_TAUCS_SUPPORT #ifdef EIGEN_TAUCS_SUPPORT
x = b; 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)"); VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (IncompleteFactorization)");
// TODO fix TAUCS with complexes // TODO fix TAUCS with complexes
x = b; 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)"); VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (SupernodalMultifrontal)");
x = b; 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)"); VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (SupernodalLeftLooking)");
#endif #endif
} }
@ -161,7 +160,7 @@ template<typename Scalar> void sparse_solvers(int rows, int cols)
refMat2.diagonal() *= 0.5; refMat2.diagonal() *= 0.5;
refX = refMat2.llt().solve(b); // FIXME use LLT to compute the reference because LDLT seems to fail with large matrices 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; x = b;
SparseLDLT<SparseSelfAdjointMatrix> ldlt(m2); SparseLDLT<SparseSelfAdjointMatrix> ldlt(m2);
if (ldlt.succeeded()) if (ldlt.succeeded())

130
test/triangular.cpp
View File

@ -53,14 +53,14 @@ template<typename MatrixType> void triangular_square(const MatrixType& m)
v2 = VectorType::Random(rows), v2 = VectorType::Random(rows),
vzero = VectorType::Zero(rows); vzero = VectorType::Zero(rows);
MatrixType m1up = m1.template triangularView<UpperTriangular>(); MatrixType m1up = m1.template triangularView<Upper>();
MatrixType m2up = m2.template triangularView<UpperTriangular>(); MatrixType m2up = m2.template triangularView<Upper>();
if (rows*cols>1) if (rows*cols>1)
{ {
VERIFY(m1up.isUpperTriangular()); VERIFY(m1up.isUpper());
VERIFY(m2up.transpose().isLowerTriangular()); VERIFY(m2up.transpose().isLower());
VERIFY(!m2.isLowerTriangular()); VERIFY(!m2.isLower());
} }
// VERIFY_IS_APPROX(m1up.transpose() * m2, m1.upper().transpose().lower() * m2); // 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+= // test overloaded operator+=
r1.setZero(); r1.setZero();
r2.setZero(); r2.setZero();
r1.template triangularView<UpperTriangular>() += m1; r1.template triangularView<Upper>() += m1;
r2 += m1up; r2 += m1up;
VERIFY_IS_APPROX(r1,r2); VERIFY_IS_APPROX(r1,r2);
// test overloaded operator= // test overloaded operator=
m1.setZero(); m1.setZero();
m1.template triangularView<UpperTriangular>() = m2.transpose() + m2; m1.template triangularView<Upper>() = m2.transpose() + m2;
m3 = 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= // test overloaded operator=
m1.setZero(); m1.setZero();
m1.template triangularView<LowerTriangular>() = m2.transpose() + m2; m1.template triangularView<Lower>() = m2.transpose() + m2;
VERIFY_IS_APPROX(m3.template triangularView<LowerTriangular>().toDenseMatrix(), m1); VERIFY_IS_APPROX(m3.template triangularView<Lower>().toDenseMatrix(), m1);
m1 = MatrixType::Random(rows, cols); m1 = MatrixType::Random(rows, cols);
for (int i=0; i<rows; ++i) for (int i=0; i<rows; ++i)
@ -89,49 +89,49 @@ template<typename MatrixType> void triangular_square(const MatrixType& m)
Transpose<MatrixType> trm4(m4); Transpose<MatrixType> trm4(m4);
// test back and forward subsitution with a vector as the rhs // test back and forward subsitution with a vector as the rhs
m3 = m1.template triangularView<UpperTriangular>(); m3 = m1.template triangularView<Upper>();
VERIFY(v2.isApprox(m3.adjoint() * (m1.adjoint().template triangularView<LowerTriangular>().solve(v2)), largerEps)); VERIFY(v2.isApprox(m3.adjoint() * (m1.adjoint().template triangularView<Lower>().solve(v2)), largerEps));
m3 = m1.template triangularView<LowerTriangular>(); m3 = m1.template triangularView<Lower>();
VERIFY(v2.isApprox(m3.transpose() * (m1.transpose().template triangularView<UpperTriangular>().solve(v2)), largerEps)); VERIFY(v2.isApprox(m3.transpose() * (m1.transpose().template triangularView<Upper>().solve(v2)), largerEps));
m3 = m1.template triangularView<UpperTriangular>(); m3 = m1.template triangularView<Upper>();
VERIFY(v2.isApprox(m3 * (m1.template triangularView<UpperTriangular>().solve(v2)), largerEps)); VERIFY(v2.isApprox(m3 * (m1.template triangularView<Upper>().solve(v2)), largerEps));
m3 = m1.template triangularView<LowerTriangular>(); m3 = m1.template triangularView<Lower>();
VERIFY(v2.isApprox(m3.conjugate() * (m1.conjugate().template triangularView<LowerTriangular>().solve(v2)), largerEps)); VERIFY(v2.isApprox(m3.conjugate() * (m1.conjugate().template triangularView<Lower>().solve(v2)), largerEps));
// test back and forward subsitution with a matrix as the rhs // test back and forward subsitution with a matrix as the rhs
m3 = m1.template triangularView<UpperTriangular>(); m3 = m1.template triangularView<Upper>();
VERIFY(m2.isApprox(m3.adjoint() * (m1.adjoint().template triangularView<LowerTriangular>().solve(m2)), largerEps)); VERIFY(m2.isApprox(m3.adjoint() * (m1.adjoint().template triangularView<Lower>().solve(m2)), largerEps));
m3 = m1.template triangularView<LowerTriangular>(); m3 = m1.template triangularView<Lower>();
VERIFY(m2.isApprox(m3.transpose() * (m1.transpose().template triangularView<UpperTriangular>().solve(m2)), largerEps)); VERIFY(m2.isApprox(m3.transpose() * (m1.transpose().template triangularView<Upper>().solve(m2)), largerEps));
m3 = m1.template triangularView<UpperTriangular>(); m3 = m1.template triangularView<Upper>();
VERIFY(m2.isApprox(m3 * (m1.template triangularView<UpperTriangular>().solve(m2)), largerEps)); VERIFY(m2.isApprox(m3 * (m1.template triangularView<Upper>().solve(m2)), largerEps));
m3 = m1.template triangularView<LowerTriangular>(); m3 = m1.template triangularView<Lower>();
VERIFY(m2.isApprox(m3.conjugate() * (m1.conjugate().template triangularView<LowerTriangular>().solve(m2)), largerEps)); VERIFY(m2.isApprox(m3.conjugate() * (m1.conjugate().template triangularView<Lower>().solve(m2)), largerEps));
// check M * inv(L) using in place API // check M * inv(L) using in place API
m4 = m3; 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>())); VERIFY(m4.cwiseAbs().isIdentity(test_precision<RealScalar>()));
// check M * inv(U) using in place API // check M * inv(U) using in place API
m3 = m1.template triangularView<UpperTriangular>(); m3 = m1.template triangularView<Upper>();
m4 = m3; 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>())); VERIFY(m4.cwiseAbs().isIdentity(test_precision<RealScalar>()));
// check solve with unit diagonal // check solve with unit diagonal
m3 = m1.template triangularView<UnitUpperTriangular>(); m3 = m1.template triangularView<UnitUpper>();
VERIFY(m2.isApprox(m3 * (m1.template triangularView<UnitUpperTriangular>().solve(m2)), largerEps)); VERIFY(m2.isApprox(m3 * (m1.template triangularView<UnitUpper>().solve(m2)), largerEps));
// VERIFY(( m1.template triangularView<UpperTriangular>() // VERIFY(( m1.template triangularView<Upper>()
// * m2.template triangularView<UpperTriangular>()).isUpperTriangular()); // * m2.template triangularView<Upper>()).isUpper());
// test swap // test swap
m1.setOnes(); m1.setOnes();
m2.setZero(); m2.setZero();
m2.template triangularView<UpperTriangular>().swap(m1); m2.template triangularView<Upper>().swap(m1);
m3.setZero(); m3.setZero();
m3.template triangularView<UpperTriangular>().setOnes(); m3.template triangularView<Upper>().setOnes();
VERIFY_IS_APPROX(m2,m3); VERIFY_IS_APPROX(m2,m3);
} }
@ -165,69 +165,69 @@ template<typename MatrixType> void triangular_rect(const MatrixType& m)
v2 = VectorType::Random(rows), v2 = VectorType::Random(rows),
vzero = VectorType::Zero(rows); vzero = VectorType::Zero(rows);
MatrixType m1up = m1.template triangularView<UpperTriangular>(); MatrixType m1up = m1.template triangularView<Upper>();
MatrixType m2up = m2.template triangularView<UpperTriangular>(); MatrixType m2up = m2.template triangularView<Upper>();
if (rows*cols>1) if (rows*cols>1)
{ {
VERIFY(m1up.isUpperTriangular()); VERIFY(m1up.isUpper());
VERIFY(m2up.transpose().isLowerTriangular()); VERIFY(m2up.transpose().isLower());
VERIFY(!m2.isLowerTriangular()); VERIFY(!m2.isLower());
} }
// test overloaded operator+= // test overloaded operator+=
r1.setZero(); r1.setZero();
r2.setZero(); r2.setZero();
r1.template triangularView<UpperTriangular>() += m1; r1.template triangularView<Upper>() += m1;
r2 += m1up; r2 += m1up;
VERIFY_IS_APPROX(r1,r2); VERIFY_IS_APPROX(r1,r2);
// test overloaded operator= // test overloaded operator=
m1.setZero(); m1.setZero();
m1.template triangularView<UpperTriangular>() = 3 * m2; m1.template triangularView<Upper>() = 3 * m2;
m3 = 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.setZero();
m1.template triangularView<LowerTriangular>() = 3 * m2; m1.template triangularView<Lower>() = 3 * m2;
VERIFY_IS_APPROX(m3.template triangularView<LowerTriangular>().toDenseMatrix(), m1); VERIFY_IS_APPROX(m3.template triangularView<Lower>().toDenseMatrix(), m1);
m1.setZero(); m1.setZero();
m1.template triangularView<StrictlyUpperTriangular>() = 3 * m2; m1.template triangularView<StrictlyUpper>() = 3 * m2;
VERIFY_IS_APPROX(m3.template triangularView<StrictlyUpperTriangular>().toDenseMatrix(), m1); VERIFY_IS_APPROX(m3.template triangularView<StrictlyUpper>().toDenseMatrix(), m1);
m1.setZero(); m1.setZero();
m1.template triangularView<StrictlyLowerTriangular>() = 3 * m2; m1.template triangularView<StrictlyLower>() = 3 * m2;
VERIFY_IS_APPROX(m3.template triangularView<StrictlyLowerTriangular>().toDenseMatrix(), m1); VERIFY_IS_APPROX(m3.template triangularView<StrictlyLower>().toDenseMatrix(), m1);
m1.setRandom(); m1.setRandom();
m2 = m1.template triangularView<UpperTriangular>(); m2 = m1.template triangularView<Upper>();
VERIFY(m2.isUpperTriangular()); VERIFY(m2.isUpper());
VERIFY(!m2.isLowerTriangular()); VERIFY(!m2.isLower());
m2 = m1.template triangularView<StrictlyUpperTriangular>(); m2 = m1.template triangularView<StrictlyUpper>();
VERIFY(m2.isUpperTriangular()); VERIFY(m2.isUpper());
VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1))); VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1)));
m2 = m1.template triangularView<UnitUpperTriangular>(); m2 = m1.template triangularView<UnitUpper>();
VERIFY(m2.isUpperTriangular()); VERIFY(m2.isUpper());
m2.diagonal().array() -= Scalar(1); m2.diagonal().array() -= Scalar(1);
VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1))); VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1)));
m2 = m1.template triangularView<LowerTriangular>(); m2 = m1.template triangularView<Lower>();
VERIFY(m2.isLowerTriangular()); VERIFY(m2.isLower());
VERIFY(!m2.isUpperTriangular()); VERIFY(!m2.isUpper());
m2 = m1.template triangularView<StrictlyLowerTriangular>(); m2 = m1.template triangularView<StrictlyLower>();
VERIFY(m2.isLowerTriangular()); VERIFY(m2.isLower());
VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1))); VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1)));
m2 = m1.template triangularView<UnitLowerTriangular>(); m2 = m1.template triangularView<UnitLower>();
VERIFY(m2.isLowerTriangular()); VERIFY(m2.isLower());
m2.diagonal().array() -= Scalar(1); m2.diagonal().array() -= Scalar(1);
VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1))); VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1)));
// test swap // test swap
m1.setOnes(); m1.setOnes();
m2.setZero(); m2.setZero();
m2.template triangularView<UpperTriangular>().swap(m1); m2.template triangularView<Upper>().swap(m1);
m3.setZero(); m3.setZero();
m3.template triangularView<UpperTriangular>().setOnes(); m3.template triangularView<Upper>().setOnes();
VERIFY_IS_APPROX(m2,m3); VERIFY_IS_APPROX(m2,m3);
} }