diff --git a/Eigen/Core b/Eigen/Core
index 67d55fee3..1f71dfdb5 100644
--- a/Eigen/Core
+++ b/Eigen/Core
@@ -23,10 +23,6 @@
// License and a copy of the GNU General Public License along with
// Eigen. If not, see .
-#ifdef EIGEN2_SUPPORT
-#include "Eigen2Support"
-#endif
-
#ifndef EIGEN_CORE_H
#define EIGEN_CORE_H
@@ -225,5 +221,9 @@ struct Dense {};
#include "src/Core/util/EnableMSVCWarnings.h"
+#ifdef EIGEN2_SUPPORT
+#include "Eigen2Support"
+#endif
+
#endif // EIGEN_CORE_H
/* vim: set filetype=cpp et sw=2 ts=2 ai: */
diff --git a/Eigen/Eigen2Support b/Eigen/Eigen2Support
index b95b51dec..ccba4ff30 100644
--- a/Eigen/Eigen2Support
+++ b/Eigen/Eigen2Support
@@ -25,14 +25,13 @@
#ifndef EIGEN2SUPPORT_H
#define EIGEN2SUPPORT_H
-#include "src/Core/util/DisableMSVCWarnings.h"
+#include "Array"
-#ifndef EIGEN2_SUPPORT
-#define EIGEN2_SUPPORT
+#if (!defined(EIGEN2_SUPPORT)) || (!defined(EIGEN_CORE_H))
+#error Eigen2 support must be enabled by defining EIGEN2_SUPPORT before any other Eigen header
#endif
-#include "Core"
-#include "Array"
+#include "src/Core/util/DisableMSVCWarnings.h"
namespace Eigen {
@@ -40,8 +39,9 @@ namespace Eigen {
* This module provides a couple of deprecated functions improving the compatibility with Eigen2.
*
* \code
- * #include
+ * #define EIGEN2_SUPPORT
* \endcode
+ *
*/
#include "src/Eigen2Support/Flagged.h"
diff --git a/Eigen/src/Array/Reverse.h b/Eigen/src/Array/Reverse.h
index a8d8310a9..7d2c34816 100644
--- a/Eigen/src/Array/Reverse.h
+++ b/Eigen/src/Array/Reverse.h
@@ -59,9 +59,7 @@ struct ei_traits >
LinearAccess = ( (Direction==BothDirections) && (int(_MatrixTypeNested::Flags)&PacketAccessBit) )
? LinearAccessBit : 0,
- Flags = (int(_MatrixTypeNested::Flags) & (HereditaryBits | PacketAccessBit | LinearAccess))
- | (int(_MatrixTypeNested::Flags)&UpperTriangularBit ? LowerTriangularBit : 0)
- | (int(_MatrixTypeNested::Flags)&LowerTriangularBit ? UpperTriangularBit : 0),
+ Flags = (int(_MatrixTypeNested::Flags) & (HereditaryBits | PacketAccessBit | LinearAccess)),
CoeffReadCost = _MatrixTypeNested::CoeffReadCost
};
diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h
index 4fb6d3d2c..8b13d8e2e 100644
--- a/Eigen/src/Cholesky/LDLT.h
+++ b/Eigen/src/Cholesky/LDLT.h
@@ -80,7 +80,7 @@ template class LDLT
}
/** \returns the lower triangular matrix L */
- inline TriangularView matrixL(void) const
+ inline TriangularView matrixL(void) const
{
ei_assert(m_isInitialized && "LDLT is not initialized.");
return m_matrix;
@@ -301,13 +301,13 @@ bool LDLT::solveInPlace(MatrixBase &bAndX) const
// y = L^-1 z
//matrixL().solveInPlace(bAndX);
- m_matrix.template triangularView().solveInPlace(bAndX);
+ m_matrix.template triangularView().solveInPlace(bAndX);
// w = D^-1 y
bAndX = m_matrix.diagonal().asDiagonal().inverse() * bAndX;
// u = L^-T w
- m_matrix.adjoint().template triangularView().solveInPlace(bAndX);
+ m_matrix.adjoint().template triangularView().solveInPlace(bAndX);
// x = P^T u
for (int i = size-1; i >= 0; --i) bAndX.row(m_transpositions.coeff(i)).swap(bAndX.row(i));
diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h
index 02645b23f..8a149a316 100644
--- a/Eigen/src/Cholesky/LLT.h
+++ b/Eigen/src/Cholesky/LLT.h
@@ -148,7 +148,7 @@ template class LLT
// forward declaration (defined at the end of this file)
template struct ei_llt_inplace;
-template<> struct ei_llt_inplace
+template<> struct ei_llt_inplace
{
template
static bool unblocked(MatrixType& mat)
@@ -198,47 +198,47 @@ template<> struct ei_llt_inplace
Block A22(m,k+bs,k+bs,rs,rs);
if(!unblocked(A11)) return false;
- if(rs>0) A11.adjoint().template triangularView().template solveInPlace(A21);
- if(rs>0) A22.template selfadjointView().rankUpdate(A21,-1); // bottleneck
+ if(rs>0) A11.adjoint().template triangularView().template solveInPlace(A21);
+ if(rs>0) A22.template selfadjointView().rankUpdate(A21,-1); // bottleneck
}
return true;
}
};
-template<> struct ei_llt_inplace
+template<> struct ei_llt_inplace
{
template
static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat)
{
Transpose matt(mat);
- return ei_llt_inplace::unblocked(matt);
+ return ei_llt_inplace::unblocked(matt);
}
template
static EIGEN_STRONG_INLINE bool blocked(MatrixType& mat)
{
Transpose matt(mat);
- return ei_llt_inplace::blocked(matt);
+ return ei_llt_inplace::blocked(matt);
}
};
-template struct LLT_Traits
+template struct LLT_Traits
{
- typedef TriangularView MatrixL;
- typedef TriangularView MatrixU;
+ typedef TriangularView MatrixL;
+ typedef TriangularView MatrixU;
inline static MatrixL getL(const MatrixType& m) { return m; }
inline static MatrixU getU(const MatrixType& m) { return m.adjoint(); }
static bool inplace_decomposition(MatrixType& m)
- { return ei_llt_inplace::blocked(m); }
+ { return ei_llt_inplace::blocked(m); }
};
-template struct LLT_Traits
+template struct LLT_Traits
{
- typedef TriangularView MatrixL;
- typedef TriangularView MatrixU;
+ typedef TriangularView MatrixL;
+ typedef TriangularView MatrixU;
inline static MatrixL getL(const MatrixType& m) { return m.adjoint(); }
inline static MatrixU getU(const MatrixType& m) { return m; }
static bool inplace_decomposition(MatrixType& m)
- { return ei_llt_inplace::blocked(m); }
+ { return ei_llt_inplace::blocked(m); }
};
/** Computes / recomputes the Cholesky decomposition A = LL^* = U^*U of \a matrix
diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h
index 4c30f30ad..7569dab58 100644
--- a/Eigen/src/Core/MatrixBase.h
+++ b/Eigen/src/Core/MatrixBase.h
@@ -405,13 +405,6 @@ template class MatrixBase
template const VectorBlock start() const;
template VectorBlock end();
template const VectorBlock end() const;
-
- template
- typename ei_plain_matrix_type_column_major::type
- solveTriangular(const MatrixBase& other) const;
-
- template
- void solveTriangularInPlace(const MatrixBase& other) const;
#endif
protected:
diff --git a/Eigen/src/Core/SelfAdjointView.h b/Eigen/src/Core/SelfAdjointView.h
index 9ab8cb2c1..6d01ee495 100644
--- a/Eigen/src/Core/SelfAdjointView.h
+++ b/Eigen/src/Core/SelfAdjointView.h
@@ -31,7 +31,7 @@
* \brief Expression of a selfadjoint matrix from a triangular part of a dense matrix
*
* \param MatrixType the type of the dense matrix storing the coefficients
- * \param TriangularPart can be either \c LowerTriangular or \c UpperTriangular
+ * \param TriangularPart can be either \c Lower or \c Upper
*
* This class is an expression of a sefladjoint matrix from a triangular part of a matrix
* with given dense storage of the coefficients. It is the return type of MatrixBase::selfadjointView()
@@ -46,7 +46,7 @@ struct ei_traits > : ei_traits::type _MatrixTypeNested;
typedef MatrixType ExpressionType;
enum {
- Mode = TriangularPart | SelfAdjointBit,
+ Mode = TriangularPart | SelfAdjoint,
Flags = _MatrixTypeNested::Flags & (HereditaryBits)
& (~(PacketAccessBit | DirectAccessBit | LinearAccessBit)), // FIXME these flags should be preserved
CoeffReadCost = _MatrixTypeNested::CoeffReadCost
diff --git a/Eigen/src/Core/SolveTriangular.h b/Eigen/src/Core/SolveTriangular.h
index 9dc019d17..73cf77387 100644
--- a/Eigen/src/Core/SolveTriangular.h
+++ b/Eigen/src/Core/SolveTriangular.h
@@ -57,7 +57,7 @@ struct ei_triangular_solver_selector LhsProductTraits;
typedef typename LhsProductTraits::ExtractType ActualLhsType;
enum {
- IsLowerTriangular = ((Mode&LowerTriangularBit)==LowerTriangularBit)
+ IsLower = ((Mode&Lower)==Lower)
};
static void run(const Lhs& lhs, Rhs& other)
{
@@ -65,20 +65,20 @@ struct ei_triangular_solver_selector0;
- IsLowerTriangular ? pi+=PanelWidth : pi-=PanelWidth)
+ for(int pi=IsLower ? 0 : size;
+ IsLower ? pi0;
+ IsLower ? pi+=PanelWidth : pi-=PanelWidth)
{
- int actualPanelWidth = std::min(IsLowerTriangular ? size - pi : pi, PanelWidth);
+ int actualPanelWidth = std::min(IsLower ? size - pi : pi, PanelWidth);
- int r = IsLowerTriangular ? pi : size - pi; // remaining size
+ int r = IsLower ? pi : size - pi; // remaining size
if (r > 0)
{
// let's directly call the low level product function because:
// 1 - it is faster to compile
// 2 - it is slighlty faster at runtime
- int startRow = IsLowerTriangular ? pi : pi-actualPanelWidth;
- int startCol = IsLowerTriangular ? 0 : pi;
+ int startRow = IsLower ? pi : pi-actualPanelWidth;
+ int startCol = IsLower ? 0 : pi;
VectorBlock target(other,startRow,actualPanelWidth);
ei_cache_friendly_product_rowmajor_times_vector(
@@ -89,12 +89,12 @@ struct ei_triangular_solver_selector0)
other.coeffRef(i) -= (lhs.row(i).segment(s,k).transpose().cwiseProduct(other.segment(s,k))).sum();
- if(!(Mode & UnitDiagBit))
+ if(!(Mode & UnitDiag))
other.coeffRef(i) /= lhs.coeff(i,i);
}
}
@@ -111,7 +111,7 @@ struct ei_triangular_solver_selector::size,
- IsLowerTriangular = ((Mode&LowerTriangularBit)==LowerTriangularBit)
+ IsLower = ((Mode&Lower)==Lower)
};
static void run(const Lhs& lhs, Rhs& other)
@@ -120,26 +120,26 @@ struct ei_triangular_solver_selector0;
- IsLowerTriangular ? pi+=PanelWidth : pi-=PanelWidth)
+ for(int pi=IsLower ? 0 : size;
+ IsLower ? pi0;
+ IsLower ? pi+=PanelWidth : pi-=PanelWidth)
{
- int actualPanelWidth = std::min(IsLowerTriangular ? size - pi : pi, PanelWidth);
- int startBlock = IsLowerTriangular ? pi : pi-actualPanelWidth;
- int endBlock = IsLowerTriangular ? pi + actualPanelWidth : 0;
+ int actualPanelWidth = std::min(IsLower ? size - pi : pi, PanelWidth);
+ int startBlock = IsLower ? pi : pi-actualPanelWidth;
+ int endBlock = IsLower ? pi + actualPanelWidth : 0;
for(int k=0; k0)
other.segment(s,r) -= other.coeffRef(i) * Block(lhs, s, i, r, 1);
}
- int r = IsLowerTriangular ? size - endBlock : startBlock; // remaining size
+ int r = IsLower ? size - endBlock : startBlock; // remaining size
if (r > 0)
{
// let's directly call the low level product function because:
@@ -198,16 +198,16 @@ struct ei_triangular_solver_unroller;
template
struct ei_triangular_solver_unroller {
enum {
- IsLowerTriangular = ((Mode&LowerTriangularBit)==LowerTriangularBit),
- I = IsLowerTriangular ? Index : Size - Index - 1,
- S = IsLowerTriangular ? 0 : I+1
+ IsLower = ((Mode&Lower)==Lower),
+ I = IsLower ? Index : Size - Index - 1,
+ S = IsLower ? 0 : I+1
};
static void run(const Lhs& lhs, Rhs& rhs)
{
if (Index>0)
rhs.coeffRef(I) -= ((lhs.row(I).template segment(S).transpose()).cwiseProduct(rhs.template segment(S))).sum();
- if(!(Mode & UnitDiagBit))
+ if(!(Mode & UnitDiag))
rhs.coeffRef(I) /= lhs.coeff(I,I);
ei_triangular_solver_unroller::run(lhs,rhs);
@@ -245,8 +245,8 @@ void TriangularView::solveInPlace(const MatrixBase::Flags & RowMajorBit && OtherDerived::IsVectorAtCompileTime };
typedef typename ei_meta_if > : ei_traits
ColsAtCompileTime = MatrixType::RowsAtCompileTime,
MaxRowsAtCompileTime = MatrixType::MaxColsAtCompileTime,
MaxColsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
- Flags = ((int(_MatrixTypeNested::Flags) ^ RowMajorBit)
- & ~(LowerTriangularBit | UpperTriangularBit))
- | (int(_MatrixTypeNested::Flags)&UpperTriangularBit ? LowerTriangularBit : 0)
- | (int(_MatrixTypeNested::Flags)&LowerTriangularBit ? UpperTriangularBit : 0),
+ Flags = (int(_MatrixTypeNested::Flags) ^ RowMajorBit),
CoeffReadCost = _MatrixTypeNested::CoeffReadCost
};
};
@@ -229,7 +226,7 @@ struct ei_inplace_transpose_selector;
template
struct ei_inplace_transpose_selector { // square matrix
static void run(MatrixType& m) {
- m.template triangularView().swap(m.transpose());
+ m.template triangularView().swap(m.transpose());
}
};
@@ -237,7 +234,7 @@ template
struct ei_inplace_transpose_selector { // non square matrix
static void run(MatrixType& m) {
if (m.rows()==m.cols())
- m.template triangularView().swap(m.transpose());
+ m.template triangularView().swap(m.transpose());
else
m = m.transpose().eval();
}
diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h
index 62d800fef..8411b0546 100644
--- a/Eigen/src/Core/TriangularMatrix.h
+++ b/Eigen/src/Core/TriangularMatrix.h
@@ -46,7 +46,7 @@ template class TriangularBase : public AnyMatrixBase
};
typedef typename ei_traits::Scalar Scalar;
- inline TriangularBase() { ei_assert(ei_are_flags_consistent::ret); }
+ inline TriangularBase() { ei_assert(!((Mode&UnitDiag) && (Mode&ZeroDiag))); }
inline int rows() const { return derived().rows(); }
inline int cols() const { return derived().cols(); }
@@ -89,10 +89,10 @@ template class TriangularBase : public AnyMatrixBase
void check_coordinates(int row, int col)
{
ei_assert(col>=0 && col=0 && row=row)
- || (Mode==LowerTriangular && col<=row)
- || ((Mode==StrictlyUpperTriangular || Mode==UnitUpperTriangular) && col>row)
- || ((Mode==StrictlyLowerTriangular || Mode==UnitLowerTriangular) && col=row)
+ || (Mode==Lower && col<=row)
+ || ((Mode==StrictlyUpper || Mode==UnitUpper) && col>row)
+ || ((Mode==StrictlyLower || Mode==UnitLower) && col class TriangularBase : public AnyMatrixBase
* \brief Base class for triangular part in a matrix
*
* \param MatrixType the type of the object in which we are taking the triangular part
- * \param Mode the kind of triangular matrix expression to construct. Can be UpperTriangular,
- * LowerTriangular, UpperSelfadjoint, or LowerSelfadjoint. This is in fact a bit field;
- * it must have either UpperBit or LowerBit, and additionnaly it may have either
- * TraingularBit or SelfadjointBit.
+ * \param Mode the kind of triangular matrix expression to construct. Can be Upper,
+ * Lower, UpperSelfadjoint, or LowerSelfadjoint. This is in fact a bit field;
+ * it must have either Upper or Lower, and additionnaly it may have either
+ * UnitDiag or Selfadjoint.
*
* This class represents a triangular part of a matrix, not necessarily square. Strictly speaking, for rectangular
* matrices one should speak ok "trapezoid" parts. This class is the return type
@@ -154,9 +154,10 @@ template class TriangularView
enum {
Mode = _Mode,
- TransposeMode = (Mode & UpperTriangularBit ? LowerTriangularBit : 0)
- | (Mode & LowerTriangularBit ? UpperTriangularBit : 0)
- | (Mode & (ZeroDiagBit | UnitDiagBit))
+ TransposeMode = (Mode & Upper ? Lower : 0)
+ | (Mode & Lower ? Upper : 0)
+ | (Mode & (UnitDiag))
+ | (Mode & (ZeroDiag))
};
inline TriangularView(const MatrixType& matrix) : m_matrix(matrix)
@@ -283,12 +284,12 @@ template class TriangularView
const SelfAdjointView<_MatrixTypeNested,Mode> selfadjointView() const
{
- EIGEN_STATIC_ASSERT((Mode&UnitDiagBit)==0,PROGRAMMING_ERROR);
+ EIGEN_STATIC_ASSERT((Mode&UnitDiag)==0,PROGRAMMING_ERROR);
return SelfAdjointView<_MatrixTypeNested,Mode>(m_matrix);
}
SelfAdjointView<_MatrixTypeNested,Mode> selfadjointView()
{
- EIGEN_STATIC_ASSERT((Mode&UnitDiagBit)==0,PROGRAMMING_ERROR);
+ EIGEN_STATIC_ASSERT((Mode&UnitDiag)==0,PROGRAMMING_ERROR);
return SelfAdjointView<_MatrixTypeNested,Mode>(m_matrix);
}
@@ -304,6 +305,16 @@ template class TriangularView
TriangularView,Mode>(const_cast(m_matrix)).lazyAssign(other.derived());
}
+ Scalar determinant() const
+ {
+ if (Mode & UnitDiag)
+ return 1;
+ else if (Mode & ZeroDiag)
+ return 0;
+ else
+ return m_matrix.diagonal().prod();
+ }
+
protected:
const MatrixTypeNested m_matrix;
@@ -325,19 +336,19 @@ struct ei_triangular_assignment_selector
{
ei_triangular_assignment_selector::run(dst, src);
- ei_assert( Mode == UpperTriangular || Mode == LowerTriangular
- || Mode == StrictlyUpperTriangular || Mode == StrictlyLowerTriangular
- || Mode == UnitUpperTriangular || Mode == UnitLowerTriangular);
- if((Mode == UpperTriangular && row <= col)
- || (Mode == LowerTriangular && row >= col)
- || (Mode == StrictlyUpperTriangular && row < col)
- || (Mode == StrictlyLowerTriangular && row > col)
- || (Mode == UnitUpperTriangular && row < col)
- || (Mode == UnitLowerTriangular && row > col))
+ ei_assert( Mode == Upper || Mode == Lower
+ || Mode == StrictlyUpper || Mode == StrictlyLower
+ || Mode == UnitUpper || Mode == UnitLower);
+ if((Mode == Upper && row <= col)
+ || (Mode == Lower && row >= col)
+ || (Mode == StrictlyUpper && row < col)
+ || (Mode == StrictlyLower && row > col)
+ || (Mode == UnitUpper && row < col)
+ || (Mode == UnitLower && row > col))
dst.copyCoeff(row, col, src);
else if(ClearOpposite)
{
- if (Mode&UnitDiagBit && row==col)
+ if (Mode&UnitDiag && row==col)
dst.coeffRef(row, col) = 1;
else
dst.coeffRef(row, col) = 0;
@@ -353,7 +364,7 @@ struct ei_triangular_assignment_selector
-struct ei_triangular_assignment_selector
+struct ei_triangular_assignment_selector
{
inline static void run(Derived1 &dst, const Derived2 &src)
{
@@ -370,7 +381,7 @@ struct ei_triangular_assignment_selector
-struct ei_triangular_assignment_selector
+struct ei_triangular_assignment_selector
{
inline static void run(Derived1 &dst, const Derived2 &src)
{
@@ -387,7 +398,7 @@ struct ei_triangular_assignment_selector
-struct ei_triangular_assignment_selector
+struct ei_triangular_assignment_selector
{
inline static void run(Derived1 &dst, const Derived2 &src)
{
@@ -404,7 +415,7 @@ struct ei_triangular_assignment_selector
-struct ei_triangular_assignment_selector
+struct ei_triangular_assignment_selector
{
inline static void run(Derived1 &dst, const Derived2 &src)
{
@@ -421,7 +432,7 @@ struct ei_triangular_assignment_selector
-struct ei_triangular_assignment_selector
+struct ei_triangular_assignment_selector
{
inline static void run(Derived1 &dst, const Derived2 &src)
{
@@ -440,7 +451,7 @@ struct ei_triangular_assignment_selector
-struct ei_triangular_assignment_selector
+struct ei_triangular_assignment_selector
{
inline static void run(Derived1 &dst, const Derived2 &src)
{
@@ -590,8 +601,8 @@ EIGEN_DEPRECATED TriangularView MatrixBase::part()
/** \nonstableyet
* \returns an expression of a triangular view extracted from the current matrix
*
- * The parameter \a Mode can have the following values: \c UpperTriangular, \c StrictlyUpperTriangular, \c UnitUpperTriangular,
- * \c LowerTriangular, \c StrictlyLowerTriangular, \c UnitLowerTriangular.
+ * The parameter \a Mode can have the following values: \c Upper, \c StrictlyUpper, \c UnitUpper,
+ * \c Lower, \c StrictlyLower, \c UnitLower.
*
* Example: \include MatrixBase_extract.cpp
* Output: \verbinclude MatrixBase_extract.out
@@ -616,22 +627,22 @@ const TriangularView MatrixBase::triangularView() const
/** \returns true if *this is approximately equal to an upper triangular matrix,
* within the precision given by \a prec.
*
- * \sa isLowerTriangular(), extract(), part(), marked()
+ * \sa isLower(), extract(), part(), marked()
*/
template
bool MatrixBase::isUpperTriangular(RealScalar prec) const
{
- RealScalar maxAbsOnUpperTriangularPart = static_cast(-1);
+ RealScalar maxAbsOnUpperPart = static_cast(-1);
for(int j = 0; j < cols(); ++j)
{
int maxi = std::min(j, rows()-1);
for(int i = 0; i <= maxi; ++i)
{
RealScalar absValue = ei_abs(coeff(i,j));
- if(absValue > maxAbsOnUpperTriangularPart) maxAbsOnUpperTriangularPart = absValue;
+ if(absValue > maxAbsOnUpperPart) maxAbsOnUpperPart = absValue;
}
}
- RealScalar threshold = maxAbsOnUpperTriangularPart * prec;
+ RealScalar threshold = maxAbsOnUpperPart * prec;
for(int j = 0; j < cols(); ++j)
for(int i = j+1; i < rows(); ++i)
if(ei_abs(coeff(i, j)) > threshold) return false;
@@ -641,19 +652,19 @@ bool MatrixBase::isUpperTriangular(RealScalar prec) const
/** \returns true if *this is approximately equal to a lower triangular matrix,
* within the precision given by \a prec.
*
- * \sa isUpperTriangular(), extract(), part(), marked()
+ * \sa isUpper(), extract(), part(), marked()
*/
template
bool MatrixBase::isLowerTriangular(RealScalar prec) const
{
- RealScalar maxAbsOnLowerTriangularPart = static_cast(-1);
+ RealScalar maxAbsOnLowerPart = static_cast(-1);
for(int j = 0; j < cols(); ++j)
for(int i = j; i < rows(); ++i)
{
RealScalar absValue = ei_abs(coeff(i,j));
- if(absValue > maxAbsOnLowerTriangularPart) maxAbsOnLowerTriangularPart = absValue;
+ if(absValue > maxAbsOnLowerPart) maxAbsOnLowerPart = absValue;
}
- RealScalar threshold = maxAbsOnLowerTriangularPart * prec;
+ RealScalar threshold = maxAbsOnLowerPart * prec;
for(int j = 1; j < cols(); ++j)
{
int maxi = std::min(j, rows()-1);
diff --git a/Eigen/src/Core/products/SelfadjointMatrixMatrix.h b/Eigen/src/Core/products/SelfadjointMatrixMatrix.h
index 35efa752e..ac072e189 100644
--- a/Eigen/src/Core/products/SelfadjointMatrixMatrix.h
+++ b/Eigen/src/Core/products/SelfadjointMatrixMatrix.h
@@ -368,10 +368,10 @@ struct SelfadjointProductMatrix
SelfadjointProductMatrix(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) {}
enum {
- LhsUpLo = LhsMode&(UpperTriangularBit|LowerTriangularBit),
- LhsIsSelfAdjoint = (LhsMode&SelfAdjointBit)==SelfAdjointBit,
- RhsUpLo = RhsMode&(UpperTriangularBit|LowerTriangularBit),
- RhsIsSelfAdjoint = (RhsMode&SelfAdjointBit)==SelfAdjointBit
+ LhsUpLo = LhsMode&(Upper|Lower),
+ LhsIsSelfAdjoint = (LhsMode&SelfAdjoint)==SelfAdjoint,
+ RhsUpLo = RhsMode&(Upper|Lower),
+ RhsIsSelfAdjoint = (RhsMode&SelfAdjoint)==SelfAdjoint
};
template void scaleAndAddTo(Dest& dst, Scalar alpha) const
@@ -385,12 +385,12 @@ struct SelfadjointProductMatrix
* RhsBlasTraits::extractScalarFactor(m_rhs);
ei_product_selfadjoint_matrix::Flags &RowMajorBit) ? RowMajor : ColMajor, LhsIsSelfAdjoint,
- NumTraits::IsComplex && EIGEN_LOGICAL_XOR(LhsUpLo==UpperTriangular,bool(LhsBlasTraits::NeedToConjugate)),
- EIGEN_LOGICAL_XOR(RhsUpLo==UpperTriangular,
+ NumTraits::IsComplex && EIGEN_LOGICAL_XOR(LhsUpLo==Upper,bool(LhsBlasTraits::NeedToConjugate)),
+ EIGEN_LOGICAL_XOR(RhsUpLo==Upper,
ei_traits::Flags &RowMajorBit) ? RowMajor : ColMajor, RhsIsSelfAdjoint,
- NumTraits::IsComplex && EIGEN_LOGICAL_XOR(RhsUpLo==UpperTriangular,bool(RhsBlasTraits::NeedToConjugate)),
+ NumTraits::IsComplex && EIGEN_LOGICAL_XOR(RhsUpLo==Upper,bool(RhsBlasTraits::NeedToConjugate)),
ei_traits::Flags&RowMajorBit ? RowMajor : ColMajor>
::run(
lhs.rows(), rhs.cols(), // sizes
diff --git a/Eigen/src/Core/products/SelfadjointMatrixVector.h b/Eigen/src/Core/products/SelfadjointMatrixVector.h
index 32b7f220e..1c48208b3 100644
--- a/Eigen/src/Core/products/SelfadjointMatrixVector.h
+++ b/Eigen/src/Core/products/SelfadjointMatrixVector.h
@@ -42,7 +42,7 @@ static EIGEN_DONT_INLINE void ei_product_selfadjoint_vector(
enum {
IsRowMajor = StorageOrder==RowMajor ? 1 : 0,
- IsLower = UpLo == LowerTriangularBit ? 1 : 0,
+ IsLower = UpLo == Lower ? 1 : 0,
FirstTriangular = IsRowMajor == IsLower
};
@@ -170,7 +170,7 @@ struct SelfadjointProductMatrix
EIGEN_PRODUCT_PUBLIC_INTERFACE(SelfadjointProductMatrix)
enum {
- LhsUpLo = LhsMode&(UpperTriangularBit|LowerTriangularBit)
+ LhsUpLo = LhsMode&(Upper|Lower)
};
SelfadjointProductMatrix(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) {}
diff --git a/Eigen/src/Core/products/SelfadjointProduct.h b/Eigen/src/Core/products/SelfadjointProduct.h
index f9cdda637..8967f62be 100644
--- a/Eigen/src/Core/products/SelfadjointProduct.h
+++ b/Eigen/src/Core/products/SelfadjointProduct.h
@@ -47,7 +47,7 @@ struct ei_selfadjoint_product
{
static EIGEN_STRONG_INLINE void run(int size, int depth, const Scalar* mat, int matStride, Scalar* res, int resStride, Scalar alpha)
{
- ei_selfadjoint_product
+ ei_selfadjoint_product
::run(size, depth, mat, matStride, res, resStride, alpha);
}
};
@@ -100,13 +100,13 @@ struct ei_selfadjoint_product
// 1 - before the diagonal => processed with gebp or skipped
// 2 - the actual_mc x actual_mc symmetric block => processed with a special kernel
// 3 - after the diagonal => processed with gebp or skipped
- if (UpLo==LowerTriangular)
+ if (UpLo==Lower)
gebp_kernel(res+i2, resStride, blockA, blockB, actual_mc, actual_kc, std::min(size,i2));
ei_sybb_kernel()
(res+resStride*i2 + i2, resStride, blockA, blockB + actual_kc*Blocking::PacketSize*i2, actual_mc, actual_kc);
- if (UpLo==UpperTriangular)
+ if (UpLo==Upper)
{
int j2 = i2+actual_mc;
gebp_kernel(res+resStride*j2+i2, resStride, blockA, blockB+actual_kc*Blocking::PacketSize*j2, actual_mc, actual_kc, std::max(0,size-j2));
@@ -173,7 +173,7 @@ struct ei_sybb_kernel
int actualBlockSize = std::min(BlockSize,size - j);
const Scalar* actual_b = blockB+j*depth*PacketSize;
- if(UpLo==UpperTriangular)
+ if(UpLo==Upper)
gebp_kernel(res+j*resStride, resStride, blockA, actual_b, j, depth, actualBlockSize);
// selfadjoint micro block
@@ -186,13 +186,13 @@ struct ei_sybb_kernel
for(int j1=0; j1
struct ei_selfadjoint_rank2_update_selector;
template
-struct ei_selfadjoint_rank2_update_selector
+struct ei_selfadjoint_rank2_update_selector
{
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
};
template
-struct ei_selfadjoint_rank2_update_selector
+struct ei_selfadjoint_rank2_update_selector
{
static void run(Scalar* mat, int stride, const UType& u, const VType& v, Scalar alpha)
{
@@ -87,7 +87,7 @@ SelfAdjointView& SelfAdjointView
ei_selfadjoint_rank2_update_selector::ret>::type,
typename ei_cleantype::ret>::type,
- (IsRowMajor ? (UpLo==UpperTriangular ? LowerTriangular : UpperTriangular) : UpLo)>
+ (IsRowMajor ? (UpLo==Upper ? Lower : Upper) : UpLo)>
::run(const_cast(_expression().data()),_expression().stride(),actualU,actualV,actualAlpha);
return *this;
diff --git a/Eigen/src/Core/products/TriangularMatrixMatrix.h b/Eigen/src/Core/products/TriangularMatrixMatrix.h
index 701ccb644..37617a915 100644
--- a/Eigen/src/Core/products/TriangularMatrixMatrix.h
+++ b/Eigen/src/Core/products/TriangularMatrixMatrix.h
@@ -75,7 +75,7 @@ struct ei_product_triangular_matrix_matrix Blocking;
enum {
SmallPanelWidth = EIGEN_ENUM_MAX(Blocking::mr,Blocking::nr),
- IsLowerTriangular = (Mode&LowerTriangular) == LowerTriangular
+ IsLower = (Mode&Lower) == Lower
};
int kc = std::min(Blocking::Max_kc/4,size); // cache block size along the K direction
@@ -130,12 +130,12 @@ struct ei_product_triangular_matrix_matrix pack_lhs;
ei_gemm_pack_rhs pack_rhs;
- for(int k2=IsLowerTriangular ? size : 0;
- IsLowerTriangular ? k2>0 : k20 : k2(actual_kc-k1, SmallPanelWidth);
- int lengthTarget = IsLowerTriangular ? actual_kc-k1-actualPanelWidth : k1;
+ int lengthTarget = IsLower ? actual_kc-k1-actualPanelWidth : k1;
int startBlock = actual_k2+k1;
int blockBOffset = k1;
@@ -158,9 +158,9 @@ struct ei_product_triangular_matrix_matrix0)
{
- int startTarget = IsLowerTriangular ? actual_k2+k1+actualPanelWidth : actual_k2;
+ int startTarget = IsLower ? actual_k2+k1+actualPanelWidth : actual_k2;
pack_lhs(blockA, &lhs(startTarget,startBlock), lhsStride, actualPanelWidth, lengthTarget);
@@ -182,8 +182,8 @@ struct ei_product_triangular_matrix_matrix GEPP
{
- int start = IsLowerTriangular ? k2 : 0;
- int end = IsLowerTriangular ? size : actual_k2;
+ int start = IsLower ? k2 : 0;
+ int end = IsLower ? size : actual_k2;
for(int i2=start; i2 Blocking;
enum {
SmallPanelWidth = EIGEN_ENUM_MAX(Blocking::mr,Blocking::nr),
- IsLowerTriangular = (Mode&LowerTriangular) == LowerTriangular
+ IsLower = (Mode&Lower) == Lower
};
int kc = std::min(Blocking::Max_kc/4,size); // cache block size along the K direction
@@ -245,16 +245,16 @@ struct ei_product_triangular_matrix_matrix pack_rhs;
ei_gemm_pack_rhs pack_rhs_panel;
- for(int k2=IsLowerTriangular ? 0 : size;
- IsLowerTriangular ? k20;
- IsLowerTriangular ? k2+=kc : k2-=kc)
+ for(int k2=IsLower ? 0 : size;
+ IsLower ? k20;
+ IsLower ? k2+=kc : k2-=kc)
{
- const int actual_kc = std::min(IsLowerTriangular ? size-k2 : k2, kc);
- int actual_k2 = IsLowerTriangular ? k2 : k2-actual_kc;
- int rs = IsLowerTriangular ? actual_k2 : size - k2;
+ const int actual_kc = std::min(IsLower ? size-k2 : k2, kc);
+ int actual_k2 = IsLower ? k2 : k2-actual_kc;
+ int rs = IsLower ? actual_k2 : size - k2;
Scalar* geb = blockB+actual_kc*actual_kc*Blocking::PacketSize;
- pack_rhs(geb, &rhs(actual_k2,IsLowerTriangular ? 0 : k2), rhsStride, alpha, actual_kc, rs);
+ pack_rhs(geb, &rhs(actual_k2,IsLower ? 0 : k2), rhsStride, alpha, actual_kc, rs);
// pack the triangular part of the rhs padding the unrolled blocks with zeros
{
@@ -262,8 +262,8 @@ struct ei_product_triangular_matrix_matrix(actual_kc-j2, SmallPanelWidth);
int actual_j2 = actual_k2 + j2;
- int panelOffset = IsLowerTriangular ? j2+actualPanelWidth : 0;
- int panelLength = IsLowerTriangular ? actual_kc-j2-actualPanelWidth : j2;
+ int panelOffset = IsLower ? j2+actualPanelWidth : 0;
+ int panelLength = IsLower ? actual_kc-j2-actualPanelWidth : j2;
// general part
pack_rhs_panel(blockB+j2*actual_kc*Blocking::PacketSize,
&rhs(actual_k2+panelOffset, actual_j2), rhsStride, alpha,
@@ -273,9 +273,9 @@ struct ei_product_triangular_matrix_matrix(actual_kc-j2, SmallPanelWidth);
- int panelLength = IsLowerTriangular ? actual_kc-j2 : j2+actualPanelWidth;
- int blockOffset = IsLowerTriangular ? j2 : 0;
+ int panelLength = IsLower ? actual_kc-j2 : j2+actualPanelWidth;
+ int blockOffset = IsLower ? j2 : 0;
gebp_kernel(res+i2+(actual_k2+j2)*resStride, resStride,
blockA, blockB+j2*actual_kc*Blocking::PacketSize,
@@ -306,7 +306,7 @@ struct ei_product_triangular_matrix_matrix::Scalar alpha)
{
@@ -50,17 +50,17 @@ struct ei_product_triangular_vector_selector0)
res.segment(s,r) += (alpha * cjRhs.coeff(i)) * cjLhs.col(i).segment(s,r);
if (HasUnitDiag)
res.coeffRef(i) += alpha * cjRhs.coeff(i);
}
- int r = IsLowerTriangular ? size - pi - actualPanelWidth : pi;
+ int r = IsLower ? size - pi - actualPanelWidth : pi;
if (r>0)
{
- int s = IsLowerTriangular ? pi+actualPanelWidth : 0;
+ int s = IsLower ? pi+actualPanelWidth : 0;
ei_cache_friendly_product_colmajor_times_vector(
r,
&(lhs.const_cast_derived().coeffRef(s,pi)), lhs.stride(),
@@ -77,8 +77,8 @@ struct ei_product_triangular_vector_selector::Scalar alpha)
{
@@ -92,17 +92,17 @@ struct ei_product_triangular_vector_selector0)
res.coeffRef(i) += alpha * (cjLhs.row(i).segment(s,r).cwiseProduct(cjRhs.segment(s,r).transpose())).sum();
if (HasUnitDiag)
res.coeffRef(i) += alpha * cjRhs.coeff(i);
}
- int r = IsLowerTriangular ? pi : size - pi - actualPanelWidth;
+ int r = IsLower ? pi : size - pi - actualPanelWidth;
if (r>0)
{
- int s = IsLowerTriangular ? 0 : pi + actualPanelWidth;
+ int s = IsLower ? 0 : pi + actualPanelWidth;
Block target(res,pi,0,actualPanelWidth,1);
ei_cache_friendly_product_rowmajor_times_vector(
&(lhs.const_cast_derived().coeffRef(pi,s)), lhs.stride(),
diff --git a/Eigen/src/Core/products/TriangularSolverMatrix.h b/Eigen/src/Core/products/TriangularSolverMatrix.h
index e49fac956..23a645d7c 100644
--- a/Eigen/src/Core/products/TriangularSolverMatrix.h
+++ b/Eigen/src/Core/products/TriangularSolverMatrix.h
@@ -36,7 +36,7 @@ struct ei_triangular_solve_matrix::IsComplex && Conjugate,
TriStorageOrder==RowMajor ? ColMajor : RowMajor, ColMajor>
::run(size, cols, tri, triStride, _other, otherStride);
@@ -60,7 +60,7 @@ struct ei_triangular_solve_matrix Blocking;
enum {
SmallPanelWidth = EIGEN_ENUM_MAX(Blocking::mr,Blocking::nr),
- IsLowerTriangular = (Mode&LowerTriangular) == LowerTriangular
+ IsLower = (Mode&Lower) == Lower
};
int kc = std::min(Blocking::Max_kc/4,size); // cache block size along the K direction
@@ -73,11 +73,11 @@ struct ei_triangular_solve_matrix > gebp_kernel;
ei_gemm_pack_lhs pack_lhs;
- for(int k2=IsLowerTriangular ? 0 : size;
- IsLowerTriangular ? k20;
- IsLowerTriangular ? k2+=kc : k2-=kc)
+ for(int k2=IsLower ? 0 : size;
+ IsLower ? k20;
+ IsLower ? k2+=kc : k2-=kc)
{
- const int actual_kc = std::min(IsLowerTriangular ? size-k2 : k2, kc);
+ const int actual_kc = std::min(IsLower ? size-k2 : k2, kc);
// We have selected and packed a big horizontal panel R1 of rhs. Let B be the packed copy of this panel,
// and R2 the remaining part of rhs. The corresponding vertical panel of lhs is split into
@@ -101,11 +101,11 @@ struct ei_triangular_solve_matrix()
@@ -141,7 +141,7 @@ struct ei_triangular_solve_matrix0)
{
- int startTarget = IsLowerTriangular ? k2+k1+actualPanelWidth : k2-actual_kc;
+ int startTarget = IsLower ? k2+k1+actualPanelWidth : k2-actual_kc;
pack_lhs(blockA, &tri(startTarget,startBlock), triStride, actualPanelWidth, lengthTarget);
@@ -153,14 +153,14 @@ struct ei_triangular_solve_matrix GEPP
{
- int start = IsLowerTriangular ? k2+kc : 0;
- int end = IsLowerTriangular ? size : k2-kc;
+ int start = IsLower ? k2+kc : 0;
+ int end = IsLower ? size : k2-kc;
for(int i2=start; i20)
{
- pack_lhs(blockA, &tri(i2, IsLowerTriangular ? k2 : k2-kc), triStride, actual_kc, actual_mc);
+ pack_lhs(blockA, &tri(i2, IsLower ? k2 : k2-kc), triStride, actual_kc, actual_mc);
gebp_kernel(_other+i2, otherStride, blockA, blockB, actual_mc, actual_kc, cols);
}
@@ -191,7 +191,7 @@ struct ei_triangular_solve_matrix(Blocking::Max_kc/4,size); // cache block size along the K direction
@@ -206,15 +206,15 @@ struct ei_triangular_solve_matrix pack_rhs_panel;
ei_gemm_pack_lhs pack_lhs_panel;
- for(int k2=IsLowerTriangular ? size : 0;
- IsLowerTriangular ? k2>0 : k20 : k20) pack_rhs(geb, &rhs(actual_k2,startPanel), triStride, -1, actual_kc, rs);
@@ -226,8 +226,8 @@ struct ei_triangular_solve_matrix(actual_kc-j2, SmallPanelWidth);
int actual_j2 = actual_k2 + j2;
- int panelOffset = IsLowerTriangular ? j2+actualPanelWidth : 0;
- int panelLength = IsLowerTriangular ? actual_kc-j2-actualPanelWidth : j2;
+ int panelOffset = IsLower ? j2+actualPanelWidth : 0;
+ int panelLength = IsLower ? actual_kc-j2-actualPanelWidth : j2;
if (panelLength>0)
pack_rhs_panel(blockB+j2*actual_kc*Blocking::PacketSize,
@@ -244,17 +244,17 @@ struct ei_triangular_solve_matrix vertical panels of rhs)
- for (int j2 = IsLowerTriangular
+ for (int j2 = IsLower
? (actual_kc - ((actual_kc%SmallPanelWidth) ? (actual_kc%SmallPanelWidth)
: SmallPanelWidth))
: 0;
- IsLowerTriangular ? j2>=0 : j2=0 : j2(actual_kc-j2, SmallPanelWidth);
int absolute_j2 = actual_k2 + j2;
- int panelOffset = IsLowerTriangular ? j2+actualPanelWidth : 0;
- int panelLength = IsLowerTriangular ? actual_kc - j2 - actualPanelWidth : j2;
+ int panelOffset = IsLower ? j2+actualPanelWidth : 0;
+ int panelLength = IsLower ? actual_kc - j2 - actualPanelWidth : j2;
// GEBP
if(panelLength>0)
@@ -269,17 +269,17 @@ struct ei_triangular_solve_matrix class ColPivHouseholderQR;
template class FullPivHouseholderQR;
template class SVD;
template class JacobiSVD;
-template class LLT;
+template class LLT;
template class LDLT;
template class HouseholderSequence;
template class PlanarRotation;
diff --git a/Eigen/src/Core/util/XprHelper.h b/Eigen/src/Core/util/XprHelper.h
index 4bff09252..9b17a2b0e 100644
--- a/Eigen/src/Core/util/XprHelper.h
+++ b/Eigen/src/Core/util/XprHelper.h
@@ -263,8 +263,7 @@ template::ty
template struct ei_are_flags_consistent
{
- enum { ret = !( (Flags&UnitDiagBit && Flags&ZeroDiagBit) )
- };
+ enum { ret = true };
};
/** \internal Helper base class to add a scalar multiple operator
diff --git a/Eigen/src/Eigen2Support/Flagged.h b/Eigen/src/Eigen2Support/Flagged.h
index 94e7f9119..bed110b64 100644
--- a/Eigen/src/Eigen2Support/Flagged.h
+++ b/Eigen/src/Eigen2Support/Flagged.h
@@ -109,6 +109,12 @@ template clas
const ExpressionType& _expression() const { return m_matrix; }
+ template
+ typename ExpressionType::PlainMatrixType solveTriangular(const MatrixBase& other) const;
+
+ template
+ void solveTriangularInPlace(const MatrixBase& other) const;
+
protected:
ExpressionTypeNested m_matrix;
};
diff --git a/Eigen/src/Eigen2Support/TriangularSolver.h b/Eigen/src/Eigen2Support/TriangularSolver.h
index e69de29bb..94b92577e 100644
--- a/Eigen/src/Eigen2Support/TriangularSolver.h
+++ b/Eigen/src/Eigen2Support/TriangularSolver.h
@@ -0,0 +1,53 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2010 Gael Guennebaud
+//
+// 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 .
+
+#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
+template
+typename ExpressionType::PlainMatrixType
+Flagged::solveTriangular(const MatrixBase& other) const
+{
+ return m_matrix.template triangularView.solve(other.derived());
+}
+
+template
+template
+void Flagged::solveTriangularInPlace(const MatrixBase& other) const
+{
+ m_matrix.template triangularView.solveInPlace(other.derived());
+}
+
+#endif // EIGEN_TRIANGULAR_SOLVER2_H
diff --git a/Eigen/src/Eigenvalues/HessenbergDecomposition.h b/Eigen/src/Eigenvalues/HessenbergDecomposition.h
index 636b2f4f7..e3f64b807 100644
--- a/Eigen/src/Eigenvalues/HessenbergDecomposition.h
+++ b/Eigen/src/Eigenvalues/HessenbergDecomposition.h
@@ -199,7 +199,7 @@ HessenbergDecomposition::matrixH() const
int n = m_matrix.rows();
MatrixType matH = m_matrix;
if (n>2)
- matH.corner(BottomLeft,n-2, n-2).template triangularView().setZero();
+ matH.corner(BottomLeft,n-2, n-2).template triangularView().setZero();
return matH;
}
diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
index 336976517..960e9b417 100644
--- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
+++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
@@ -247,11 +247,11 @@ compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors
matC.adjointInPlace();
// this version works too:
// matC = matC.transpose();
-// cholB.matrixL().conjugate().template marked().solveTriangularInPlace(matC);
+// cholB.matrixL().conjugate().template marked().solveTriangularInPlace(matC);
// matC = matC.transpose();
// FIXME: this should work: (currently it only does for small matrices)
// Transpose trMatC(matC);
-// cholB.matrixL().conjugate().eval().template marked().solveTriangularInPlace(trMatC);
+// cholB.matrixL().conjugate().eval().template marked().solveTriangularInPlace(trMatC);
compute(matC, computeEigenvectors);
@@ -275,7 +275,7 @@ template
inline Matrix::Scalar>::Real, ei_traits