devtools/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-04-10 11:34:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

TriangularMatrix: extend to rectangular matrices

Browse Source
This commit is contained in:
Benoit Jacob
2009-11-19 17:07:55 -05:00
parent 2275f98d7b
commit a20a744adc
2 changed files with 193 additions and 83 deletions
Download Patch File Download Diff File Expand all files Collapse all files

93
Eigen/src/Core/TriangularMatrix.h
View File

@@ -26,22 +26,11 @@
#ifndef EIGEN_TRIANGULARMATRIX_H
#define EIGEN_TRIANGULARMATRIX_H
/** \nonstableyet
/** \internal
*
* \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()
* \brief Base class for triangular part in a matrix
*/
template<typename Derived> class TriangularBase : public AnyMatrixBase<Derived>
{
@@ -115,19 +104,21 @@ template<typename Derived> class TriangularBase : public AnyMatrixBase<Derived>
};
/** \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<typename MatrixType, unsigned int _Mode>
struct ei_traits<TriangularView<MatrixType, _Mode> > : ei_traits<MatrixType>
@@ -155,7 +146,7 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView
typedef TriangularBase<TriangularView> Base;
typedef typename ei_traits<TriangularView>::Scalar Scalar;
typedef _MatrixType MatrixType;
typedef typename MatrixType::PlainMatrixType PlainMatrixType;
typedef typename MatrixType::PlainMatrixType DenseMatrixType;
typedef typename MatrixType::Nested MatrixTypeNested;
typedef typename ei_cleantype<MatrixTypeNested>::type _MatrixTypeNested;
@@ -244,9 +235,9 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView
inline const TriangularView<NestByValue<Transpose<MatrixType> >,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<typename Derived1, typename Derived2, unsigned int Mode, bool ClearOpposite>
struct ei_triangular_assignment_selector<Derived1, Derived2, Mode, 1, ClearOpposite>
{
@@ -365,6 +357,7 @@ struct ei_triangular_assignment_selector<Derived1, Derived2, Mode, 1, ClearOppos
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_triangular_assignment_selector<Derived1, Derived2, Mode, 0, ClearOpposite>
@@ -379,14 +372,16 @@ struct ei_triangular_assignment_selector<Derived1, Derived2, UpperTriangular, Dy
{
for(int j = 0; j < dst.cols(); ++j)
{
for(int i = 0; i <= j; ++i)
int maxi = std::min(j, dst.rows()-1);
for(int i = 0; i <= maxi; ++i)
dst.copyCoeff(i, j, src);
if (ClearOpposite)
for(int i = j+1; i < dst.rows(); ++i)
for(int i = maxi+1; i < dst.rows(); ++i)
dst.coeffRef(i, j) = 0;
}
}
};
template<typename Derived1, typename Derived2, bool ClearOpposite>
struct ei_triangular_assignment_selector<Derived1, Derived2, LowerTriangular, Dynamic, ClearOpposite>
{
@@ -396,8 +391,9 @@ struct ei_triangular_assignment_selector<Derived1, Derived2, LowerTriangular, Dy
{
for(int i = j; i < dst.rows(); ++i)
dst.copyCoeff(i, j, src);
int maxi = std::min(j, dst.rows());
if (ClearOpposite)
for(int i = 0; i < j; ++i)
for(int i = 0; i < maxi; ++i)
dst.coeffRef(i, j) = 0;
}
}
@@ -410,14 +406,16 @@ struct ei_triangular_assignment_selector<Derived1, Derived2, StrictlyUpperTriang
{
for(int j = 0; j < dst.cols(); ++j)
{
for(int i = 0; i < j; ++i)
int maxi = std::min(j, dst.rows());
for(int i = 0; i < maxi; ++i)
dst.copyCoeff(i, j, src);
if (ClearOpposite)
for(int i = j; i < dst.rows(); ++i)
for(int i = maxi; i < dst.rows(); ++i)
dst.coeffRef(i, j) = 0;
}
}
};
template<typename Derived1, typename Derived2, bool ClearOpposite>
struct ei_triangular_assignment_selector<Derived1, Derived2, StrictlyLowerTriangular, Dynamic, ClearOpposite>
{
@@ -427,8 +425,9 @@ struct ei_triangular_assignment_selector<Derived1, Derived2, StrictlyLowerTriang
{
for(int i = j+1; i < dst.rows(); ++i)
dst.copyCoeff(i, j, src);
int maxi = std::min(j, dst.rows());
if (ClearOpposite)
for(int i = 0; i <= j; ++i)
for(int i = 0; i <= maxi; ++i)
dst.coeffRef(i, j) = 0;
}
}
@@ -441,11 +440,12 @@ struct ei_triangular_assignment_selector<Derived1, Derived2, UnitUpperTriangular
{
for(int j = 0; j < dst.cols(); ++j)
{
for(int i = 0; i < j; ++i)
int maxi = std::min(j, dst.rows());
for(int i = 0; i < maxi; ++i)
dst.copyCoeff(i, j, src);
if (ClearOpposite)
{
for(int i = j+1; i < dst.rows(); ++i)
for(int i = maxi+1; i < dst.rows(); ++i)
dst.coeffRef(i, j) = 0;
dst.coeffRef(j, j) = 1;
}
@@ -459,11 +459,12 @@ struct ei_triangular_assignment_selector<Derived1, Derived2, UnitLowerTriangular
{
for(int j = 0; j < dst.cols(); ++j)
{
for(int i = j+1; i < dst.rows(); ++i)
int maxi = std::min(j, dst.rows());
for(int i = maxi+1; i < dst.rows(); ++i)
dst.copyCoeff(i, j, src);
if (ClearOpposite)
{
for(int i = 0; i < j; ++i)
for(int i = 0; i < maxi; ++i)
dst.coeffRef(i, j) = 0;
dst.coeffRef(j, j) = 1;
}
@@ -514,7 +515,7 @@ TriangularView<MatrixType, Mode>::operator=(const TriangularBase<OtherDerived>&
ei_assert(Mode == OtherDerived::Mode);
if(ei_traits<OtherDerived>::Flags & EvalBeforeAssigningBit)
{
typename OtherDerived::PlainMatrixType other_evaluated(other.rows(), other.cols());
typename OtherDerived::DenseMatrixType other_evaluated(other.rows(), other.cols());
other_evaluated.template triangularView<Mode>().lazyAssign(other.derived());
lazyAssign(other_evaluated);
}
@@ -633,17 +634,20 @@ const TriangularView<Derived, Mode> MatrixBase<Derived>::triangularView() const
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)
{
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<Derived>::isUpperTriangular(RealScalar prec) const
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)
@@ -663,9 +666,13 @@ bool MatrixBase<Derived>::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;
}