mirror of
https://gitlab.com/libeigen/eigen.git
synced 2026-04-10 11:34:33 +08:00
Added Triangular expression to extract upper or lower (strictly or not)
part of a matrix. Triangular also provide an optimised method for forward and backward substitution. Further optimizations regarding assignments and products might come later. Updated determinant() to take into account triangular matrices. Started the QR module with a QR decompostion algorithm. Help needed to build a QR algorithm (eigen solver) based on it.
This commit is contained in:
Gael Guennebaud
347
Eigen/src/Core/Triangular.h
Executable file
347
Eigen/src/Core/Triangular.h
Executable file
@@ -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 <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_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 = (<any expression>).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() = <any expression>;
|
||||
* \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 <any \c expression> has to be a squared matrix.
|
||||
*
|
||||
* \sa MatrixBase::upper(), MatrixBase::lower(), class TriangularProduct
|
||||
*/
|
||||
template<int Mode, typename MatrixType>
|
||||
struct ei_traits<Triangular<Mode, MatrixType> >
|
||||
{
|
||||
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<int Mode, typename MatrixType> class Triangular
|
||||
: public MatrixBase<Triangular<Mode,MatrixType> >
|
||||
{
|
||||
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<MatrixType> > transpose()
|
||||
{
|
||||
return Triangular<(Upper | Lower) xor Mode, Transpose<MatrixType> >((m_matrix.transpose()));
|
||||
}
|
||||
|
||||
/** Overloaded to keep a Triangular expression */
|
||||
const Triangular<(Upper | Lower) xor Mode, Transpose<MatrixType> > transpose() const
|
||||
{
|
||||
return Triangular<(Upper | Lower) xor Mode, Transpose<MatrixType> >((m_matrix.transpose()));
|
||||
}
|
||||
|
||||
#if 0
|
||||
|
||||
template<typename OtherDerived>
|
||||
Triangular& operator=(const MatrixBase<OtherDerived>& other);
|
||||
|
||||
/** Overloaded to provide optimal evaluation loops */
|
||||
template<typename OtherDerived>
|
||||
Triangular& operator +=(const MatrixBase<OtherDerived>& other)
|
||||
{
|
||||
return *this = m_matrix + other;
|
||||
}
|
||||
|
||||
/** Overloaded to provide optimal evaluation loops */
|
||||
template<typename OtherDerived>
|
||||
Triangular& operator *=(const MatrixBase<OtherDerived>& other)
|
||||
{
|
||||
return *this = this->lazyProduct(other).eval();
|
||||
}
|
||||
|
||||
/** Optimized triangular matrix - matrix product */
|
||||
template<typename OtherDerived>
|
||||
TriangularProduct<Mode, MatrixType, OtherDerived> lazyProduct(const MatrixBase<Scalar, OtherDerived>& other) const
|
||||
{
|
||||
return TriangularProduct<Mode,MatrixType,OtherDerived>(m_matrix, other.ref());
|
||||
}
|
||||
|
||||
/** Optimized triangular matrix - matrix product */
|
||||
template<typename OtherDerived>
|
||||
Eval<TriangularProduct<Mode, MatrixType, OtherDerived> > operator * (const MatrixBase<Scalar, OtherDerived>& other) const
|
||||
{
|
||||
return this->lazyProduct(other).eval();
|
||||
}
|
||||
|
||||
/** Optimized matrix - triangular matrix product */
|
||||
template<typename OtherDerived>
|
||||
friend Eval<Transpose<TriangularProduct<0x1 xor Mode, Transpose<MatRef>, Transpose<OtherDerived> > > >
|
||||
operator * (const MatrixBase<Scalar, OtherDerived>& other, const Triangular<Mode,MatrixType>& 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>
|
||||
typename OtherDerived::Eval inverseProduct(const MatrixBase<OtherDerived>& other) const
|
||||
{
|
||||
assert(_cols() == other.rows());
|
||||
assert(!(Flags & NullDiagBit));
|
||||
|
||||
typename OtherDerived::Eval res(other.rows(), other.cols());
|
||||
|
||||
for (int c=0 ; c<other.cols() ; ++c)
|
||||
{
|
||||
if (Flags & Lower)
|
||||
{
|
||||
// forward substitution
|
||||
if (Flags & UnitDiagBit)
|
||||
res.col(c)[0] = other.col(c)[0];
|
||||
else
|
||||
res.col(c)[0] = other.col(c)[0]/_coeff(0, 0);
|
||||
for (int i=1 ; i<_rows() ; ++i)
|
||||
{
|
||||
Scalar tmp = other.col(c)[i];
|
||||
for (int j = 0 ; j < i ; ++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);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
// backward substitution
|
||||
if (Flags & UnitDiagBit)
|
||||
res.col(c)[_cols()-1] = other.col(c)[_cols()-1];
|
||||
else
|
||||
res.col(c)[_cols()-1] = other.col(c)[_cols()-1]/_coeff(_rows()-1, _cols()-1);
|
||||
for (int i=_rows()-2 ; i>=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<typename Derived>
|
||||
Triangular<Upper, Derived> MatrixBase<Derived>::upper(void)
|
||||
{
|
||||
return Triangular<Upper,Derived>(derived());
|
||||
}
|
||||
|
||||
/** This is the const version of upper(). */
|
||||
template<typename Derived>
|
||||
const Triangular<Upper, Derived> MatrixBase<Derived>::upper(void) const
|
||||
{
|
||||
return Triangular<Upper,Derived>(derived());
|
||||
}
|
||||
|
||||
/** \returns an expression of a lower triangular matrix
|
||||
*
|
||||
* \sa isLower(), lowerWithUnitDiag(), lowerWithNullDiag(), upper()
|
||||
*/
|
||||
template<typename Derived>
|
||||
Triangular<Lower, Derived> MatrixBase<Derived>::lower(void)
|
||||
{
|
||||
return Triangular<Lower,Derived>(derived());
|
||||
}
|
||||
|
||||
/** This is the const version of lower().*/
|
||||
template<typename Derived>
|
||||
const Triangular<Lower, Derived> MatrixBase<Derived>::lower(void) const
|
||||
{
|
||||
return Triangular<Lower,Derived>(derived());
|
||||
}
|
||||
|
||||
/** \returns an expression of a upper triangular matrix with a unit diagonal
|
||||
*
|
||||
* \sa upper(), lowerWithUnitDiagBit()
|
||||
*/
|
||||
template<typename Derived>
|
||||
const Triangular<Upper|UnitDiagBit, Derived> MatrixBase<Derived>::upperWithUnitDiag(void) const
|
||||
{
|
||||
return Triangular<Upper|UnitDiagBit, Derived>(derived());
|
||||
}
|
||||
|
||||
/** \returns an expression of a strictly upper triangular matrix (diagonal==zero)
|
||||
* FIXME could also be called strictlyUpper() or upperStrict()
|
||||
*
|
||||
* \sa upper(), lowerWithNullDiag()
|
||||
*/
|
||||
template<typename Derived>
|
||||
const Triangular<Upper|NullDiagBit, Derived> MatrixBase<Derived>::upperWithNullDiag(void) const
|
||||
{
|
||||
return Triangular<Upper|NullDiagBit, Derived>(derived());
|
||||
}
|
||||
|
||||
/** \returns an expression of a lower triangular matrix with a unit diagonal
|
||||
*
|
||||
* \sa lower(), upperWithUnitDiag()
|
||||
*/
|
||||
template<typename Derived>
|
||||
const Triangular<Lower|UnitDiagBit, Derived> MatrixBase<Derived>::lowerWithUnitDiag(void) const
|
||||
{
|
||||
return Triangular<Lower|UnitDiagBit, Derived>(derived());
|
||||
}
|
||||
|
||||
/** \returns an expression of a strictly lower triangular matrix (diagonal==zero)
|
||||
* FIXME could also be called strictlyLower() or lowerStrict()
|
||||
*
|
||||
* \sa lower(), upperWithNullDiag()
|
||||
*/
|
||||
template<typename Derived>
|
||||
const Triangular<Lower|NullDiagBit, Derived> MatrixBase<Derived>::lowerWithNullDiag(void) const
|
||||
{
|
||||
return Triangular<Lower|NullDiagBit, Derived>(derived());
|
||||
}
|
||||
|
||||
/** \returns true if *this is approximately equal to an upper triangular matrix,
|
||||
* within the precision given by \a prec.
|
||||
*
|
||||
* \sa isLower(), upper()
|
||||
*/
|
||||
template<typename Derived>
|
||||
bool MatrixBase<Derived>::isUpper(RealScalar prec) const
|
||||
{
|
||||
if(cols() != rows()) return false;
|
||||
RealScalar maxAbsOnUpperPart = 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 > 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<typename Derived>
|
||||
bool MatrixBase<Derived>::isLower(RealScalar prec) const
|
||||
{
|
||||
if(cols() != rows()) return false;
|
||||
RealScalar maxAbsOnLowerPart = 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 > 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
|
||||
Reference in New Issue
Block a user