diff options
author | Gael Guennebaud <g.gael@free.fr> | 2009-07-06 23:43:20 +0200 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2009-07-06 23:43:20 +0200 |
commit | 1aea45335ff60c6a5ed6dc06cd798d050eff661a (patch) | |
tree | d4dfb36ebe76b45f25dd984ac681ac9f92894515 /Eigen | |
parent | 08e419dcb151da41f169304f751e5467cf0c7b4a (diff) |
* bybye Part, welcome TriangularView and SelfAdjointView.
* move solveTriangular*() to TriangularView::solve*()
* move .llt() to SelfAdjointView
* add a high level wrapper to the efficient selfadjoint * vector product
* improve LLT so that we can specify which triangular part is meaningless
=> there are still many things to do (doc, cleaning, improve the matrix products, etc.)
Diffstat (limited to 'Eigen')
-rw-r--r-- | Eigen/Core | 5 | ||||
-rw-r--r-- | Eigen/src/Cholesky/LDLT.h | 7 | ||||
-rw-r--r-- | Eigen/src/Cholesky/LLT.h | 197 | ||||
-rw-r--r-- | Eigen/src/Core/Matrix.h | 15 | ||||
-rw-r--r-- | Eigen/src/Core/MatrixBase.h | 16 | ||||
-rw-r--r-- | Eigen/src/Core/Part.h | 375 | ||||
-rw-r--r-- | Eigen/src/Core/Product.h | 2 | ||||
-rw-r--r-- | Eigen/src/Core/SelfAdjointView.h | 171 | ||||
-rw-r--r-- | Eigen/src/Core/SolveTriangular.h | 94 | ||||
-rw-r--r-- | Eigen/src/Core/Transpose.h | 4 | ||||
-rw-r--r-- | Eigen/src/Core/TriangularMatrix.h | 637 | ||||
-rw-r--r-- | Eigen/src/Core/util/ForwardDeclarations.h | 7 | ||||
-rw-r--r-- | Eigen/src/QR/QR.h | 2 |
13 files changed, 1057 insertions, 475 deletions
diff --git a/Eigen/Core b/Eigen/Core index 792fc0135..89b18d201 100644 --- a/Eigen/Core +++ b/Eigen/Core @@ -159,7 +159,6 @@ namespace Eigen { #include "src/Core/CwiseNullaryOp.h" #include "src/Core/CwiseUnaryView.h" #include "src/Core/Dot.h" -#include "src/Core/SolveTriangular.h" #include "src/Core/MapBase.h" #include "src/Core/Map.h" #include "src/Core/Block.h" @@ -174,11 +173,13 @@ namespace Eigen { #include "src/Core/IO.h" #include "src/Core/Swap.h" #include "src/Core/CommaInitializer.h" -#include "src/Core/Part.h" #include "src/Core/Product.h" #include "src/Core/products/GeneralMatrixMatrix.h" #include "src/Core/products/GeneralMatrixVector.h" #include "src/Core/products/SelfadjointMatrixVector.h" +#include "src/Core/TriangularMatrix.h" +#include "src/Core/SelfAdjointView.h" +#include "src/Core/SolveTriangular.h" } // namespace Eigen diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index 94660245e..18029bab8 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -80,7 +80,7 @@ template<typename MatrixType> class LDLT } /** \returns the lower triangular matrix L */ - inline Part<MatrixType, UnitLowerTriangular> matrixL(void) const + inline TriangularView<MatrixType, UnitLowerTriangular> matrixL(void) const { ei_assert(m_isInitialized && "LDLT is not initialized."); return m_matrix; @@ -282,13 +282,14 @@ bool LDLT<MatrixType>::solveInPlace(MatrixBase<Derived> &bAndX) const for(int i = 0; i < size; ++i) bAndX.row(m_transpositions.coeff(i)).swap(bAndX.row(i)); // y = L^-1 z - matrixL().solveTriangularInPlace(bAndX); + //matrixL().solveInPlace(bAndX); + m_matrix.template triangularView<UnitLowerTriangular>().solveInPlace(bAndX); // w = D^-1 y bAndX = (m_matrix.diagonal().cwise().inverse().asDiagonal() * bAndX).lazy(); // u = L^-T w - m_matrix.adjoint().template part<UnitUpperTriangular>().solveTriangularInPlace(bAndX); + m_matrix.adjoint().template triangularView<UnitUpperTriangular>().solveInPlace(bAndX); // x = P^T u for (int i = size-1; i >= 0; --i) bAndX.row(m_transpositions.coeff(i)).swap(bAndX.row(i)); diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h index 215cda346..bad360579 100644 --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h @@ -25,6 +25,8 @@ #ifndef EIGEN_LLT_H #define EIGEN_LLT_H +template<typename MatrixType, int UpLo> struct LLT_Traits; + /** \ingroup cholesky_Module * * \class LLT @@ -51,7 +53,7 @@ * Note that during the decomposition, only the upper triangular part of A is considered. Therefore, * the strict lower part does not have to store correct values. */ -template<typename MatrixType> class LLT +template<typename MatrixType, int _UpLo> class LLT { private: typedef typename MatrixType::Scalar Scalar; @@ -60,9 +62,12 @@ template<typename MatrixType> class LLT enum { PacketSize = ei_packet_traits<Scalar>::size, - AlignmentMask = int(PacketSize)-1 + AlignmentMask = int(PacketSize)-1, + UpLo = _UpLo }; + typedef LLT_Traits<MatrixType,UpLo> Traits; + public: /** @@ -80,11 +85,18 @@ template<typename MatrixType> class LLT compute(matrix); } - /** \returns the lower triangular matrix L */ - inline Part<MatrixType, LowerTriangular> matrixL(void) const + /** \returns a view of the upper triangular matrix U */ + inline typename Traits::MatrixU matrixU() const + { + ei_assert(m_isInitialized && "LLT is not initialized."); + return Traits::getU(m_matrix); + } + + /** \returns a view of the lower triangular matrix L */ + inline typename Traits::MatrixL matrixL() const { ei_assert(m_isInitialized && "LLT is not initialized."); - return m_matrix; + return Traits::getL(m_matrix); } template<typename RhsDerived, typename ResultType> @@ -104,51 +116,162 @@ template<typename MatrixType> class LLT bool m_isInitialized; }; -/** Computes / recomputes the Cholesky decomposition A = LL^* = U^*U of \a matrix - */ -template<typename MatrixType> -void LLT<MatrixType>::compute(const MatrixType& a) +template<typename MatrixType/*, int UpLo*/> +bool ei_inplace_llt_lo(MatrixType& mat) { - assert(a.rows()==a.cols()); - const int size = a.rows(); - m_matrix.resize(size, size); + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + assert(mat.rows()==mat.cols()); + const int size = mat.rows(); + Matrix<Scalar,Dynamic,1> aux(size); // The biggest overall is the point of reference to which further diagonals // are compared; if any diagonal is negligible compared // to the largest overall, the algorithm bails. This cutoff is suggested // in "Analysis of the Cholesky Decomposition of a Semi-definite Matrix" by // Nicholas J. Higham. Also see "Accuracy and Stability of Numerical // Algorithms" page 217, also by Higham. - const RealScalar cutoff = machine_epsilon<Scalar>() * size * a.diagonal().cwise().abs().maxCoeff(); + const RealScalar cutoff = machine_epsilon<Scalar>() * size * mat.diagonal().cwise().abs().maxCoeff(); RealScalar x; - x = ei_real(a.coeff(0,0)); - m_matrix.coeffRef(0,0) = ei_sqrt(x); + x = ei_real(mat.coeff(0,0)); + mat.coeffRef(0,0) = ei_sqrt(x); if(size==1) { - m_isInitialized = true; - return; + return true; } - m_matrix.col(0).end(size-1) = a.row(0).end(size-1).adjoint() / ei_real(m_matrix.coeff(0,0)); + mat.col(0).end(size-1) = mat.col(0).end(size-1) / ei_real(mat.coeff(0,0)); for (int j = 1; j < size; ++j) { - x = ei_real(a.coeff(j,j)) - m_matrix.row(j).start(j).squaredNorm(); + x = ei_real(mat.coeff(j,j)) - mat.row(j).start(j).squaredNorm(); if (ei_abs(x) < cutoff) continue; - m_matrix.coeffRef(j,j) = x = ei_sqrt(x); + mat.coeffRef(j,j) = x = ei_sqrt(x); int endSize = size-j-1; if (endSize>0) { // Note that when all matrix columns have good alignment, then the following // product is guaranteed to be optimal with respect to alignment. - m_matrix.col(j).end(endSize) = - (m_matrix.block(j+1, 0, endSize, j) * m_matrix.row(j).start(j).adjoint()).lazy(); + aux.end(endSize) = + (mat.block(j+1, 0, endSize, j) * mat.row(j).start(j).adjoint()).lazy(); + + mat.col(j).end(endSize) = (mat.col(j).end(endSize) - aux.end(endSize) ) / x; + // TODO improve the products so that the following is efficient: +// mat.col(j).end(endSize) -= (mat.block(j+1, 0, endSize, j) * mat.row(j).start(j).adjoint()); +// mat.col(j).end(endSize) *= Scalar(1)/x; + } + } + + return true; +} + +template<typename MatrixType/*, int UpLo*/> +bool ei_inplace_llt_up(MatrixType& mat) +{ + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + assert(mat.rows()==mat.cols()); + const int size = mat.rows(); + Matrix<Scalar,Dynamic,1> aux(size); + + const RealScalar cutoff = machine_epsilon<Scalar>() * size * mat.diagonal().cwise().abs().maxCoeff(); + RealScalar x; + x = ei_real(mat.coeff(0,0)); + mat.coeffRef(0,0) = ei_sqrt(x); + if(size==1) + { + return true; + } + mat.row(0).end(size-1) = mat.row(0).end(size-1) / ei_real(mat.coeff(0,0)); + for (int j = 1; j < size; ++j) + { + x = ei_real(mat.coeff(j,j)) - mat.col(j).start(j).squaredNorm(); + if (ei_abs(x) < cutoff) continue; + + mat.coeffRef(j,j) = x = ei_sqrt(x); - // FIXME could use a.col instead of a.row - m_matrix.col(j).end(endSize) = (a.row(j).end(endSize).adjoint() - - m_matrix.col(j).end(endSize) ) / x; + int endSize = size-j-1; + if (endSize>0) { + aux.start(endSize) = + (mat.block(0, j+1, j, endSize).adjoint() * mat.col(j).start(j)).lazy(); + + mat.row(j).end(endSize) = (mat.row(j).end(endSize) - aux.start(endSize).adjoint() ) / x; } } - m_isInitialized = true; + return true; +} + +template<typename MatrixType> struct LLT_Traits<MatrixType,LowerTriangular> +{ + typedef TriangularView<MatrixType, LowerTriangular> MatrixL; + typedef TriangularView<NestByValue<typename MatrixType::AdjointReturnType>, UpperTriangular> MatrixU; + inline static MatrixL getL(const MatrixType& m) { return m; } + inline static MatrixU getU(const MatrixType& m) { return m.adjoint().nestByValue(); } + static bool inplace_decomposition(MatrixType& m) + { return ei_inplace_llt_lo(m); } +}; + +template<typename MatrixType> struct LLT_Traits<MatrixType,UpperTriangular> +{ + typedef TriangularView<NestByValue<typename MatrixType::AdjointReturnType>, LowerTriangular> MatrixL; + typedef TriangularView<MatrixType, UpperTriangular> MatrixU; + inline static MatrixL getL(const MatrixType& m) { return m.adjoint().nestByValue(); } + inline static MatrixU getU(const MatrixType& m) { return m; } + static bool inplace_decomposition(MatrixType& m) + { return ei_inplace_llt_up(m); } +}; + +/** Computes / recomputes the Cholesky decomposition A = LL^* = U^*U of \a matrix + */ +template<typename MatrixType, int _UpLo> +void LLT<MatrixType,_UpLo>::compute(const MatrixType& a) +{ +// assert(a.rows()==a.cols()); +// const int size = a.rows(); +// m_matrix.resize(size, size); +// // The biggest overall is the point of reference to which further diagonals +// // are compared; if any diagonal is negligible compared +// // to the largest overall, the algorithm bails. This cutoff is suggested +// // in "Analysis of the Cholesky Decomposition of a Semi-definite Matrix" by +// // Nicholas J. Higham. Also see "Accuracy and Stability of Numerical +// // Algorithms" page 217, also by Higham. +// const RealScalar cutoff = machine_epsilon<Scalar>() * size * a.diagonal().cwise().abs().maxCoeff(); +// RealScalar x; +// x = ei_real(a.coeff(0,0)); +// m_matrix.coeffRef(0,0) = ei_sqrt(x); +// if(size==1) +// { +// m_isInitialized = true; +// return; +// } +// m_matrix.col(0).end(size-1) = a.row(0).end(size-1).adjoint() / ei_real(m_matrix.coeff(0,0)); +// for (int j = 1; j < size; ++j) +// { +// x = ei_real(a.coeff(j,j)) - m_matrix.row(j).start(j).squaredNorm(); +// if (ei_abs(x) < cutoff) continue; +// +// m_matrix.coeffRef(j,j) = x = ei_sqrt(x); +// +// int endSize = size-j-1; +// if (endSize>0) { +// // Note that when all matrix columns have good alignment, then the following +// // product is guaranteed to be optimal with respect to alignment. +// m_matrix.col(j).end(endSize) = +// (m_matrix.block(j+1, 0, endSize, j) * m_matrix.row(j).start(j).adjoint()).lazy(); +// +// // FIXME could use a.col instead of a.row +// m_matrix.col(j).end(endSize) = (a.row(j).end(endSize).adjoint() +// - m_matrix.col(j).end(endSize) ) / x; +// } +// } +// +// m_isInitialized = true; + + assert(a.rows()==a.cols()); + const int size = a.rows(); + m_matrix.resize(size, size); + m_matrix = a; + + m_isInitialized = Traits::inplace_decomposition(m_matrix); } /** Computes the solution x of \f$ A x = b \f$ using the current decomposition of A. @@ -164,9 +287,9 @@ void LLT<MatrixType>::compute(const MatrixType& a) * * \sa LLT::solveInPlace(), MatrixBase::llt() */ -template<typename MatrixType> +template<typename MatrixType, int _UpLo> template<typename RhsDerived, typename ResultType> -bool LLT<MatrixType>::solve(const MatrixBase<RhsDerived> &b, ResultType *result) const +bool LLT<MatrixType,_UpLo>::solve(const MatrixBase<RhsDerived> &b, ResultType *result) const { ei_assert(m_isInitialized && "LLT is not initialized."); const int size = m_matrix.rows(); @@ -185,15 +308,15 @@ bool LLT<MatrixType>::solve(const MatrixBase<RhsDerived> &b, ResultType *result) * * \sa LLT::solve(), MatrixBase::llt() */ -template<typename MatrixType> +template<typename MatrixType, int _UpLo> template<typename Derived> -bool LLT<MatrixType>::solveInPlace(MatrixBase<Derived> &bAndX) const +bool LLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const { ei_assert(m_isInitialized && "LLT is not initialized."); const int size = m_matrix.rows(); ei_assert(size==bAndX.rows()); - matrixL().solveTriangularInPlace(bAndX); - m_matrix.adjoint().template part<UpperTriangular>().solveTriangularInPlace(bAndX); + matrixL().solveInPlace(bAndX); + matrixU().solveInPlace(bAndX); return true; } @@ -207,4 +330,14 @@ MatrixBase<Derived>::llt() const return LLT<PlainMatrixType>(derived()); } +/** \cholesky_module + * \returns the LLT decomposition of \c *this + */ +template<typename MatrixType, unsigned int UpLo> +inline const LLT<typename SelfAdjointView<MatrixType, UpLo>::PlainMatrixType, UpLo> +SelfAdjointView<MatrixType, UpLo>::llt() const +{ + return LLT<PlainMatrixType,UpLo>(m_matrix); +} + #endif // EIGEN_LLT_H diff --git a/Eigen/src/Core/Matrix.h b/Eigen/src/Core/Matrix.h index 8cf5bc64b..b7a8ebf21 100644 --- a/Eigen/src/Core/Matrix.h +++ b/Eigen/src/Core/Matrix.h @@ -456,6 +456,21 @@ class Matrix *this = other; } + template<typename TriangularDerived> + EIGEN_STRONG_INLINE Matrix& operator=(const TriangularBase<TriangularDerived> &other) + { + resize(other.rows(), other.cols()); + Base::operator=(other.derived()); + return *this; + } + + template<typename TriangularDerived> + EIGEN_STRONG_INLINE Matrix(const TriangularBase<TriangularDerived> &other) + : m_storage(other.rows() * other.cols(), other.rows(), other.cols()) + { + *this = other; + } + /** Override MatrixBase::swap() since for dynamic-sized matrices of same type it is enough to swap the * data pointers. */ diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index b0ff09531..82417a144 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -297,6 +297,10 @@ template<typename Derived> class MatrixBase template<typename OtherDerived> Derived& lazyAssign(const Flagged<OtherDerived, 0, EvalBeforeNestingBit | EvalBeforeAssigningBit>& other) { return lazyAssign(other._expression()); } + + /** Overloaded for fast triangular part to dense matrix evaluation */ + template<typename TriangularDerived> + Derived& lazyAssign(const TriangularBase<TriangularDerived> &other); #endif // not EIGEN_PARSED_BY_DOXYGEN CommaInitializer<Derived> operator<< (const Scalar& s); @@ -402,6 +406,9 @@ template<typename Derived> class MatrixBase template<typename DiagonalDerived> Derived& operator=(const DiagonalBase<DiagonalDerived> &other); + template<typename TriangularDerived> + Derived& operator=(const TriangularBase<TriangularDerived> &other); + template<typename OtherDerived> typename ei_plain_matrix_type_column_major<OtherDerived>::type solveTriangular(const MatrixBase<OtherDerived>& other) const; @@ -477,9 +484,14 @@ template<typename Derived> class MatrixBase Diagonal<Derived, Dynamic> diagonal(int index); const Diagonal<Derived, Dynamic> diagonal(int index) const; - template<unsigned int Mode> Part<Derived, Mode> part(); - template<unsigned int Mode> const Part<Derived, Mode> part() const; + template<unsigned int Mode> TriangularView<Derived, Mode> part(); + template<unsigned int Mode> const TriangularView<Derived, Mode> part() const; + + template<unsigned int Mode> TriangularView<Derived, Mode> triangularView(); + template<unsigned int Mode> const TriangularView<Derived, Mode> triangularView() const; + template<unsigned int UpLo> SelfAdjointView<Derived, UpLo> selfadjointView(); + template<unsigned int UpLo> const SelfAdjointView<Derived, UpLo> selfadjointView() const; static const ConstantReturnType Constant(int rows, int cols, const Scalar& value); diff --git a/Eigen/src/Core/Part.h b/Eigen/src/Core/Part.h deleted file mode 100644 index a831c2a42..000000000 --- a/Eigen/src/Core/Part.h +++ /dev/null @@ -1,375 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2008 Benoit Jacob <jacob.benoit.1@gmail.com> -// 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_PART_H -#define EIGEN_PART_H - -/** \nonstableyet - * \class Part - * - * \brief Expression of a triangular matrix extracted from a given matrix - * - * \param MatrixType the type of the object in which we are taking the triangular part - * \param Mode the kind of triangular matrix expression to construct. Can be UpperTriangular, StrictlyUpperTriangular, - * UnitUpperTriangular, LowerTriangular, StrictlyLowerTriangular, UnitLowerTriangular. This is in fact a bit field; it must have either - * UpperTriangularBit or LowerTriangularBit, and additionnaly it may have either ZeroDiagBit or - * UnitDiagBit. - * - * This class represents an expression of the upper or lower triangular part of - * a square matrix, possibly with a further assumption on the diagonal. It is the return type - * of MatrixBase::part() and most of the time this is the only way it is used. - * - * \sa MatrixBase::part() - */ -template<typename MatrixType, unsigned int Mode> -struct ei_traits<Part<MatrixType, Mode> > : ei_traits<MatrixType> -{ - typedef typename ei_nested<MatrixType>::type MatrixTypeNested; - typedef typename ei_unref<MatrixTypeNested>::type _MatrixTypeNested; - enum { - Flags = (_MatrixTypeNested::Flags & (HereditaryBits) & (~(PacketAccessBit | DirectAccessBit | LinearAccessBit))) | Mode, - CoeffReadCost = _MatrixTypeNested::CoeffReadCost - }; -}; - -template<typename MatrixType, unsigned int Mode> class Part - : public MatrixBase<Part<MatrixType, Mode> > -{ - public: - - EIGEN_GENERIC_PUBLIC_INTERFACE(Part) - - inline Part(const MatrixType& matrix) : m_matrix(matrix) - { ei_assert(ei_are_flags_consistent<Mode>::ret); } - - /** \sa MatrixBase::operator+=() */ - template<typename Other> Part& operator+=(const Other& other); - /** \sa MatrixBase::operator-=() */ - template<typename Other> Part& operator-=(const Other& other); - /** \sa MatrixBase::operator*=() */ - Part& operator*=(const typename ei_traits<MatrixType>::Scalar& other); - /** \sa MatrixBase::operator/=() */ - Part& operator/=(const typename ei_traits<MatrixType>::Scalar& other); - - /** \sa operator=(), MatrixBase::lazyAssign() */ - template<typename Other> void lazyAssign(const Other& other); - /** \sa MatrixBase::operator=() */ - template<typename Other> Part& operator=(const Other& other); - - inline int rows() const { return m_matrix.rows(); } - inline int cols() const { return m_matrix.cols(); } - inline int stride() const { return m_matrix.stride(); } - - inline Scalar coeff(int row, int col) const - { - // SelfAdjointBit doesn't play any role here: just because a matrix is selfadjoint doesn't say anything about - // each individual coefficient, except for the not-very-useful-here fact that diagonal coefficients are real. - if( ((Flags & LowerTriangularBit) && (col>row)) || ((Flags & UpperTriangularBit) && (row>col)) ) - return (Scalar)0; - if(Flags & UnitDiagBit) - return col==row ? (Scalar)1 : m_matrix.coeff(row, col); - else if(Flags & ZeroDiagBit) - return col==row ? (Scalar)0 : m_matrix.coeff(row, col); - else - return m_matrix.coeff(row, col); - } - - inline Scalar& coeffRef(int row, int col) - { - EIGEN_STATIC_ASSERT(!(Flags & UnitDiagBit), WRITING_TO_TRIANGULAR_PART_WITH_UNIT_DIAGONAL_IS_NOT_SUPPORTED) - EIGEN_STATIC_ASSERT(!(Flags & SelfAdjointBit), COEFFICIENT_WRITE_ACCESS_TO_SELFADJOINT_NOT_SUPPORTED) - ei_assert( (Mode==UpperTriangular && col>=row) - || (Mode==LowerTriangular && col<=row) - || (Mode==StrictlyUpperTriangular && col>row) - || (Mode==StrictlyLowerTriangular && col<row)); - return m_matrix.const_cast_derived().coeffRef(row, col); - } - - /** \internal */ - const MatrixType& _expression() const { return m_matrix; } - - /** discard any writes to a row */ - const Block<Part, 1, ColsAtCompileTime> row(int i) { return Base::row(i); } - const Block<Part, 1, ColsAtCompileTime> row(int i) const { return Base::row(i); } - /** discard any writes to a column */ - const Block<Part, RowsAtCompileTime, 1> col(int i) { return Base::col(i); } - const Block<Part, RowsAtCompileTime, 1> col(int i) const { return Base::col(i); } - - template<typename OtherDerived> - void swap(const MatrixBase<OtherDerived>& other) - { - Part<SwapWrapper<MatrixType>,Mode>(const_cast<MatrixType&>(m_matrix)).lazyAssign(other.derived()); - } - - protected: - - const typename MatrixType::Nested m_matrix; -}; - -/** \nonstableyet - * \returns an expression of a triangular matrix extracted from the current matrix - * - * The parameter \a Mode can have the following values: \c UpperTriangular, \c StrictlyUpperTriangular, \c UnitUpperTriangular, - * \c LowerTriangular, \c StrictlyLowerTriangular, \c UnitLowerTriangular. - * - * \addexample PartExample \label How to extract a triangular part of an arbitrary matrix - * - * Example: \include MatrixBase_extract.cpp - * Output: \verbinclude MatrixBase_extract.out - * - * \sa class Part, part(), marked() - */ -template<typename Derived> -template<unsigned int Mode> -const Part<Derived, Mode> MatrixBase<Derived>::part() const -{ - return derived(); -} - -template<typename MatrixType, unsigned int Mode> -template<typename Other> -inline Part<MatrixType, Mode>& Part<MatrixType, Mode>::operator=(const Other& other) -{ - if(Other::Flags & EvalBeforeAssigningBit) - { - typename MatrixBase<Other>::PlainMatrixType other_evaluated(other.rows(), other.cols()); - other_evaluated.template part<Mode>().lazyAssign(other); - lazyAssign(other_evaluated); - } - else - lazyAssign(other.derived()); - return *this; -} - -template<typename Derived1, typename Derived2, unsigned int Mode, int UnrollCount> -struct ei_part_assignment_impl -{ - enum { - col = (UnrollCount-1) / Derived1::RowsAtCompileTime, - row = (UnrollCount-1) % Derived1::RowsAtCompileTime - }; - - inline static void run(Derived1 &dst, const Derived2 &src) - { - ei_part_assignment_impl<Derived1, Derived2, Mode, UnrollCount-1>::run(dst, src); - - if(Mode == SelfAdjoint) - { - if(row == col) - dst.coeffRef(row, col) = ei_real(src.coeff(row, col)); - else if(row < col) - dst.coeffRef(col, row) = ei_conj(dst.coeffRef(row, col) = src.coeff(row, col)); - } - else - { - ei_assert(Mode == UpperTriangular || Mode == LowerTriangular || Mode == StrictlyUpperTriangular || Mode == StrictlyLowerTriangular); - if((Mode == UpperTriangular && row <= col) - || (Mode == LowerTriangular && row >= col) - || (Mode == StrictlyUpperTriangular && row < col) - || (Mode == StrictlyLowerTriangular && row > col)) - dst.copyCoeff(row, col, src); - } - } -}; - -template<typename Derived1, typename Derived2, unsigned int Mode> -struct ei_part_assignment_impl<Derived1, Derived2, Mode, 1> -{ - inline static void run(Derived1 &dst, const Derived2 &src) - { - if(!(Mode & ZeroDiagBit)) - dst.copyCoeff(0, 0, src); - } -}; - -// prevent buggy user code from causing an infinite recursion -template<typename Derived1, typename Derived2, unsigned int Mode> -struct ei_part_assignment_impl<Derived1, Derived2, Mode, 0> -{ - inline static void run(Derived1 &, const Derived2 &) {} -}; - -template<typename Derived1, typename Derived2> -struct ei_part_assignment_impl<Derived1, Derived2, UpperTriangular, Dynamic> -{ - inline static void run(Derived1 &dst, const Derived2 &src) - { - for(int j = 0; j < dst.cols(); ++j) - for(int i = 0; i <= j; ++i) - dst.copyCoeff(i, j, src); - } -}; - -template<typename Derived1, typename Derived2> -struct ei_part_assignment_impl<Derived1, Derived2, LowerTriangular, Dynamic> -{ - inline static void run(Derived1 &dst, const Derived2 &src) - { - for(int j = 0; j < dst.cols(); ++j) - for(int i = j; i < dst.rows(); ++i) - dst.copyCoeff(i, j, src); - } -}; - -template<typename Derived1, typename Derived2> -struct ei_part_assignment_impl<Derived1, Derived2, StrictlyUpperTriangular, Dynamic> -{ - inline static void run(Derived1 &dst, const Derived2 &src) - { - for(int j = 0; j < dst.cols(); ++j) - for(int i = 0; i < j; ++i) - dst.copyCoeff(i, j, src); - } -}; -template<typename Derived1, typename Derived2> -struct ei_part_assignment_impl<Derived1, Derived2, StrictlyLowerTriangular, Dynamic> -{ - inline static void run(Derived1 &dst, const Derived2 &src) - { - for(int j = 0; j < dst.cols(); ++j) - for(int i = j+1; i < dst.rows(); ++i) - dst.copyCoeff(i, j, src); - } -}; -template<typename Derived1, typename Derived2> -struct ei_part_assignment_impl<Derived1, Derived2, SelfAdjoint, Dynamic> -{ - inline static void run(Derived1 &dst, const Derived2 &src) - { - for(int j = 0; j < dst.cols(); ++j) - { - for(int i = 0; i < j; ++i) - dst.coeffRef(j, i) = ei_conj(dst.coeffRef(i, j) = src.coeff(i, j)); - dst.coeffRef(j, j) = ei_real(src.coeff(j, j)); - } - } -}; - -template<typename MatrixType, unsigned int Mode> -template<typename Other> -void Part<MatrixType, Mode>::lazyAssign(const Other& other) -{ - const bool unroll = MatrixType::SizeAtCompileTime * Other::CoeffReadCost / 2 <= EIGEN_UNROLLING_LIMIT; - ei_assert(m_matrix.rows() == other.rows() && m_matrix.cols() == other.cols()); - - ei_part_assignment_impl - <MatrixType, Other, Mode, - unroll ? int(MatrixType::SizeAtCompileTime) : Dynamic - >::run(m_matrix.const_cast_derived(), other.derived()); -} - -/** \nonstableyet - * \returns a lvalue pseudo-expression allowing to perform special operations on \c *this. - * - * The \a Mode parameter can have the following values: \c UpperTriangular, \c StrictlyUpperTriangular, \c LowerTriangular, - * \c StrictlyLowerTriangular, \c SelfAdjoint. - * - * \addexample PartExample \label How to write to a triangular part of a matrix - * - * Example: \include MatrixBase_part.cpp - * Output: \verbinclude MatrixBase_part.out - * - * \sa class Part, MatrixBase::extract(), MatrixBase::marked() - */ -template<typename Derived> -template<unsigned int Mode> -inline Part<Derived, Mode> MatrixBase<Derived>::part() -{ - return Part<Derived, Mode>(derived()); -} - -/** \returns true if *this is approximately equal to an upper triangular matrix, - * within the precision given by \a prec. - * - * \sa isLowerTriangular(), extract(), part(), marked() - */ -template<typename Derived> -bool MatrixBase<Derived>::isUpperTriangular(RealScalar prec) const -{ - if(cols() != rows()) return false; - RealScalar maxAbsOnUpperTriangularPart = 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 > maxAbsOnUpperTriangularPart) maxAbsOnUpperTriangularPart = absValue; - } - for(int j = 0; j < cols()-1; ++j) - for(int i = j+1; i < rows(); ++i) - if(!ei_isMuchSmallerThan(coeff(i, j), maxAbsOnUpperTriangularPart, 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 isUpperTriangular(), extract(), part(), marked() - */ -template<typename Derived> -bool MatrixBase<Derived>::isLowerTriangular(RealScalar prec) const -{ - if(cols() != rows()) return false; - RealScalar maxAbsOnLowerTriangularPart = 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 > maxAbsOnLowerTriangularPart) maxAbsOnLowerTriangularPart = absValue; - } - for(int j = 1; j < cols(); ++j) - for(int i = 0; i < j; ++i) - if(!ei_isMuchSmallerThan(coeff(i, j), maxAbsOnLowerTriangularPart, prec)) return false; - return true; -} - -template<typename MatrixType, unsigned int Mode> -template<typename Other> -inline Part<MatrixType, Mode>& Part<MatrixType, Mode>::operator+=(const Other& other) -{ - return *this = m_matrix + other; -} - -template<typename MatrixType, unsigned int Mode> -template<typename Other> -inline Part<MatrixType, Mode>& Part<MatrixType, Mode>::operator-=(const Other& other) -{ - return *this = m_matrix - other; -} - -template<typename MatrixType, unsigned int Mode> -inline Part<MatrixType, Mode>& Part<MatrixType, Mode>::operator*= -(const typename ei_traits<MatrixType>::Scalar& other) -{ - return *this = m_matrix * other; -} - -template<typename MatrixType, unsigned int Mode> -inline Part<MatrixType, Mode>& Part<MatrixType, Mode>::operator/= -(const typename ei_traits<MatrixType>::Scalar& other) -{ - return *this = m_matrix / other; -} - -#endif // EIGEN_PART_H diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h index fe6c5027d..75cfadfa1 100644 --- a/Eigen/src/Core/Product.h +++ b/Eigen/src/Core/Product.h @@ -708,7 +708,7 @@ inline Derived& MatrixBase<Derived>::lazyAssign(const Product<Lhs,Rhs,CacheFrien } else { - lazyAssign<Product<Lhs,Rhs,CacheFriendlyProduct> >(product); + lazyAssign(static_cast<const MatrixBase<Product<Lhs,Rhs,CacheFriendlyProduct> > &>(product)); } return derived(); } diff --git a/Eigen/src/Core/SelfAdjointView.h b/Eigen/src/Core/SelfAdjointView.h new file mode 100644 index 000000000..2f66cfa45 --- /dev/null +++ b/Eigen/src/Core/SelfAdjointView.h @@ -0,0 +1,171 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 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_SELFADJOINTMATRIX_H +#define EIGEN_SELFADJOINTMATRIX_H + +/** \class SelfAdjointView + * \nonstableyet + * + * \brief Expression of a selfadjoint matrix from a triangular part of a dense matrix + * + * \param MatrixType the type of the dense matrix storing the coefficients + * \param TriangularPart can be either \c LowerTriangular or \c UpperTriangular + * + * This class is an expression of a sefladjoint matrix from a triangular part of a matrix + * with given dense storage of the coefficients. It is the return type of MatrixBase::selfadjointView() + * and most of the time this is the only way that it is used. + * + * \sa class TriangularBase, MatrixBase::selfAdjointView() + */ +template<typename MatrixType, unsigned int TriangularPart> +struct ei_traits<SelfAdjointView<MatrixType, TriangularPart> > : ei_traits<MatrixType> +{ + typedef typename ei_nested<MatrixType>::type MatrixTypeNested; + typedef typename ei_unref<MatrixTypeNested>::type _MatrixTypeNested; + typedef MatrixType ExpressionType; + enum { + Mode = TriangularPart | SelfAdjointBit, + Flags = _MatrixTypeNested::Flags & (HereditaryBits) + & (~(PacketAccessBit | DirectAccessBit | LinearAccessBit)), // FIXME these flags should be preserved + CoeffReadCost = _MatrixTypeNested::CoeffReadCost + }; +}; + +template<typename Lhs,typename Rhs> +struct ei_selfadjoint_vector_product_returntype; + +// FIXME could also be called SelfAdjointWrapper to be consistent with DiagonalWrapper ?? +template<typename MatrixType, unsigned int UpLo> class SelfAdjointView + : public TriangularBase<SelfAdjointView<MatrixType, UpLo> > +{ + public: + + typedef TriangularBase<SelfAdjointView> Base; + typedef typename ei_traits<SelfAdjointView>::Scalar Scalar; + enum { + Mode = ei_traits<SelfAdjointView>::Mode + }; + typedef typename MatrixType::PlainMatrixType PlainMatrixType; + + inline SelfAdjointView(const MatrixType& matrix) : m_matrix(matrix) + { ei_assert(ei_are_flags_consistent<Mode>::ret); } + + inline int rows() const { return m_matrix.rows(); } + inline int cols() const { return m_matrix.cols(); } + inline int stride() const { return m_matrix.stride(); } + + /** \sa MatrixBase::coeff() + * \warning the coordinates must fit into the referenced triangular part + */ + inline Scalar coeff(int row, int col) const + { + Base::check_coordinates_internal(row, col); + return m_matrix.coeff(row, col); + } + + /** \sa MatrixBase::coeffRef() + * \warning the coordinates must fit into the referenced triangular part + */ + inline Scalar& coeffRef(int row, int col) + { + Base::check_coordinates_internal(row, col); + return m_matrix.const_cast_derived().coeffRef(row, col); + } + + /** \internal */ + const MatrixType& _expression() const { return m_matrix; } + + /** Efficient self-adjoint matrix times vector product */ + // TODO this product is far to be ready + template<typename OtherDerived> + ei_selfadjoint_vector_product_returntype<SelfAdjointView,OtherDerived> + operator*(const MatrixBase<OtherDerived>& rhs) const + { + return ei_selfadjoint_vector_product_returntype<SelfAdjointView,OtherDerived>(*this, rhs.derived()); + } + +/////////// Cholesky module /////////// + + const LLT<PlainMatrixType, UpLo> llt() const; + const LDLT<PlainMatrixType> ldlt() const; + + protected: + + const typename MatrixType::Nested m_matrix; +}; + +/*************************************************************************** +* Wrapper to ei_product_selfadjoint_vector +***************************************************************************/ + +template<typename Lhs,typename Rhs> +struct ei_selfadjoint_vector_product_returntype + : public ReturnByValue<ei_selfadjoint_vector_product_returntype<Lhs,Rhs>, + Matrix<typename ei_traits<Rhs>::Scalar, + Rhs::RowsAtCompileTime,Rhs::ColsAtCompileTime> > +{ + typedef typename ei_cleantype<typename Rhs::Nested>::type RhsNested; + ei_selfadjoint_vector_product_returntype(const Lhs& lhs, const Rhs& rhs) + : m_lhs(lhs), m_rhs(rhs) + {} + + template<typename Dest> void evalTo(Dest& dst) const + { + dst.resize(m_rhs.rows(), m_rhs.cols()); + ei_product_selfadjoint_vector<typename Lhs::Scalar,ei_traits<Lhs>::Flags&RowMajorBit, + Lhs::Mode&(UpperTriangularBit|LowerTriangularBit)> + ( + m_lhs.rows(), // size + m_lhs._expression().data(), // lhs + m_lhs.stride(), // lhsStride, + m_rhs.data(), // rhs + // int rhsIncr, + dst.data() // res + ); + } + + const Lhs m_lhs; + const typename Rhs::Nested m_rhs; +}; + +/*************************************************************************** +* Implementation of MatrixBase methods +***************************************************************************/ + +template<typename Derived> +template<unsigned int Mode> +const SelfAdjointView<Derived, Mode> MatrixBase<Derived>::selfadjointView() const +{ + return derived(); +} + +template<typename Derived> +template<unsigned int Mode> +SelfAdjointView<Derived, Mode> MatrixBase<Derived>::selfadjointView() +{ + return derived(); +} + +#endif // EIGEN_SELFADJOINTMATRIX_H diff --git a/Eigen/src/Core/SolveTriangular.h b/Eigen/src/Core/SolveTriangular.h index 19a1dd2c1..f717a3875 100644 --- a/Eigen/src/Core/SolveTriangular.h +++ b/Eigen/src/Core/SolveTriangular.h @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> +// Copyright (C) 2008-2009 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 @@ -25,33 +25,19 @@ #ifndef EIGEN_SOLVETRIANGULAR_H #define EIGEN_SOLVETRIANGULAR_H -template<typename XprType> struct ei_is_part { enum {value=false}; }; -template<typename XprType, unsigned int Mode> struct ei_is_part<Part<XprType,Mode> > { enum {value=true}; }; - -template<typename Lhs, typename Rhs, - int TriangularPart = (int(Lhs::Flags) & LowerTriangularBit) - ? LowerTriangular - : (int(Lhs::Flags) & UpperTriangularBit) - ? UpperTriangular - : 0xffffff, - int StorageOrder = ei_is_part<Lhs>::value ? 0xffffff // this is to solve ambiguous specializations - : int(Lhs::Flags) & (RowMajorBit|SparseBit) +template<typename Lhs, typename Rhs, int Mode, + int UpLo = (Mode & LowerTriangularBit) + ? LowerTriangular + : (Mode & UpperTriangularBit) + ? UpperTriangular + : -1, + int StorageOrder = int(Lhs::Flags) & RowMajorBit > struct ei_solve_triangular_selector; -// transform a Part xpr to a Flagged xpr -template<typename Lhs, unsigned int LhsMode, typename Rhs, int UpLo, int StorageOrder> -struct ei_solve_triangular_selector<Part<Lhs,LhsMode>,Rhs,UpLo,StorageOrder> -{ - static void run(const Part<Lhs,LhsMode>& lhs, Rhs& other) - { - ei_solve_triangular_selector<Flagged<Lhs,LhsMode,0>,Rhs>::run(lhs._expression(), other); - } -}; - // forward substitution, row-major -template<typename Lhs, typename Rhs, int UpLo> -struct ei_solve_triangular_selector<Lhs,Rhs,UpLo,RowMajor|IsDense> +template<typename Lhs, typename Rhs, int Mode, int UpLo> +struct ei_solve_triangular_selector<Lhs,Rhs,Mode,UpLo,RowMajor> { typedef typename Rhs::Scalar Scalar; static void run(const Lhs& lhs, Rhs& other) @@ -70,7 +56,7 @@ struct ei_solve_triangular_selector<Lhs,Rhs,UpLo,RowMajor|IsDense> for(int c=0 ; c<other.cols() ; ++c) { // process first rows using the non block version - if(!(Lhs::Flags & UnitDiagBit)) + if(!(Mode & UnitDiagBit)) { if (IsLowerTriangular) other.coeffRef(0,c) = other.coeff(0,c)/lhs.coeff(0, 0); @@ -82,7 +68,7 @@ struct ei_solve_triangular_selector<Lhs,Rhs,UpLo,RowMajor|IsDense> Scalar tmp = other.coeff(i,c) - (IsLowerTriangular ? ((lhs.row(i).start(i)) * other.col(c).start(i)).coeff(0,0) : ((lhs.row(i).end(size-i-1)) * other.col(c).end(size-i-1)).coeff(0,0)); - if (Lhs::Flags & UnitDiagBit) + if (Mode & UnitDiagBit) other.coeffRef(i,c) = tmp; else other.coeffRef(i,c) = tmp/lhs.coeff(i,i); @@ -107,7 +93,7 @@ struct ei_solve_triangular_selector<Lhs,Rhs,UpLo,RowMajor|IsDense> */ { Scalar tmp = other.coeff(startBlock,c)-btmp.coeff(IsLowerTriangular?0:3); - if (Lhs::Flags & UnitDiagBit) + if (Mode & UnitDiagBit) other.coeffRef(i,c) = tmp; else other.coeffRef(i,c) = tmp/lhs.coeff(i,i); @@ -122,7 +108,7 @@ struct ei_solve_triangular_selector<Lhs,Rhs,UpLo,RowMajor|IsDense> - ( lhs.row(i).segment(IsLowerTriangular ? startBlock : i+1, remainingSize) * other.col(c).segment(IsLowerTriangular ? startBlock : i+1, remainingSize)).coeff(0,0); - if (Lhs::Flags & UnitDiagBit) + if (Mode & UnitDiagBit) other.coeffRef(i,c) = tmp; else other.coeffRef(i,c) = tmp/lhs.coeff(i,i); @@ -137,8 +123,8 @@ struct ei_solve_triangular_selector<Lhs,Rhs,UpLo,RowMajor|IsDense> // - inv(LowerTriangular,UnitDiag,ColMajor) * Column vector // - inv(UpperTriangular, ColMajor) * Column vector // - inv(UpperTriangular,UnitDiag,ColMajor) * Column vector -template<typename Lhs, typename Rhs, int UpLo> -struct ei_solve_triangular_selector<Lhs,Rhs,UpLo,ColMajor|IsDense> +template<typename Lhs, typename Rhs, int Mode, int UpLo> +struct ei_solve_triangular_selector<Lhs,Rhs,Mode,UpLo,ColMajor> { typedef typename Rhs::Scalar Scalar; typedef typename ei_packet_traits<Scalar>::type Packet; @@ -168,7 +154,7 @@ struct ei_solve_triangular_selector<Lhs,Rhs,UpLo,ColMajor|IsDense> for (;IsLowerTriangular ? i<endBlock : i>endBlock; i += IsLowerTriangular ? 1 : -1) { - if(!(Lhs::Flags & UnitDiagBit)) + if(!(Mode & UnitDiagBit)) other.coeffRef(i,c) /= lhs.coeff(i,i); int remainingSize = IsLowerTriangular ? endBlock-i-1 : i-endBlock-1; if (remainingSize>0) @@ -203,7 +189,7 @@ struct ei_solve_triangular_selector<Lhs,Rhs,UpLo,ColMajor|IsDense> int i; for(i=blockyEnd; IsLowerTriangular ? i<size-1 : i>0; i += (IsLowerTriangular ? 1 : -1) ) { - if(!(Lhs::Flags & UnitDiagBit)) + if(!(Mode & UnitDiagBit)) other.coeffRef(i,c) /= lhs.coeff(i,i); /* NOTE we cannot use lhs.col(i).end(size-i-1) because Part::coeffRef gets called by .col() to @@ -214,7 +200,7 @@ struct ei_solve_triangular_selector<Lhs,Rhs,UpLo,ColMajor|IsDense> else other.col(c).start(i) -= other.coeffRef(i,c) * Block<Lhs,Dynamic,1>(lhs, 0,i, i, 1); } - if(!(Lhs::Flags & UnitDiagBit)) + if(!(Mode & UnitDiagBit)) other.coeffRef(i,c) /= lhs.coeff(i,i); } } @@ -224,31 +210,31 @@ struct ei_solve_triangular_selector<Lhs,Rhs,UpLo,ColMajor|IsDense> * * \nonstableyet * - * The parameter is only marked 'const' to make the C++ compiler accept a temporary expression here. + * \warning The parameter is only marked 'const' to make the C++ compiler accept a temporary expression here. * This function will const_cast it, so constness isn't honored here. * * See MatrixBase:solveTriangular() for the details. */ -template<typename Derived> -template<typename OtherDerived> -void MatrixBase<Derived>::solveTriangularInPlace(const MatrixBase<OtherDerived>& _other) const +template<typename MatrixType, unsigned int Mode> +template<typename RhsDerived> +void TriangularView<MatrixType,Mode>::solveInPlace(const MatrixBase<RhsDerived>& _rhs) const { - MatrixBase<OtherDerived>& other = _other.const_cast_derived(); - ei_assert(derived().cols() == derived().rows()); - ei_assert(derived().cols() == other.rows()); - ei_assert(!(Flags & ZeroDiagBit)); - ei_assert(Flags & (UpperTriangularBit|LowerTriangularBit)); + RhsDerived& rhs = _rhs.const_cast_derived(); + ei_assert(cols() == rows()); + ei_assert(cols() == rhs.rows()); + ei_assert(!(Mode & ZeroDiagBit)); + ei_assert(Mode & (UpperTriangularBit|LowerTriangularBit)); - enum { copy = ei_traits<OtherDerived>::Flags & RowMajorBit }; + enum { copy = ei_traits<RhsDerived>::Flags & RowMajorBit }; typedef typename ei_meta_if<copy, - typename ei_plain_matrix_type_column_major<OtherDerived>::type, OtherDerived&>::ret OtherCopy; - OtherCopy otherCopy(other.derived()); + typename ei_plain_matrix_type_column_major<RhsDerived>::type, RhsDerived&>::ret RhsCopy; + RhsCopy rhsCopy(rhs); - ei_solve_triangular_selector<Derived, typename ei_unref<OtherCopy>::type>::run(derived(), otherCopy); + ei_solve_triangular_selector<MatrixType, typename ei_unref<RhsCopy>::type, Mode>::run(_expression(), rhsCopy); if (copy) - other = otherCopy; + rhs = rhsCopy; } /** \returns the product of the inverse of \c *this with \a other, \a *this being triangular. @@ -282,15 +268,15 @@ void MatrixBase<Derived>::solveTriangularInPlace(const MatrixBase<OtherDerived>& * M * T^1 <=> T.transpose().solveTriangularInPlace(M.transpose()); * \endcode * - * \sa solveTriangularInPlace(), marked(), extract() + * \sa solveTriangularInPlace() */ -template<typename Derived> -template<typename OtherDerived> -typename ei_plain_matrix_type_column_major<OtherDerived>::type -MatrixBase<Derived>::solveTriangular(const MatrixBase<OtherDerived>& other) const +template<typename Derived, unsigned int Mode> +template<typename RhsDerived> +typename ei_plain_matrix_type_column_major<RhsDerived>::type +TriangularView<Derived,Mode>::solve(const MatrixBase<RhsDerived>& rhs) const { - typename ei_plain_matrix_type_column_major<OtherDerived>::type res(other); - solveTriangularInPlace(res); + typename ei_plain_matrix_type_column_major<RhsDerived>::type res(rhs); + solveInPlace(res); return res; } diff --git a/Eigen/src/Core/Transpose.h b/Eigen/src/Core/Transpose.h index b1e8068a8..c5e52d1dd 100644 --- a/Eigen/src/Core/Transpose.h +++ b/Eigen/src/Core/Transpose.h @@ -195,7 +195,7 @@ struct ei_inplace_transpose_selector; template<typename MatrixType> struct ei_inplace_transpose_selector<MatrixType,true> { // square matrix static void run(MatrixType& m) { - m.template part<StrictlyUpperTriangular>().swap(m.transpose()); + m.template triangularView<StrictlyUpperTriangular>().swap(m.transpose()); } }; @@ -203,7 +203,7 @@ template<typename MatrixType> struct ei_inplace_transpose_selector<MatrixType,false> { // non square matrix static void run(MatrixType& m) { if (m.rows()==m.cols()) - m.template part<StrictlyUpperTriangular>().swap(m.transpose()); + m.template triangularView<StrictlyUpperTriangular>().swap(m.transpose()); else m = m.transpose().eval(); } diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h new file mode 100644 index 000000000..8131ef323 --- /dev/null +++ b/Eigen/src/Core/TriangularMatrix.h @@ -0,0 +1,637 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008 Benoit Jacob <jacob.benoit.1@gmail.com> +// Copyright (C) 2008-2009 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_TRIANGULARMATRIX_H +#define EIGEN_TRIANGULARMATRIX_H + +/** \nonstableyet + * \class TriangularBase + * + * \brief Expression of a triangular matrix extracted from a given matrix + * + * \param MatrixType the type of the object in which we are taking the triangular part + * \param Mode the kind of triangular matrix expression to construct. Can be UpperTriangular, + * LowerTriangular, UpperSelfadjoint, or LowerSelfadjoint. This is in fact a bit field; + * it must have either UpperBit or LowerBit, and additionnaly it may have either + * TraingularBit or SelfadjointBit. + * + * This class represents an expression of the upper or lower triangular part of + * a square matrix, possibly with a further assumption on the diagonal. It is the return type + * of MatrixBase::part() and most of the time this is the only way it is used. + * + * \sa MatrixBase::part() + */ +template<typename Derived> class TriangularBase : public MultiplierBase<Derived> +{ + public: + + enum { + Mode = ei_traits<Derived>::Mode, + CoeffReadCost = ei_traits<Derived>::CoeffReadCost, + RowsAtCompileTime = ei_traits<Derived>::RowsAtCompileTime, + ColsAtCompileTime = ei_traits<Derived>::ColsAtCompileTime, + MaxRowsAtCompileTime = ei_traits<Derived>::MaxRowsAtCompileTime, + MaxColsAtCompileTime = ei_traits<Derived>::MaxColsAtCompileTime + }; + typedef typename ei_traits<Derived>::Scalar Scalar; + + inline TriangularBase() { ei_assert(ei_are_flags_consistent<Mode>::ret); } + + inline int rows() const { return derived().rows(); } + inline int cols() const { return derived().cols(); } + inline int stride() const { return derived().stride(); } + + inline Scalar coeff(int row, int col) const { return derived().coeff(row,col); } + inline Scalar& coeffRef(int row, int col) { return derived().coeffRef(row,col); } + + /** \see MatrixBase::copyCoeff(row,col) + */ + template<typename Other> + EIGEN_STRONG_INLINE void copyCoeff(int row, int col, Other& other) + { + derived().coeffRef(row, col) = other.coeff(row, col); + } + + inline Scalar operator()(int row, int col) const + { + check_coordinates(row, col); + return coeff(row,col); + } + inline Scalar& operator()(int row, int col) + { + check_coordinates(row, col); + return coeffRef(row,col); + } + + #ifndef EIGEN_PARSED_BY_DOXYGEN + inline const Derived& derived() const { return *static_cast<const Derived*>(this); } + inline Derived& derived() { return *static_cast<Derived*>(this); } + #endif // not EIGEN_PARSED_BY_DOXYGEN + + protected: + + void check_coordinates(int row, int col) + { + ei_assert(col>0 && col<cols() && row>0 && row<rows()); + ei_assert( (Mode==UpperTriangular && col>=row) + || (Mode==LowerTriangular && col<=row) + || (Mode==StrictlyUpperTriangular && col>row) + || (Mode==StrictlyLowerTriangular && col<row)); + } + + void check_coordinates_internal(int row, int col) + { + #ifdef EIGEN_INTERNAL_DEBUGGING + check_coordinates(row, col); + #endif + } + +}; + + +/** \class TriangularView + * \nonstableyet + * + * \brief Expression of a triangular part of a dense matrix + * + * \param MatrixType the type of the dense matrix storing the coefficients + * + * This class is an expression of a triangular part of a matrix with given dense + * storage of the coefficients. It is the return type of MatrixBase::triangularPart() + * and most of the time this is the only way that it is used. + * + * \sa class TriangularBase, MatrixBase::triangularPart(), class DiagonalWrapper + */ +template<typename MatrixType, unsigned int _Mode> +struct ei_traits<TriangularView<MatrixType, _Mode> > : ei_traits<MatrixType> +{ + typedef typename ei_nested<MatrixType>::type MatrixTypeNested; + typedef typename ei_unref<MatrixTypeNested>::type _MatrixTypeNested; + typedef MatrixType ExpressionType; + enum { + Mode = _Mode, + Flags = (_MatrixTypeNested::Flags & (HereditaryBits) & (~(PacketAccessBit | DirectAccessBit | LinearAccessBit))) | Mode, + CoeffReadCost = _MatrixTypeNested::CoeffReadCost + }; +}; + +template<typename MatrixType, unsigned int _Mode> class TriangularView + : public TriangularBase<TriangularView<MatrixType, _Mode> > +{ + public: + + typedef TriangularBase<TriangularView> Base; + typedef typename ei_traits<TriangularView>::Scalar Scalar; + typedef typename MatrixType::PlainMatrixType PlainMatrixType; + + enum { + Mode = _Mode, + TransposeMode = (Mode & UpperTriangularBit ? LowerTriangularBit : 0) + | (Mode & LowerTriangularBit ? UpperTriangularBit : 0) + | (Mode & (ZeroDiagBit | UnitDiagBit)) + }; + + inline TriangularView(const MatrixType& matrix) : m_matrix(matrix) + { ei_assert(ei_are_flags_consistent<Mode>::ret); } + + inline int rows() const { return m_matrix.rows(); } + inline int cols() const { return m_matrix.cols(); } + inline int stride() const { return m_matrix.stride(); } + + /** \sa MatrixBase::operator+=() */ + template<typename Other> TriangularView& operator+=(const Other& other) { return *this = m_matrix + other; } + /** \sa MatrixBase::operator-=() */ + template<typename Other> TriangularView& operator-=(const Other& other) { return *this = m_matrix - other; } + /** \sa MatrixBase::operator*=() */ + TriangularView& operator*=(const typename ei_traits<MatrixType>::Scalar& other) { return *this = m_matrix * other; } + /** \sa MatrixBase::operator/=() */ + TriangularView& operator/=(const typename ei_traits<MatrixType>::Scalar& other) { return *this = m_matrix / other; } + + /** \sa MatrixBase::fill() */ + void fill(const Scalar& value) { setConstant(value); } + /** \sa MatrixBase::setConstant() */ + TriangularView& setConstant(const Scalar& value) + { return *this = MatrixType::Constant(rows(), cols(), value); } + /** \sa MatrixBase::setZero() */ + TriangularView& setZero() { return setConstant(Scalar(0)); } + /** \sa MatrixBase::setOnes() */ + TriangularView& setOnes() { return setConstant(Scalar(1)); } + + /** \sa MatrixBase::coeff() + * \warning the coordinates must fit into the referenced triangular part + */ + inline Scalar coeff(int row, int col) const + { + Base::check_coordinates_internal(row, col); + return m_matrix.coeff(row, col); + } + + /** \sa MatrixBase::coeffRef() + * \warning the coordinates must fit into the referenced triangular part + */ + inline Scalar& coeffRef(int row, int col) + { + Base::check_coordinates_internal(row, col); + return m_matrix.const_cast_derived().coeffRef(row, col); + } + + /** \internal */ + const MatrixType& _expression() const { return m_matrix; } + + /** Assigns a triangular matrix to a triangular part of a dense matrix */ + template<typename OtherDerived> + TriangularView& operator=(const TriangularBase<OtherDerived>& other); + + template<typename OtherDerived> + TriangularView& operator=(const MatrixBase<OtherDerived>& other); + + TriangularView& operator=(const TriangularView& other) + { return *this = other._expression(); } + + template<typename OtherDerived> + void lazyAssign(const TriangularBase<OtherDerived>& other); + + template<typename OtherDerived> + void lazyAssign(const MatrixBase<OtherDerived>& other); + + + /** \sa MatrixBase::adjoint() */ + inline TriangularView<NestByValue<typename MatrixType::AdjointReturnType>,TransposeMode> adjoint() + { return m_matrix.adjoint().nestByValue(); } + /** \sa MatrixBase::adjoint() const */ + const inline TriangularView<NestByValue<typename MatrixType::AdjointReturnType>,TransposeMode> adjoint() const + { return m_matrix.adjoint().nestByValue(); } + + /** \sa MatrixBase::transpose() */ + inline TriangularView<NestByValue<Transpose<MatrixType> >,TransposeMode> transpose() + { return m_matrix.transpose().nestByValue(); } + /** \sa MatrixBase::transpose() const */ + const inline TriangularView<NestByValue<Transpose<MatrixType> >,TransposeMode> transpose() const + { return m_matrix.transpose().nestByValue(); } + + PlainMatrixType toDense() const + { + PlainMatrixType res(rows(), cols()); + res = *this; + return res; + } + + template<typename OtherDerived> + typename ei_plain_matrix_type_column_major<OtherDerived>::type + solve(const MatrixBase<OtherDerived>& other) const; + + template<typename OtherDerived> + void solveInPlace(const MatrixBase<OtherDerived>& other) const; + + template<typename OtherDerived> + void swap(const TriangularBase<OtherDerived>& other) + { + TriangularView<SwapWrapper<MatrixType>,Mode>(const_cast<MatrixType&>(m_matrix)).lazyAssign(other.derived()); + } + + template<typename OtherDerived> + void swap(const MatrixBase<OtherDerived>& other) + { + TriangularView<SwapWrapper<MatrixType>,Mode>(const_cast<MatrixType&>(m_matrix)).lazyAssign(other.derived()); + } + + protected: + + const typename MatrixType::Nested m_matrix; +}; + +/*************************************************************************** +* Implementation of triangular evaluation/assignment +***************************************************************************/ + +template<typename Derived1, typename Derived2, unsigned int Mode, int UnrollCount, bool ClearOpposite> +struct ei_part_assignment_impl +{ + enum { + col = (UnrollCount-1) / Derived1::RowsAtCompileTime, + row = (UnrollCount-1) % Derived1::RowsAtCompileTime + }; + + inline static void run(Derived1 &dst, const Derived2 &src) + { + ei_part_assignment_impl<Derived1, Derived2, Mode, UnrollCount-1, ClearOpposite>::run(dst, src); + + if(Mode == SelfAdjoint) + { + if(row == col) + dst.coeffRef(row, col) = ei_real(src.coeff(row, col)); + else if(row < col) + dst.coeffRef(col, row) = ei_conj(dst.coeffRef(row, col) = src.coeff(row, col)); + } + else + { + ei_assert( Mode == UpperTriangular || Mode == LowerTriangular + || Mode == StrictlyUpperTriangular || Mode == StrictlyLowerTriangular + || Mode == UnitUpperTriangular || Mode == UnitLowerTriangular); + if((Mode == UpperTriangular && row <= col) + || (Mode == LowerTriangular && row >= col) + || (Mode == StrictlyUpperTriangular && row < col) + || (Mode == StrictlyLowerTriangular && row > col) + || (Mode == UnitUpperTriangular && row < col) + || (Mode == UnitLowerTriangular && row > col)) + dst.copyCoeff(row, col, src); + else if(ClearOpposite) + { + if (Mode&UnitDiagBit && row==col) + dst.coeffRef(row, col) = 1; + else + dst.coeffRef(row, col) = 0; + } + } + } +}; +template<typename Derived1, typename Derived2, unsigned int Mode, bool ClearOpposite> +struct ei_part_assignment_impl<Derived1, Derived2, Mode, 1, ClearOpposite> +{ + inline static void run(Derived1 &dst, const Derived2 &src) + { + if(Mode&UnitDiagBit) + { + if(ClearOpposite) + dst.coeffRef(0, 0) = 1; + } + else if(!(Mode & ZeroDiagBit)) + dst.copyCoeff(0, 0, src); + } +}; +// prevent buggy user code from causing an infinite recursion +template<typename Derived1, typename Derived2, unsigned int Mode, bool ClearOpposite> +struct ei_part_assignment_impl<Derived1, Derived2, Mode, 0, ClearOpposite> +{ + inline static void run(Derived1 &, const Derived2 &) {} +}; + +template<typename Derived1, typename Derived2, bool ClearOpposite> +struct ei_part_assignment_impl<Derived1, Derived2, UpperTriangular, Dynamic, ClearOpposite> +{ + inline static void run(Derived1 &dst, const Derived2 &src) + { + for(int j = 0; j < dst.cols(); ++j) + { + for(int i = 0; i <= j; ++i) + dst.copyCoeff(i, j, src); + if (ClearOpposite) + for(int i = j+1; i < dst.rows(); ++i) + dst.coeffRef(i, j) = 0; + } + } +}; +template<typename Derived1, typename Derived2, bool ClearOpposite> +struct ei_part_assignment_impl<Derived1, Derived2, LowerTriangular, Dynamic, ClearOpposite> +{ + inline static void run(Derived1 &dst, const Derived2 &src) + { + for(int j = 0; j < dst.cols(); ++j) + { + for(int i = j; i < dst.rows(); ++i) + dst.copyCoeff(i, j, src); + if (ClearOpposite) + for(int i = 0; i < j; ++i) + dst.coeffRef(i, j) = 0; + } + } +}; + +template<typename Derived1, typename Derived2, bool ClearOpposite> +struct ei_part_assignment_impl<Derived1, Derived2, StrictlyUpperTriangular, Dynamic, ClearOpposite> +{ + inline static void run(Derived1 &dst, const Derived2 &src) + { + for(int j = 0; j < dst.cols(); ++j) + { + for(int i = 0; i < j; ++i) + dst.copyCoeff(i, j, src); + if (ClearOpposite) + for(int i = j; i < dst.rows(); ++i) + dst.coeffRef(i, j) = 0; + } + } +}; +template<typename Derived1, typename Derived2, bool ClearOpposite> +struct ei_part_assignment_impl<Derived1, Derived2, StrictlyLowerTriangular, Dynamic, ClearOpposite> +{ + inline static void run(Derived1 &dst, const Derived2 &src) + { + for(int j = 0; j < dst.cols(); ++j) + { + for(int i = j+1; i < dst.rows(); ++i) + dst.copyCoeff(i, j, src); + if (ClearOpposite) + for(int i = 0; i <= j; ++i) + dst.coeffRef(i, j) = 0; + } + } +}; + +template<typename Derived1, typename Derived2, bool ClearOpposite> +struct ei_part_assignment_impl<Derived1, Derived2, UnitUpperTriangular, Dynamic, ClearOpposite> +{ + inline static void run(Derived1 &dst, const Derived2 &src) + { + for(int j = 0; j < dst.cols(); ++j) + { + for(int i = 0; i < j; ++i) + dst.copyCoeff(i, j, src); + if (ClearOpposite) + { + for(int i = j+1; i < dst.rows(); ++i) + dst.coeffRef(i, j) = 0; + dst.coeffRef(j, j) = 1; + } + } + } +}; +template<typename Derived1, typename Derived2, bool ClearOpposite> +struct ei_part_assignment_impl<Derived1, Derived2, UnitLowerTriangular, Dynamic, ClearOpposite> +{ + inline static void run(Derived1 &dst, const Derived2 &src) + { + for(int j = 0; j < dst.cols(); ++j) + { + for(int i = j+1; i < dst.rows(); ++i) + dst.copyCoeff(i, j, src); + if (ClearOpposite) + { + for(int i = 0; i < j; ++i) + dst.coeffRef(i, j) = 0; + dst.coeffRef(j, j) = 1; + } + } + } +}; + +// selfadjoint to dense matrix +template<typename Derived1, typename Derived2, bool ClearOpposite> +struct ei_part_assignment_impl<Derived1, Derived2, SelfAdjoint, Dynamic, ClearOpposite> +{ + inline static void run(Derived1 &dst, const Derived2 &src) + { + for(int j = 0; j < dst.cols(); ++j) + { + for(int i = 0; i < j; ++i) + dst.coeffRef(j, i) = ei_conj(dst.coeffRef(i, j) = src.coeff(i, j)); + dst.coeffRef(j, j) = ei_real(src.coeff(j, j)); + } + } +}; + +// FIXME should we keep that possibility +template<typename MatrixType, unsigned int Mode> +template<typename OtherDerived> +inline TriangularView<MatrixType, Mode>& +TriangularView<MatrixType, Mode>::operator=(const MatrixBase<OtherDerived>& other) +{ + if(OtherDerived::Flags & EvalBeforeAssigningBit) + { + typename OtherDerived::PlainMatrixType other_evaluated(other.rows(), other.cols()); + other_evaluated.template triangularView<Mode>().lazyAssign(other.derived()); + lazyAssign(other_evaluated); + } + else + lazyAssign(other.derived()); + return *this; +} + +// FIXME should we keep that possibility +template<typename MatrixType, unsigned int Mode> +template<typename OtherDerived> +void TriangularView<MatrixType, Mode>::lazyAssign(const MatrixBase<OtherDerived>& other) +{ + const bool unroll = MatrixType::SizeAtCompileTime * ei_traits<OtherDerived>::CoeffReadCost / 2 + <= EIGEN_UNROLLING_LIMIT; + ei_assert(m_matrix.rows() == other.rows() && m_matrix.cols() == other.cols()); + + ei_part_assignment_impl + <MatrixType, OtherDerived, int(Mode), + unroll ? int(MatrixType::SizeAtCompileTime) : Dynamic, + false // do not change the opposite triangular part + >::run(m_matrix.const_cast_derived(), other.derived()); +} + + + +template<typename MatrixType, unsigned int Mode> +template<typename OtherDerived> +inline TriangularView<MatrixType, Mode>& +TriangularView<MatrixType, Mode>::operator=(const TriangularBase<OtherDerived>& other) +{ + ei_assert(Mode == OtherDerived::Mode); + if(ei_traits<OtherDerived>::Flags & EvalBeforeAssigningBit) + { + typename OtherDerived::PlainMatrixType other_evaluated(other.rows(), other.cols()); + other_evaluated.template triangularView<Mode>().lazyAssign(other.derived()); + lazyAssign(other_evaluated); + } + else + lazyAssign(other.derived()); + return *this; +} + +template<typename MatrixType, unsigned int Mode> +template<typename OtherDerived> +void TriangularView<MatrixType, Mode>::lazyAssign(const TriangularBase<OtherDerived>& other) +{ + const bool unroll = MatrixType::SizeAtCompileTime * ei_traits<OtherDerived>::CoeffReadCost / 2 + <= EIGEN_UNROLLING_LIMIT; + ei_assert(m_matrix.rows() == other.rows() && m_matrix.cols() == other.cols()); + + ei_part_assignment_impl + <MatrixType, OtherDerived, int(Mode), + unroll ? int(MatrixType::SizeAtCompileTime) : Dynamic, + false // preserve the opposite triangular part + >::run(m_matrix.const_cast_derived(), other.derived()._expression()); +} + +/*************************************************************************** +* Implementation of MatrixBase methods +***************************************************************************/ + +/** Assigns a triangular or selfadjoint matrix to a dense matrix. + * If the matrix is triangular, the opposite part is set to zero. */ +template<typename Derived> +template<typename TriangularDerived> +Derived& MatrixBase<Derived>::operator=(const TriangularBase<TriangularDerived> &other) +{ + if(ei_traits<TriangularDerived>::Flags & EvalBeforeAssigningBit) + { + typename TriangularDerived::PlainMatrixType other_evaluated(other.rows(), other.cols()); + other_evaluated.lazyAssign(other); + this->swap(other_evaluated); + } + else + lazyAssign(other.derived()); + return derived(); +} + +/** Assigns a triangular or selfadjoint matrix to a dense matrix. + * If the matrix is triangular, the opposite part is set to zero. */ +template<typename Derived> +template<typename TriangularDerived> +Derived& MatrixBase<Derived>::lazyAssign(const TriangularBase<TriangularDerived> &other) +{ + const bool unroll = Derived::SizeAtCompileTime * TriangularDerived::CoeffReadCost / 2 + <= EIGEN_UNROLLING_LIMIT; + ei_assert(this->rows() == other.rows() && this->cols() == other.cols()); + + ei_part_assignment_impl + <Derived, typename ei_traits<TriangularDerived>::ExpressionType, TriangularDerived::Mode, + unroll ? int(Derived::SizeAtCompileTime) : Dynamic, + true // clear the opposite triangular part + >::run(this->const_cast_derived(), other.derived()._expression()); + return derived(); +} + +/** \deprecated use MatrixBase::triangularView() */ +template<typename Derived> +template<unsigned int Mode> +EIGEN_DEPRECATED const TriangularView<Derived, Mode> MatrixBase<Derived>::part() const +{ + return derived(); +} + +/** \deprecated use MatrixBase::triangularView() */ +template<typename Derived> +template<unsigned int Mode> +EIGEN_DEPRECATED TriangularView<Derived, Mode> MatrixBase<Derived>::part() +{ + return derived(); +} + +/** \nonstableyet + * \returns an expression of a triangular view extracted from the current matrix + * + * The parameter \a Mode can have the following values: \c UpperTriangular, \c StrictlyUpperTriangular, \c UnitUpperTriangular, + * \c LowerTriangular, \c StrictlyLowerTriangular, \c UnitLowerTriangular. + * + * Example: \include MatrixBase_extract.cpp + * Output: \verbinclude MatrixBase_extract.out + * + * \sa class TriangularView + */ +template<typename Derived> +template<unsigned int Mode> +TriangularView<Derived, Mode> MatrixBase<Derived>::triangularView() +{ + return derived(); +} + +/** This is the const version of MatrixBase::triangularView() */ +template<typename Derived> +template<unsigned int Mode> +const TriangularView<Derived, Mode> MatrixBase<Derived>::triangularView() const +{ + return derived(); +} + +/** \returns true if *this is approximately equal to an upper triangular matrix, + * within the precision given by \a prec. + * + * \sa isLowerTriangular(), extract(), part(), marked() + */ +template<typename Derived> +bool MatrixBase<Derived>::isUpperTriangular(RealScalar prec) const +{ + if(cols() != rows()) return false; + RealScalar maxAbsOnUpperTriangularPart = 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 > maxAbsOnUpperTriangularPart) maxAbsOnUpperTriangularPart = absValue; + } + for(int j = 0; j < cols()-1; ++j) + for(int i = j+1; i < rows(); ++i) + if(!ei_isMuchSmallerThan(coeff(i, j), maxAbsOnUpperTriangularPart, 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 isUpperTriangular(), extract(), part(), marked() + */ +template<typename Derived> +bool MatrixBase<Derived>::isLowerTriangular(RealScalar prec) const +{ + if(cols() != rows()) return false; + RealScalar maxAbsOnLowerTriangularPart = 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 > maxAbsOnLowerTriangularPart) maxAbsOnLowerTriangularPart = absValue; + } + for(int j = 1; j < cols(); ++j) + for(int i = 0; i < j; ++i) + if(!ei_isMuchSmallerThan(coeff(i, j), maxAbsOnLowerTriangularPart, prec)) return false; + return true; +} + +#endif // EIGEN_TRIANGULARMATRIX_H diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index af21c190f..db58e4982 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -57,8 +57,9 @@ template<typename MatrixType, typename DiagonalType, int ProductOrder> class Dia template<typename MatrixType, int Index> class Diagonal; template<typename MatrixType, int PacketAccess = AsRequested> class Map; -template<typename MatrixType, unsigned int Mode> class Part; -template<typename MatrixType, unsigned int Mode> class Extract; +template<typename Derived> class TriangularBase; +template<typename MatrixType, unsigned int Mode> class TriangularView; +template<typename MatrixType, unsigned int Mode> class SelfAdjointView; template<typename ExpressionType> class Cwise; template<typename ExpressionType> class WithFormat; template<typename MatrixType> struct CommaInitializer; @@ -119,7 +120,7 @@ template<typename MatrixType> class LU; template<typename MatrixType> class PartialLU; template<typename MatrixType> class QR; template<typename MatrixType> class SVD; -template<typename MatrixType> class LLT; +template<typename MatrixType, int UpLo = LowerTriangular> class LLT; template<typename MatrixType> class LDLT; // Geometry module: diff --git a/Eigen/src/QR/QR.h b/Eigen/src/QR/QR.h index 90751dd42..5883eed65 100644 --- a/Eigen/src/QR/QR.h +++ b/Eigen/src/QR/QR.h @@ -131,7 +131,7 @@ template<typename MatrixType> class QR } /** \returns a read-only expression of the matrix R of the actual the QR decomposition */ - const Part<NestByValue<MatrixRBlockType>, UpperTriangular> + const TriangularView<NestByValue<MatrixRBlockType>, UpperTriangular> matrixR(void) const { ei_assert(m_isInitialized && "QR is not initialized."); |