aboutsummaryrefslogtreecommitdiffhomepage
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Eigen/Core1
-rw-r--r--Eigen/QR12
-rw-r--r--Eigen/src/CMakeLists.txt1
-rw-r--r--Eigen/src/Core/Block.h2
-rw-r--r--Eigen/src/Core/CwiseBinaryOp.h6
-rw-r--r--Eigen/src/Core/CwiseNullaryOp.h2
-rw-r--r--Eigen/src/Core/CwiseUnaryOp.h5
-rw-r--r--Eigen/src/Core/DiagonalCoeffs.h2
-rw-r--r--Eigen/src/Core/DiagonalMatrix.h2
-rw-r--r--Eigen/src/Core/Map.h2
-rw-r--r--Eigen/src/Core/MatrixBase.h24
-rw-r--r--Eigen/src/Core/Minor.h2
-rw-r--r--Eigen/src/Core/Product.h2
-rw-r--r--Eigen/src/Core/Redux.h2
-rwxr-xr-xEigen/src/Core/Triangular.h347
-rw-r--r--Eigen/src/Core/util/Constants.h17
-rw-r--r--Eigen/src/Core/util/ForwardDeclarations.h4
-rw-r--r--Eigen/src/LU/Determinant.h11
-rw-r--r--Eigen/src/LU/Inverse.h2
-rw-r--r--Eigen/src/QR/CMakeLists.txt6
-rw-r--r--Eigen/src/QR/QR.h153
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
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-06 09:20:57 +0000