From a20a744adc5414f9a79ae0637de9c3227c3dd420 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Thu, 19 Nov 2009 17:07:55 -0500 Subject: TriangularMatrix: extend to rectangular matrices --- Eigen/src/Core/TriangularMatrix.h | 93 +++++++++++++++++++++------------------ 1 file changed, 50 insertions(+), 43 deletions(-) (limited to 'Eigen') diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h index e60d57e70..ed26683ea 100644 --- a/Eigen/src/Core/TriangularMatrix.h +++ b/Eigen/src/Core/TriangularMatrix.h @@ -26,22 +26,11 @@ #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. +/** \internal * - * 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. + * \class TriangularBase * - * \sa MatrixBase::part() + * \brief Base class for triangular part in a matrix */ template class TriangularBase : public AnyMatrixBase { @@ -115,19 +104,21 @@ template class TriangularBase : public AnyMatrixBase }; - /** \class TriangularView - * \nonstableyet * - * \brief Expression of a triangular part of a dense matrix + * \brief Base class for triangular part in a matrix * - * \param MatrixType the type of the dense matrix storing the coefficients + * \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 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. + * This class represents a triangular part of a matrix, not necessarily square. Strictly speaking, for rectangular + * matrices one should speak ok "trapezoid" parts. This class is the return type + * of MatrixBase::triangularView() and most of the time this is the only way it is used. * - * \sa class TriangularBase, MatrixBase::triangularPart(), class DiagonalWrapper + * \sa MatrixBase::triangularView() */ template struct ei_traits > : ei_traits @@ -155,7 +146,7 @@ template class TriangularView typedef TriangularBase Base; typedef typename ei_traits::Scalar Scalar; typedef _MatrixType MatrixType; - typedef typename MatrixType::PlainMatrixType PlainMatrixType; + typedef typename MatrixType::PlainMatrixType DenseMatrixType; typedef typename MatrixType::Nested MatrixTypeNested; typedef typename ei_cleantype::type _MatrixTypeNested; @@ -244,9 +235,9 @@ template class TriangularView inline const TriangularView >,TransposeMode> transpose() const { return m_matrix.transpose().nestByValue(); } - PlainMatrixType toDense() const + DenseMatrixType toDenseMatrix() const { - PlainMatrixType res(rows(), cols()); + DenseMatrixType res(rows(), cols()); res = *this; return res; } @@ -351,6 +342,7 @@ struct ei_triangular_assignment_selector } } }; + template struct ei_triangular_assignment_selector { @@ -365,6 +357,7 @@ struct ei_triangular_assignment_selector struct ei_triangular_assignment_selector @@ -379,14 +372,16 @@ struct ei_triangular_assignment_selector struct ei_triangular_assignment_selector { @@ -396,8 +391,9 @@ struct ei_triangular_assignment_selector struct ei_triangular_assignment_selector { @@ -427,8 +425,9 @@ struct ei_triangular_assignment_selector::operator=(const TriangularBase& ei_assert(Mode == OtherDerived::Mode); if(ei_traits::Flags & EvalBeforeAssigningBit) { - typename OtherDerived::PlainMatrixType other_evaluated(other.rows(), other.cols()); + typename OtherDerived::DenseMatrixType other_evaluated(other.rows(), other.cols()); other_evaluated.template triangularView().lazyAssign(other.derived()); lazyAssign(other_evaluated); } @@ -633,17 +634,20 @@ const TriangularView MatrixBase::triangularView() const template bool MatrixBase::isUpperTriangular(RealScalar prec) const { - if(cols() != rows()) return false; RealScalar maxAbsOnUpperTriangularPart = static_cast(-1); for(int j = 0; j < cols(); ++j) - for(int i = 0; i <= j; ++i) + { + int maxi = std::min(j, rows()-1); + for(int i = 0; i <= maxi; ++i) { RealScalar absValue = ei_abs(coeff(i,j)); if(absValue > maxAbsOnUpperTriangularPart) maxAbsOnUpperTriangularPart = absValue; } - for(int j = 0; j < cols()-1; ++j) + } + RealScalar threshold = maxAbsOnUpperTriangularPart * prec; + for(int j = 0; j < cols(); ++j) for(int i = j+1; i < rows(); ++i) - if(!ei_isMuchSmallerThan(coeff(i, j), maxAbsOnUpperTriangularPart, prec)) return false; + if(ei_abs(coeff(i, j)) > threshold) return false; return true; } @@ -655,7 +659,6 @@ bool MatrixBase::isUpperTriangular(RealScalar prec) const template bool MatrixBase::isLowerTriangular(RealScalar prec) const { - if(cols() != rows()) return false; RealScalar maxAbsOnLowerTriangularPart = static_cast(-1); for(int j = 0; j < cols(); ++j) for(int i = j; i < rows(); ++i) @@ -663,9 +666,13 @@ bool MatrixBase::isLowerTriangular(RealScalar prec) const RealScalar absValue = ei_abs(coeff(i,j)); if(absValue > maxAbsOnLowerTriangularPart) maxAbsOnLowerTriangularPart = absValue; } + RealScalar threshold = maxAbsOnLowerTriangularPart * prec; for(int j = 1; j < cols(); ++j) - for(int i = 0; i < j; ++i) - if(!ei_isMuchSmallerThan(coeff(i, j), maxAbsOnLowerTriangularPart, prec)) return false; + { + int maxi = std::min(j, rows()-1); + for(int i = 0; i < maxi; ++i) + if(ei_abs(coeff(i, j)) > threshold) return false; + } return true; } -- cgit v1.2.3