Added Triangular expression to extract upper or lower (strictly or not)

part of a matrix. Triangular also provide an optimised method for forward
and backward substitution. Further optimizations regarding assignments and
products might come later.

Updated determinant() to take into account triangular matrices.

Started the QR module with a QR decompostion algorithm.
Help needed to build a QR algorithm (eigen solver) based on it.

21 files changed, 586 insertions, 19 deletions

diff --git a/Eigen/Core b/Eigen/Core
index 3007899d1..e5ec8c08c 100644
--- a/Eigen/Core
+++ b/Eigen/Core
@@ -51,6 +51,7 @@ namespace Eigen {
 #include "src/Core/IO.h"
 #include "src/Core/Swap.h"
 #include "src/Core/CommaInitializer.h"
+#include "src/Core/Triangular.h"
 
 } // namespace Eigen
 
diff --git a/Eigen/QR b/Eigen/QR
new file mode 100644
index 000000000..2d166dcc6
--- /dev/null
+++ b/Eigen/QR
@@ -0,0 +1,12 @@
+#ifndef EIGEN_QR_MODULE_H
+#define EIGEN_QR_MODULE_H
+
+#include "Core"
+
+namespace Eigen {
+
+#include "Eigen/src/QR/QR.h"
+
+} // namespace Eigen
+
+#endif // EIGEN_QR_MODULE_H
diff --git a/Eigen/src/CMakeLists.txt b/Eigen/src/CMakeLists.txt
index e811772e5..ac05e77cd 100644
--- a/Eigen/src/CMakeLists.txt
+++ b/Eigen/src/CMakeLists.txt
@@ -1,2 +1,3 @@
 ADD_SUBDIRECTORY(Core)
 ADD_SUBDIRECTORY(LU)
+ADD_SUBDIRECTORY(QR)
diff --git a/Eigen/src/Core/Block.h b/Eigen/src/Core/Block.h
index 26e8fd50b..3b777a153 100644
--- a/Eigen/src/Core/Block.h
+++ b/Eigen/src/Core/Block.h
@@ -71,7 +71,7 @@ struct ei_traits<Block<MatrixType, BlockRows, BlockCols> >
                       || (ColsAtCompileTime != Dynamic && MatrixType::ColsAtCompileTime == Dynamic))
                       ? ~LargeBit
                       : ~(unsigned int)0,
-    Flags = MatrixType::Flags & FlagsMaskLargeBit & ~(VectorizableBit | Like1DArrayBit),
+    Flags = MatrixType::Flags & DefaultLostFlagMask & FlagsMaskLargeBit,
     CoeffReadCost = MatrixType::CoeffReadCost
   };
 };
diff --git a/Eigen/src/Core/CwiseBinaryOp.h b/Eigen/src/Core/CwiseBinaryOp.h
index 0cddc5a30..0032d22e2 100644
--- a/Eigen/src/Core/CwiseBinaryOp.h
+++ b/Eigen/src/Core/CwiseBinaryOp.h
@@ -60,9 +60,11 @@ struct ei_traits<CwiseBinaryOp<BinaryOp, Lhs, Rhs> >
     ColsAtCompileTime = Lhs::ColsAtCompileTime,
     MaxRowsAtCompileTime = Lhs::MaxRowsAtCompileTime,
     MaxColsAtCompileTime = Lhs::MaxColsAtCompileTime,
-    Flags = ((Lhs::Flags | Rhs::Flags) & ~VectorizableBit)
+    Flags = ((Lhs::Flags | Rhs::Flags) & (
+        DefaultLostFlagMask
+      | Like1DArrayBit
       | (ei_functor_traits<BinaryOp>::IsVectorizable && ((Lhs::Flags & RowMajorBit)==(Rhs::Flags & RowMajorBit))
-        ? (Lhs::Flags & Rhs::Flags & VectorizableBit) : 0),
+        ? VectorizableBit : 0))),
     CoeffReadCost = Lhs::CoeffReadCost + Rhs::CoeffReadCost + ei_functor_traits<BinaryOp>::Cost
   };
 };
diff --git a/Eigen/src/Core/CwiseNullaryOp.h b/Eigen/src/Core/CwiseNullaryOp.h
index 6a0b6bc0c..c93331039 100644
--- a/Eigen/src/Core/CwiseNullaryOp.h
+++ b/Eigen/src/Core/CwiseNullaryOp.h
@@ -48,7 +48,7 @@ struct ei_traits<CwiseNullaryOp<NullaryOp, MatrixType> >
     ColsAtCompileTime = MatrixType::ColsAtCompileTime,
     MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
     MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
-    Flags = (MatrixType::Flags & ~VectorizableBit)
+    Flags = (MatrixType::Flags & (DefaultLostFlagMask | Like1DArrayBit))
       | ei_functor_traits<NullaryOp>::IsVectorizable
       | (ei_functor_traits<NullaryOp>::IsRepeatable ? 0 : EvalBeforeNestingBit),
     CoeffReadCost = ei_functor_traits<NullaryOp>::Cost
diff --git a/Eigen/src/Core/CwiseUnaryOp.h b/Eigen/src/Core/CwiseUnaryOp.h
index 2a6858703..287a157c5 100644
--- a/Eigen/src/Core/CwiseUnaryOp.h
+++ b/Eigen/src/Core/CwiseUnaryOp.h
@@ -50,8 +50,9 @@ struct ei_traits<CwiseUnaryOp<UnaryOp, MatrixType> >
     ColsAtCompileTime = MatrixType::ColsAtCompileTime,
     MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
     MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
-    Flags = (MatrixType::Flags & ~VectorizableBit)
-      | (ei_functor_traits<UnaryOp>::IsVectorizable ? (MatrixType::Flags & VectorizableBit) : 0),
+    Flags = (MatrixType::Flags & (
+      DefaultLostFlagMask | Like1DArrayBit
+      | ei_functor_traits<UnaryOp>::IsVectorizable ? VectorizableBit : 0)),
     CoeffReadCost = MatrixType::CoeffReadCost + ei_functor_traits<UnaryOp>::Cost
   };
 };
diff --git a/Eigen/src/Core/DiagonalCoeffs.h b/Eigen/src/Core/DiagonalCoeffs.h
index dc35de4d8..8cdb4735e 100644
--- a/Eigen/src/Core/DiagonalCoeffs.h
+++ b/Eigen/src/Core/DiagonalCoeffs.h
@@ -54,7 +54,7 @@ struct ei_traits<DiagonalCoeffs<MatrixType> >
     MaxColsAtCompileTime = 1,
     Flags = (RowsAtCompileTime == Dynamic && ColsAtCompileTime == Dynamic
             ? (unsigned int)MatrixType::Flags
-            : (unsigned int)MatrixType::Flags &~ LargeBit) & ~(VectorizableBit | Like1DArrayBit),
+            : (unsigned int)MatrixType::Flags &~ LargeBit) & DefaultLostFlagMask,
     CoeffReadCost = MatrixType::CoeffReadCost
   };
 };
diff --git a/Eigen/src/Core/DiagonalMatrix.h b/Eigen/src/Core/DiagonalMatrix.h
index 2c8f6cdf7..94c7d2b88 100644
--- a/Eigen/src/Core/DiagonalMatrix.h
+++ b/Eigen/src/Core/DiagonalMatrix.h
@@ -47,7 +47,7 @@ struct ei_traits<DiagonalMatrix<CoeffsVectorType> >
     ColsAtCompileTime = CoeffsVectorType::SizeAtCompileTime,
     MaxRowsAtCompileTime = CoeffsVectorType::MaxSizeAtCompileTime,
     MaxColsAtCompileTime = CoeffsVectorType::MaxSizeAtCompileTime,
-    Flags = CoeffsVectorType::Flags & ~(VectorizableBit | Like1DArrayBit),
+    Flags = CoeffsVectorType::Flags & DefaultLostFlagMask,
     CoeffReadCost = CoeffsVectorType::CoeffReadCost
   };
 };
diff --git a/Eigen/src/Core/Map.h b/Eigen/src/Core/Map.h
index f17107a65..266da0cfe 100644
--- a/Eigen/src/Core/Map.h
+++ b/Eigen/src/Core/Map.h
@@ -47,7 +47,7 @@ struct ei_traits<Map<MatrixType> >
     ColsAtCompileTime = MatrixType::ColsAtCompileTime,
     MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
     MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
-    Flags = MatrixType::Flags & ~VectorizableBit,
+    Flags = MatrixType::Flags & DefaultLostFlagMask,
     CoeffReadCost = NumTraits<Scalar>::ReadCost
   };
 };
diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h
index bed6b2f93..a60ead6ba 100644
--- a/Eigen/src/Core/MatrixBase.h
+++ b/Eigen/src/Core/MatrixBase.h
@@ -519,6 +519,22 @@ template<typename Derived> class MatrixBase
     { return *static_cast<Derived*>(const_cast<MatrixBase*>(this)); }
     //@}
 
+    /// \name Triangular matrices
+    //@{
+    Triangular<Upper, Derived> upper(void);
+    const Triangular<Upper, Derived> upper(void) const;
+    const Triangular<Upper|UnitDiagBit, Derived> upperWithUnitDiag(void) const;
+    const Triangular<Upper|NullDiagBit, Derived> upperWithNullDiag(void) const;
+
+    Triangular<Lower, Derived> lower(void);
+    const Triangular<Lower, Derived> lower(void) const;
+    const Triangular<Lower|UnitDiagBit, Derived> lowerWithUnitDiag(void) const;
+    const Triangular<Lower|NullDiagBit, Derived> lowerWithNullDiag(void) const;
+
+    bool isUpper(RealScalar prec = precision<Scalar>()) const;
+    bool isLower(RealScalar prec = precision<Scalar>()) const;
+    //@}
+
     /** \name LU module
       *
       * \code #include <Eigen/LU> \endcode
@@ -529,6 +545,14 @@ template<typename Derived> class MatrixBase
     Scalar determinant() const;
     //@}
 
+    /** \name QR module
+      *
+      * \code #include <Eigen/QR> \endcode
+      */
+    //@{
+    const QR<typename ei_eval<Derived>::type> qr() const;
+    //@}
+
   private:
 
     PacketScalar _packetCoeff(int , int) const { ei_internal_assert(false && "_packetCoeff not defined"); }
diff --git a/Eigen/src/Core/Minor.h b/Eigen/src/Core/Minor.h
index 5c9afc162..b44dca6e5 100644
--- a/Eigen/src/Core/Minor.h
+++ b/Eigen/src/Core/Minor.h
@@ -50,7 +50,7 @@ struct ei_traits<Minor<MatrixType> >
                                 MatrixType::MaxRowsAtCompileTime - 1 : Dynamic,
     MaxColsAtCompileTime = (MatrixType::MaxColsAtCompileTime != Dynamic) ?
                                 MatrixType::MaxColsAtCompileTime - 1 : Dynamic,
-    Flags = MatrixType::Flags & ~VectorizableBit,
+    Flags = MatrixType::Flags & DefaultLostFlagMask,
     CoeffReadCost = MatrixType::CoeffReadCost
   };
 };
diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h
index b77fd3eae..a60c13e45 100644
--- a/Eigen/src/Core/Product.h
+++ b/Eigen/src/Core/Product.h
@@ -135,7 +135,7 @@ struct ei_traits<Product<Lhs, Rhs, EvalMode> >
           | EvalBeforeAssigningBit
           | (ei_product_eval_mode<Lhs, Rhs>::value == (int)CacheOptimalProduct ? EvalBeforeNestingBit : 0))
           & (
-              ~(RowMajorBit | VectorizableBit | Like1DArrayBit)
+              DefaultLostFlagMask & (~RowMajorBit)
               | (
                   (
                     (!(Lhs::Flags & RowMajorBit)) && (Lhs::Flags & VectorizableBit)
diff --git a/Eigen/src/Core/Redux.h b/Eigen/src/Core/Redux.h
index c9332c847..cb8f4d525 100644
--- a/Eigen/src/Core/Redux.h
+++ b/Eigen/src/Core/Redux.h
@@ -96,7 +96,7 @@ struct ei_traits<PartialRedux<Direction, BinaryOp, MatrixType> >
     MaxColsAtCompileTime = Direction==Horizontal ? 1 : MatrixType::MaxColsAtCompileTime,
     Flags = ((RowsAtCompileTime == Dynamic || ColsAtCompileTime == Dynamic)
           ? (unsigned int)_MatrixTypeNested::Flags
-          : (unsigned int)_MatrixTypeNested::Flags & ~LargeBit) & ~(VectorizableBit | Like1DArrayBit),
+          : (unsigned int)_MatrixTypeNested::Flags & ~LargeBit) & DefaultLostFlagMask,
     TraversalSize = Direction==Vertical ? RowsAtCompileTime : ColsAtCompileTime,
     CoeffReadCost = TraversalSize * _MatrixTypeNested::CoeffReadCost
                   + (TraversalSize - 1) * ei_functor_traits<BinaryOp>::Cost
diff --git a/Eigen/src/Core/Triangular.h b/Eigen/src/Core/Triangular.h
new file mode 100755
index 000000000..9bc8def38
--- /dev/null
+++ b/Eigen/src/Core/Triangular.h
@@ -0,0 +1,347 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_TRIANGULAR_H
+#define EIGEN_TRIANGULAR_H
+
+/** \class Triangular
+  *
+  * \brief Expression of a triangular matrix from a squared matrix
+  *
+  * \param Mode or-ed bit field indicating the triangular part (Upper or Lower) we are taking,
+  * and the property of the diagonal if any (UnitDiagBit or NullDiagBit).
+  * \param MatrixType the type of the object in which we are taking the triangular part
+  *
+  * This class represents an expression of the upper or lower triangular part of
+  * a squared matrix. It is the return type of MatrixBase::upper(), MatrixBase::lower(),
+  * MatrixBase::upperWithUnitDiagBit(), etc., and used to optimize operations involving
+  * triangular matrices. Most of the time this is the only way it is used.
+  *
+  * Examples of some key features:
+  * \code
+  * m1 = (<any expression>).upper();
+  * \endcode
+  * In this example, the strictly lower part of the expression is not evaluated,
+  * m1 might be resized and the strict lower part of m1 == 0.
+  *
+  * \code
+  * m1.upper() = <any expression>;
+  * \endcode
+  * This example diverge from the previous one in the sense that the strictly
+  * lower part of m1 is left unchanged, and optimal loops are employed. Note that
+  * m1 might also be resized.
+  *
+  * Of course, in both examples \c <any \c expression> has to be a squared matrix.
+  *
+  * \sa MatrixBase::upper(), MatrixBase::lower(), class TriangularProduct
+  */
+template<int Mode, typename MatrixType>
+struct ei_traits<Triangular<Mode, MatrixType> >
+{
+  typedef typename MatrixType::Scalar Scalar;
+  enum {
+    RowsAtCompileTime = MatrixType::SizeAtCompileTime,
+    ColsAtCompileTime = MatrixType::SizeAtCompileTime,
+    MaxRowsAtCompileTime = MatrixType::MaxSizeAtCompileTime,
+    MaxColsAtCompileTime = MatrixType::MaxSizeAtCompileTime,
+    Flags = MatrixType::Flags & (~(VectorizableBit | Like1DArrayBit)) | Mode,
+    CoeffReadCost = MatrixType::CoeffReadCost
+  };
+};
+
+template<int Mode, typename MatrixType> class Triangular
+  : public MatrixBase<Triangular<Mode,MatrixType> >
+{
+  public:
+
+    EIGEN_GENERIC_PUBLIC_INTERFACE(Triangular)
+
+    Triangular(const MatrixType& matrix)
+      : m_matrix(matrix)
+    {
+      assert(!( (Flags&UnitDiagBit) && (Flags&NullDiagBit)));
+      assert(matrix.rows()==matrix.cols());
+    }
+
+    /** Overloaded to keep a Triangular expression */
+    Triangular<(Upper | Lower) xor Mode, Transpose<MatrixType> > transpose()
+    {
+      return Triangular<(Upper | Lower) xor Mode, Transpose<MatrixType> >((m_matrix.transpose()));
+    }
+
+    /** Overloaded to keep a Triangular expression */
+    const Triangular<(Upper | Lower) xor Mode, Transpose<MatrixType> > transpose() const
+    {
+      return Triangular<(Upper | Lower) xor Mode, Transpose<MatrixType> >((m_matrix.transpose()));
+    }
+
+#if 0
+
+    template<typename OtherDerived>
+    Triangular& operator=(const MatrixBase<OtherDerived>& other);
+
+    /** Overloaded to provide optimal evaluation loops */
+    template<typename OtherDerived>
+    Triangular& operator +=(const MatrixBase<OtherDerived>& other)
+    {
+        return *this = m_matrix + other;
+    }
+
+    /** Overloaded to provide optimal evaluation loops */
+    template<typename OtherDerived>
+    Triangular& operator *=(const MatrixBase<OtherDerived>& other)
+    {
+        return *this = this->lazyProduct(other).eval();
+    }
+
+    /** Optimized triangular matrix - matrix product */
+    template<typename OtherDerived>
+    TriangularProduct<Mode, MatrixType, OtherDerived> lazyProduct(const MatrixBase<Scalar, OtherDerived>& other) const
+    {
+      return TriangularProduct<Mode,MatrixType,OtherDerived>(m_matrix, other.ref());
+    }
+
+    /** Optimized triangular matrix - matrix product */
+    template<typename OtherDerived>
+    Eval<TriangularProduct<Mode, MatrixType, OtherDerived> > operator * (const MatrixBase<Scalar, OtherDerived>& other) const
+    {
+      return this->lazyProduct(other).eval();
+    }
+
+    /** Optimized matrix - triangular matrix product */
+    template<typename OtherDerived>
+    friend Eval<Transpose<TriangularProduct<0x1 xor Mode, Transpose<MatRef>, Transpose<OtherDerived> > > >
+    operator * (const MatrixBase<Scalar, OtherDerived>& other, const Triangular<Mode,MatrixType>& tri)
+    {
+      return tri.transpose().lazyProduct(other.transpose()).transpose().eval();
+    }
+
+#endif
+
+    /** \returns the product of the inverse of *this with \a other.
+      *
+      * This function computes the inverse-matrix matrix product inv(*this) \a other
+      * This process is also as forward (resp. backward) substitution if *this is an upper (resp. lower)
+      * triangular matrix.
+      */
+    template<typename OtherDerived>
+    typename OtherDerived::Eval inverseProduct(const MatrixBase<OtherDerived>& other) const
+    {
+      assert(_cols() == other.rows());
+      assert(!(Flags & NullDiagBit));
+
+      typename OtherDerived::Eval res(other.rows(), other.cols());
+
+      for (int c=0 ; c<other.cols() ; ++c)
+      {
+        if (Flags & Lower)
+        {
+          // forward substitution
+          if (Flags & UnitDiagBit)
+            res.col(c)[0] = other.col(c)[0];
+          else
+            res.col(c)[0] = other.col(c)[0]/_coeff(0, 0);
+          for (int i=1 ; i<_rows() ; ++i)
+          {
+            Scalar tmp = other.col(c)[i];
+            for (int j = 0 ; j < i ; ++j)
+              tmp -= _coeff(i,j) * res.col(c)[j];
+            if (Flags & UnitDiagBit)
+              res.col(c)[i] = tmp;
+            else
+              res.col(c)[i] = tmp/_coeff(i,i);
+          }
+        }
+        else
+        {
+          // backward substitution
+          if (Flags & UnitDiagBit)
+            res.col(c)[_cols()-1] = other.col(c)[_cols()-1];
+          else
+            res.col(c)[_cols()-1] = other.col(c)[_cols()-1]/_coeff(_rows()-1, _cols()-1);
+          for (int i=_rows()-2 ; i>=0 ; --i)
+          {
+            Scalar tmp = other.col(c)[i];
+            for (int j = i+1 ; j < _cols() ; ++j)
+              tmp -= _coeff(i,j) * res.col(c)[j];
+            if (Flags & UnitDiagBit)
+              res.col(c)[i] = tmp;
+            else
+              res.col(c)[i] = tmp/_coeff(i,i);
+          }
+        }
+      }
+      return res;
+    }
+
+  private:
+
+    int _rows() const { return m_matrix.rows(); }
+    int _cols() const { return m_matrix.cols(); }
+
+    Scalar& _coeffRef(int row, int col)
+    {
+      assert( ((! Flags & Lower) && row<=col) || (Flags & Lower && col<=row));
+      return m_matrix.coeffRef(row, col);
+    }
+
+    Scalar _coeff(int row, int col) const
+    {
+      if ((Flags & Lower) ? col>row : row>col)
+        return 0;
+      if (Flags & UnitDiagBit)
+        return col==row ? 1 : m_matrix.coeff(row, col);
+      else if (Flags & NullDiagBit)
+        return col==row ? 0 : m_matrix.coeff(row, col);
+      else
+        return m_matrix.coeff(row, col);
+    }
+
+  protected:
+
+    const typename MatrixType::Nested m_matrix;
+};
+
+/** \returns an expression of a upper triangular matrix
+  *
+  * \sa isUpper(), upperWithNullDiagBit(), upperWithNullDiagBit(), lower()
+  */
+template<typename Derived>
+Triangular<Upper, Derived> MatrixBase<Derived>::upper(void)
+{
+  return Triangular<Upper,Derived>(derived());
+}
+
+/** This is the const version of upper(). */
+template<typename Derived>
+const Triangular<Upper, Derived> MatrixBase<Derived>::upper(void) const
+{
+  return Triangular<Upper,Derived>(derived());
+}
+
+/** \returns an expression of a lower triangular matrix
+  *
+  * \sa isLower(), lowerWithUnitDiag(), lowerWithNullDiag(), upper()
+  */
+template<typename Derived>
+Triangular<Lower, Derived> MatrixBase<Derived>::lower(void)
+{
+  return Triangular<Lower,Derived>(derived());
+}
+
+/** This is the const version of lower().*/
+template<typename Derived>
+const Triangular<Lower, Derived> MatrixBase<Derived>::lower(void) const
+{
+  return Triangular<Lower,Derived>(derived());
+}
+
+/** \returns an expression of a upper triangular matrix with a unit diagonal
+  *
+  * \sa upper(), lowerWithUnitDiagBit()
+  */
+template<typename Derived>
+const Triangular<Upper|UnitDiagBit, Derived> MatrixBase<Derived>::upperWithUnitDiag(void) const
+{
+  return Triangular<Upper|UnitDiagBit, Derived>(derived());
+}
+
+/** \returns an expression of a strictly upper triangular matrix (diagonal==zero)
+  * FIXME could also be called strictlyUpper() or upperStrict()
+  *
+  * \sa upper(), lowerWithNullDiag()
+  */
+template<typename Derived>
+const Triangular<Upper|NullDiagBit, Derived> MatrixBase<Derived>::upperWithNullDiag(void) const
+{
+  return Triangular<Upper|NullDiagBit, Derived>(derived());
+}
+
+/** \returns an expression of a lower triangular matrix with a unit diagonal
+  *
+  * \sa lower(), upperWithUnitDiag()
+  */
+template<typename Derived>
+const Triangular<Lower|UnitDiagBit, Derived> MatrixBase<Derived>::lowerWithUnitDiag(void) const
+{
+  return Triangular<Lower|UnitDiagBit, Derived>(derived());
+}
+
+/** \returns an expression of a strictly lower triangular matrix (diagonal==zero)
+  * FIXME could also be called strictlyLower() or lowerStrict()
+  *
+  * \sa lower(), upperWithNullDiag()
+  */
+template<typename Derived>
+const Triangular<Lower|NullDiagBit, Derived> MatrixBase<Derived>::lowerWithNullDiag(void) const
+{
+  return Triangular<Lower|NullDiagBit, Derived>(derived());
+}
+
+/** \returns true if *this is approximately equal to an upper triangular matrix,
+  *          within the precision given by \a prec.
+  *
+  * \sa isLower(), upper()
+  */
+template<typename Derived>
+bool MatrixBase<Derived>::isUpper(RealScalar prec) const
+{
+  if(cols() != rows()) return false;
+  RealScalar maxAbsOnUpperPart = static_cast<RealScalar>(-1);
+  for(int j = 0; j < cols(); j++)
+    for(int i = 0; i <= j; i++)
+    {
+      RealScalar absValue = ei_abs(coeff(i,j));
+      if(absValue > maxAbsOnUpperPart) maxAbsOnUpperPart = absValue;
+    }
+  for(int j = 0; j < cols()-1; j++)
+    for(int i = j+1; i < rows(); i++)
+      if(!ei_isMuchSmallerThan(coeff(i, j), maxAbsOnUpperPart, prec)) return false;
+  return true;
+}
+
+/** \returns true if *this is approximately equal to a lower triangular matrix,
+  *          within the precision given by \a prec.
+  *
+  * \sa isUpper(), upper()
+  */
+template<typename Derived>
+bool MatrixBase<Derived>::isLower(RealScalar prec) const
+{
+  if(cols() != rows()) return false;
+  RealScalar maxAbsOnLowerPart = static_cast<RealScalar>(-1);
+  for(int j = 0; j < cols(); j++)
+    for(int i = j; i < rows(); i++)
+    {
+      RealScalar absValue = ei_abs(coeff(i,j));
+      if(absValue > maxAbsOnLowerPart) maxAbsOnLowerPart = absValue;
+    }
+  for(int j = 1; j < cols(); j++)
+    for(int i = 0; i < j; i++)
+      if(!ei_isMuchSmallerThan(coeff(i, j), maxAbsOnLowerPart, prec)) return false;
+  return true;
+}
+
+
+#endif // EIGEN_TRIANGULAR_H
diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h
index 97d8c346a..ae0451156 100644
--- a/Eigen/src/Core/util/Constants.h
+++ b/Eigen/src/Core/util/Constants.h
@@ -30,15 +30,24 @@ const int Dynamic = 10000;
 
 // matrix/expression flags
 const unsigned int RowMajorBit = 0x1;
-const unsigned int EvalBeforeNestingBit = 0x2;
-const unsigned int EvalBeforeAssigningBit = 0x4;
+const unsigned int EvalBeforeNestingBit = 0x2;  ///< means the expression should be evaluated by the calling expression
+const unsigned int EvalBeforeAssigningBit = 0x4;///< means the expression should be evaluated before any assignement
 const unsigned int LargeBit = 0x8;
 #ifdef EIGEN_VECTORIZE
-const unsigned int VectorizableBit = 0x10;
+const unsigned int VectorizableBit = 0x10;  ///< means the expression might be vectorized
 #else
 const unsigned int VectorizableBit = 0x0;
 #endif
-const unsigned int Like1DArrayBit = 0x20;
+const unsigned int Like1DArrayBit = 0x20;   ///< means the expression can be seen as 1D vector (used for explicit vectorization)
+const unsigned int NullDiagBit = 0x40;      ///< means all diagonal coefficients are equal to 0
+const unsigned int UnitDiagBit = 0x80;      ///< means all diagonal coefficients are equal to 1
+const unsigned int NullLowerBit = 0x200;    ///< means the strictly triangular lower part is 0
+const unsigned int NullUpperBit = 0x400;    ///< means the strictly triangular upper part is 0
+
+enum { Upper=NullLowerBit, Lower=NullUpperBit };
+
+// list of flags that are lost by default
+const unsigned int DefaultLostFlagMask = ~(VectorizableBit | Like1DArrayBit | NullDiagBit | UnitDiagBit | NullLowerBit | NullUpperBit);
 
 enum { ConditionalJumpCost = 5 };
 enum CornerType { TopLeft, TopRight, BottomLeft, BottomRight };
diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h
index adbbdeff8..c18977820 100644
--- a/Eigen/src/Core/util/ForwardDeclarations.h
+++ b/Eigen/src/Core/util/ForwardDeclarations.h
@@ -47,8 +47,9 @@ template<typename Lhs, typename Rhs, int EvalMode=ei_product_eval_mode<Lhs,Rhs>:
 template<typename CoeffsVectorType> class DiagonalMatrix;
 template<typename MatrixType> class DiagonalCoeffs;
 template<typename MatrixType> class Map;
-template<typename Derived> class Eval;
+// template<typename Derived> class Eval;
 template<int Direction, typename UnaryOp, typename MatrixType> class PartialRedux;
+template<int Mode, typename MatrixType> class Triangular;
 
 template<typename Scalar> struct ei_scalar_sum_op;
 template<typename Scalar> struct ei_scalar_difference_op;
@@ -71,5 +72,6 @@ template<typename Scalar> struct ei_scalar_min_op;
 template<typename Scalar> struct ei_scalar_max_op;
 
 template<typename ExpressionType, bool CheckExistence = true> class Inverse;
+template<typename MatrixType> class QR;
 
 #endif // EIGEN_FORWARDDECLARATIONS_H
diff --git a/Eigen/src/LU/Determinant.h b/Eigen/src/LU/Determinant.h
index d3998195c..1675efc10 100644
--- a/Eigen/src/LU/Determinant.h
+++ b/Eigen/src/LU/Determinant.h
@@ -71,7 +71,16 @@ template<typename Derived>
 typename ei_traits<Derived>::Scalar MatrixBase<Derived>::determinant() const
 {
   assert(rows() == cols());
-  if(rows() <= 4) return ei_bruteforce_det(derived());
+  if (Derived::Flags & (NullLowerBit | NullUpperBit))
+  {
+    if (Derived::Flags & UnitDiagBit)
+      return 1;
+    else if (Derived::Flags & NullDiagBit)
+      return 0;
+    else
+      return derived().diagonal().redux(ei_scalar_product_op<Scalar>());
+  }
+  else if(rows() <= 4) return ei_bruteforce_det(derived());
   else assert(false); // unimplemented for now
 }
 
diff --git a/Eigen/src/LU/Inverse.h b/Eigen/src/LU/Inverse.h
index 00ffa25c8..a9c0ead96 100644
--- a/Eigen/src/LU/Inverse.h
+++ b/Eigen/src/LU/Inverse.h
@@ -98,7 +98,7 @@ template<typename MatrixType, bool CheckExistence> class Inverse : ei_no_assignm
 
   protected:
     bool m_exists;
-    typename MatrixType::Eval m_inverse;
+    MatrixType m_inverse;
 };
 
 template<typename MatrixType, bool CheckExistence>
diff --git a/Eigen/src/QR/CMakeLists.txt b/Eigen/src/QR/CMakeLists.txt
new file mode 100644
index 000000000..fca336ae6
--- /dev/null
+++ b/Eigen/src/QR/CMakeLists.txt
@@ -0,0 +1,6 @@
+FILE(GLOB Eigen_QR_SRCS "*.h")
+
+INSTALL(FILES
+  ${Eigen_QR_SRCS}
+  DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/QR
+  )
diff --git a/Eigen/src/QR/QR.h b/Eigen/src/QR/QR.h
new file mode 100644
index 000000000..d0121cc7a
--- /dev/null
+++ b/Eigen/src/QR/QR.h
@@ -0,0 +1,153 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_QR_H
+#define EIGEN_QR_H
+
+/** \class QR
+  *
+  * \brief QR decomposition of a matrix
+  *
+  * \param MatrixType the type of the matrix of which we are computing the QR decomposition
+  *
+  * This class performs a QR decomposition using Householder transformations. The result is
+  * stored in a compact way.
+  *
+  * \todo add convenient method to direclty use the result in a compact way. First need to determine
+  * typical use cases though.
+  *
+  * \todo what about complex matrices ?
+  *
+  * \sa MatrixBase::qr()
+  */
+template<typename MatrixType> class QR
+{
+  public:
+
+    typedef typename MatrixType::Scalar Scalar;
+    typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> RMatrixType;
+    typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
+
+    QR(const MatrixType& matrix)
+      : m_qr(matrix.rows(), matrix.cols()),
+        m_norms(matrix.cols())
+    {
+      _compute(matrix);
+    }
+
+    /** \returns whether or not the matrix is of full rank */
+    bool isFullRank() const { return ei_isMuchSmallerThan(m_norms.cwiseAbs().minCoeff(), Scalar(1)); }
+
+    RMatrixType matrixR(void) const;
+
+    MatrixType matrixQ(void) const;
+
+  private:
+
+    void _compute(const MatrixType& matrix);
+
+  protected:
+    MatrixType m_qr;
+    VectorType m_norms;
+};
+
+template<typename MatrixType>
+void QR<MatrixType>::_compute(const MatrixType& matrix)
+{
+  m_qr = matrix;
+  int rows = matrix.rows();
+  int cols = matrix.cols();
+
+  for (int k = 0; k < cols; k++)
+  {
+    int remainingSize = rows-k;
+
+    Scalar nrm = m_qr.col(k).end(remainingSize).norm();
+
+    if (nrm != Scalar(0))
+    {
+      // form k-th Householder vector
+      if (m_qr(k,k) < 0)
+        nrm = -nrm;
+
+      m_qr.col(k).end(rows-k) /= nrm;
+      m_qr(k,k) += 1.0;
+
+      // apply transformation to remaining columns
+      for (int j = k+1; j < cols; j++)
+      {
+        Scalar s = -(m_qr.col(k).end(remainingSize).transpose() * m_qr.col(j).end(remainingSize))(0,0) / m_qr(k,k);
+        m_qr.col(j).end(remainingSize) += s * m_qr.col(k).end(remainingSize);
+      }
+    }
+    m_norms[k] = -nrm;
+  }
+}
+
+/** \returns the matrix R */
+template<typename MatrixType>
+typename QR<MatrixType>::RMatrixType QR<MatrixType>::matrixR(void) const
+{
+  int cols = m_qr.cols();
+  RMatrixType res = m_qr.block(0,0,cols,cols).upperWithNullDiag();
+  res.diagonal() = m_norms;
+  return res;
+}
+
+/** \returns the matrix Q */
+template<typename MatrixType>
+MatrixType QR<MatrixType>::matrixQ(void) const
+{
+  int rows = m_qr.rows();
+  int cols = m_qr.cols();
+  MatrixType res = MatrixType::identity(rows, cols);
+  for (int k = cols-1; k >= 0; k--)
+  {
+    for (int j = k; j < cols; j++)
+    {
+      if (res(k,k) != Scalar(0))
+      {
+        int endLength = rows-k;
+        Scalar s = -(m_qr.col(k).end(endLength).transpose() * res.col(j).end(endLength))(0,0) / m_qr(k,k);
+
+        res.col(j).end(endLength) += s * m_qr.col(k).end(endLength);
+      }
+    }
+  }
+  return res;
+}
+
+/** \return the QR decomposition of \c *this.
+  *
+  * \sa class QR
+  */
+template<typename Derived>
+const QR<typename ei_eval<Derived>::type>
+MatrixBase<Derived>::qr() const
+{
+  return QR<typename ei_eval<Derived>::type>(derived());
+}
+
+
+#endif // EIGEN_QR_H