From 4c92150676dea30748215549ba1b94df2624e652 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Sat, 26 Apr 2008 18:26:05 +0000 Subject: [PATCH] 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. --- Eigen/Core | 1 + Eigen/QR | 12 + Eigen/src/CMakeLists.txt | 1 + Eigen/src/Core/Block.h | 2 +- Eigen/src/Core/CwiseBinaryOp.h | 6 +- Eigen/src/Core/CwiseNullaryOp.h | 2 +- Eigen/src/Core/CwiseUnaryOp.h | 5 +- Eigen/src/Core/DiagonalCoeffs.h | 2 +- Eigen/src/Core/DiagonalMatrix.h | 2 +- Eigen/src/Core/Map.h | 2 +- Eigen/src/Core/MatrixBase.h | 24 ++ Eigen/src/Core/Minor.h | 2 +- Eigen/src/Core/Product.h | 2 +- Eigen/src/Core/Redux.h | 2 +- Eigen/src/Core/Triangular.h | 347 ++++++++++++++++++++++ Eigen/src/Core/util/Constants.h | 17 +- Eigen/src/Core/util/ForwardDeclarations.h | 4 +- Eigen/src/LU/Determinant.h | 11 +- Eigen/src/LU/Inverse.h | 2 +- Eigen/src/QR/CMakeLists.txt | 6 + Eigen/src/QR/QR.h | 153 ++++++++++ 21 files changed, 586 insertions(+), 19 deletions(-) create mode 100644 Eigen/QR create mode 100755 Eigen/src/Core/Triangular.h create mode 100644 Eigen/src/QR/CMakeLists.txt create mode 100644 Eigen/src/QR/QR.h 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 > || (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 > 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::IsVectorizable && ((Lhs::Flags & RowMajorBit)==(Rhs::Flags & RowMajorBit)) - ? (Lhs::Flags & Rhs::Flags & VectorizableBit) : 0), + ? VectorizableBit : 0))), CoeffReadCost = Lhs::CoeffReadCost + Rhs::CoeffReadCost + ei_functor_traits::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 > ColsAtCompileTime = MatrixType::ColsAtCompileTime, MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, - Flags = (MatrixType::Flags & ~VectorizableBit) + Flags = (MatrixType::Flags & (DefaultLostFlagMask | Like1DArrayBit)) | ei_functor_traits::IsVectorizable | (ei_functor_traits::IsRepeatable ? 0 : EvalBeforeNestingBit), CoeffReadCost = ei_functor_traits::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 > ColsAtCompileTime = MatrixType::ColsAtCompileTime, MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, - Flags = (MatrixType::Flags & ~VectorizableBit) - | (ei_functor_traits::IsVectorizable ? (MatrixType::Flags & VectorizableBit) : 0), + Flags = (MatrixType::Flags & ( + DefaultLostFlagMask | Like1DArrayBit + | ei_functor_traits::IsVectorizable ? VectorizableBit : 0)), CoeffReadCost = MatrixType::CoeffReadCost + ei_functor_traits::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 > 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 > 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 > ColsAtCompileTime = MatrixType::ColsAtCompileTime, MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, - Flags = MatrixType::Flags & ~VectorizableBit, + Flags = MatrixType::Flags & DefaultLostFlagMask, CoeffReadCost = NumTraits::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 class MatrixBase { return *static_cast(const_cast(this)); } //@} + /// \name Triangular matrices + //@{ + Triangular upper(void); + const Triangular upper(void) const; + const Triangular upperWithUnitDiag(void) const; + const Triangular upperWithNullDiag(void) const; + + Triangular lower(void); + const Triangular lower(void) const; + const Triangular lowerWithUnitDiag(void) const; + const Triangular lowerWithNullDiag(void) const; + + bool isUpper(RealScalar prec = precision()) const; + bool isLower(RealScalar prec = precision()) const; + //@} + /** \name LU module * * \code #include \endcode @@ -529,6 +545,14 @@ template class MatrixBase Scalar determinant() const; //@} + /** \name QR module + * + * \code #include \endcode + */ + //@{ + const QR::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 > 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 > | EvalBeforeAssigningBit | (ei_product_eval_mode::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 > 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::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 +// +// 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_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 = ().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() = ; + * \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 has to be a squared matrix. + * + * \sa MatrixBase::upper(), MatrixBase::lower(), class TriangularProduct + */ +template +struct ei_traits > +{ + 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 class Triangular + : public MatrixBase > +{ + 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 > transpose() + { + return Triangular<(Upper | Lower) xor Mode, Transpose >((m_matrix.transpose())); + } + + /** Overloaded to keep a Triangular expression */ + const Triangular<(Upper | Lower) xor Mode, Transpose > transpose() const + { + return Triangular<(Upper | Lower) xor Mode, Transpose >((m_matrix.transpose())); + } + +#if 0 + + template + Triangular& operator=(const MatrixBase& other); + + /** Overloaded to provide optimal evaluation loops */ + template + Triangular& operator +=(const MatrixBase& other) + { + return *this = m_matrix + other; + } + + /** Overloaded to provide optimal evaluation loops */ + template + Triangular& operator *=(const MatrixBase& other) + { + return *this = this->lazyProduct(other).eval(); + } + + /** Optimized triangular matrix - matrix product */ + template + TriangularProduct lazyProduct(const MatrixBase& other) const + { + return TriangularProduct(m_matrix, other.ref()); + } + + /** Optimized triangular matrix - matrix product */ + template + Eval > operator * (const MatrixBase& other) const + { + return this->lazyProduct(other).eval(); + } + + /** Optimized matrix - triangular matrix product */ + template + friend Eval, Transpose > > > + operator * (const MatrixBase& other, const Triangular& 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::Eval inverseProduct(const MatrixBase& other) const + { + assert(_cols() == other.rows()); + assert(!(Flags & NullDiagBit)); + + typename OtherDerived::Eval res(other.rows(), other.cols()); + + for (int c=0 ; c=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 +Triangular MatrixBase::upper(void) +{ + return Triangular(derived()); +} + +/** This is the const version of upper(). */ +template +const Triangular MatrixBase::upper(void) const +{ + return Triangular(derived()); +} + +/** \returns an expression of a lower triangular matrix + * + * \sa isLower(), lowerWithUnitDiag(), lowerWithNullDiag(), upper() + */ +template +Triangular MatrixBase::lower(void) +{ + return Triangular(derived()); +} + +/** This is the const version of lower().*/ +template +const Triangular MatrixBase::lower(void) const +{ + return Triangular(derived()); +} + +/** \returns an expression of a upper triangular matrix with a unit diagonal + * + * \sa upper(), lowerWithUnitDiagBit() + */ +template +const Triangular MatrixBase::upperWithUnitDiag(void) const +{ + return Triangular(derived()); +} + +/** \returns an expression of a strictly upper triangular matrix (diagonal==zero) + * FIXME could also be called strictlyUpper() or upperStrict() + * + * \sa upper(), lowerWithNullDiag() + */ +template +const Triangular MatrixBase::upperWithNullDiag(void) const +{ + return Triangular(derived()); +} + +/** \returns an expression of a lower triangular matrix with a unit diagonal + * + * \sa lower(), upperWithUnitDiag() + */ +template +const Triangular MatrixBase::lowerWithUnitDiag(void) const +{ + return Triangular(derived()); +} + +/** \returns an expression of a strictly lower triangular matrix (diagonal==zero) + * FIXME could also be called strictlyLower() or lowerStrict() + * + * \sa lower(), upperWithNullDiag() + */ +template +const Triangular MatrixBase::lowerWithNullDiag(void) const +{ + return Triangular(derived()); +} + +/** \returns true if *this is approximately equal to an upper triangular matrix, + * within the precision given by \a prec. + * + * \sa isLower(), upper() + */ +template +bool MatrixBase::isUpper(RealScalar prec) const +{ + if(cols() != rows()) return false; + RealScalar maxAbsOnUpperPart = static_cast(-1); + for(int j = 0; j < cols(); j++) + for(int i = 0; i <= j; i++) + { + RealScalar absValue = ei_abs(coeff(i,j)); + if(absValue > maxAbsOnUpperPart) maxAbsOnUpperPart = absValue; + } + 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 +bool MatrixBase::isLower(RealScalar prec) const +{ + if(cols() != rows()) return false; + 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 > 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: template class DiagonalMatrix; template class DiagonalCoeffs; template class Map; -template class Eval; +// template class Eval; template class PartialRedux; +template class Triangular; template struct ei_scalar_sum_op; template struct ei_scalar_difference_op; @@ -71,5 +72,6 @@ template struct ei_scalar_min_op; template struct ei_scalar_max_op; template class Inverse; +template 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 ei_traits::Scalar MatrixBase::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()); + } + 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 class Inverse : ei_no_assignm protected: bool m_exists; - typename MatrixType::Eval m_inverse; + MatrixType m_inverse; }; template 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 +// +// 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_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 class QR +{ + public: + + typedef typename MatrixType::Scalar Scalar; + typedef Matrix RMatrixType; + typedef Matrix 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 +void QR::_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 QR::RMatrixType QR::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 +MatrixType QR::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 +const QR::type> +MatrixBase::qr() const +{ + return QR::type>(derived()); +} + + +#endif // EIGEN_QR_H