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