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

* the Upper->UpperTriangular change

Browse Source 
* finally get ei_add_test right
This commit is contained in:
Benoit Jacob
2008-12-20 13:36:12 +00:00
parent 21ab65e4b3
commit 9e00d94543
33 changed files with 181 additions and 179 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
Eigen/src/Core/MatrixBase.h
View File

@@ -470,8 +470,8 @@ template<typename Derived> class MatrixBase
bool isIdentity(RealScalar prec = precision<Scalar>()) const;
bool isDiagonal(RealScalar prec = precision<Scalar>()) const;
bool isUpper(RealScalar prec = precision<Scalar>()) const;
bool isLower(RealScalar prec = precision<Scalar>()) const;
bool isUpperTriangular(RealScalar prec = precision<Scalar>()) const;
bool isLowerTriangular(RealScalar prec = precision<Scalar>()) const;
template<typename OtherDerived>
bool isOrthogonal(const MatrixBase<OtherDerived>& other,

58
Eigen/src/Core/Part.h
View File

@@ -31,8 +31,8 @@
* \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 Upper, StrictlyUpper,
* UnitUpper, Lower, StrictlyLower, UnitLower. This is in fact a bit field; it must have either
* \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.
*
@@ -104,10 +104,10 @@ template<typename MatrixType, unsigned int Mode> class Part
{
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==Upper && col>=row)
|| (Mode==Lower && col<=row)
|| (Mode==StrictlyUpper && col>row)
|| (Mode==StrictlyLower && col<row));
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);
}
@@ -134,8 +134,8 @@ template<typename MatrixType, unsigned int Mode> class Part
/** \returns an expression of a triangular matrix extracted from the current matrix
*
* The parameter \a Mode can have the following values: \c Upper, \c StrictlyUpper, \c UnitUpper,
* \c Lower, \c StrictlyLower, \c UnitLower.
* 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
*
@@ -187,11 +187,11 @@ struct ei_part_assignment_impl
}
else
{
ei_assert(Mode == Upper || Mode == Lower || Mode == StrictlyUpper || Mode == StrictlyLower);
if((Mode == Upper && row <= col)
|| (Mode == Lower && row >= col)
|| (Mode == StrictlyUpper && row < col)
|| (Mode == StrictlyLower && row > col))
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);
}
}
@@ -215,7 +215,7 @@ struct ei_part_assignment_impl<Derived1, Derived2, Mode, 0>
};
template<typename Derived1, typename Derived2>
struct ei_part_assignment_impl<Derived1, Derived2, Upper, Dynamic>
struct ei_part_assignment_impl<Derived1, Derived2, UpperTriangular, Dynamic>
{
inline static void run(Derived1 &dst, const Derived2 &src)
{
@@ -226,7 +226,7 @@ struct ei_part_assignment_impl<Derived1, Derived2, Upper, Dynamic>
};
template<typename Derived1, typename Derived2>
struct ei_part_assignment_impl<Derived1, Derived2, Lower, Dynamic>
struct ei_part_assignment_impl<Derived1, Derived2, LowerTriangular, Dynamic>
{
inline static void run(Derived1 &dst, const Derived2 &src)
{
@@ -237,7 +237,7 @@ struct ei_part_assignment_impl<Derived1, Derived2, Lower, Dynamic>
};
template<typename Derived1, typename Derived2>
struct ei_part_assignment_impl<Derived1, Derived2, StrictlyUpper, Dynamic>
struct ei_part_assignment_impl<Derived1, Derived2, StrictlyUpperTriangular, Dynamic>
{
inline static void run(Derived1 &dst, const Derived2 &src)
{
@@ -247,7 +247,7 @@ struct ei_part_assignment_impl<Derived1, Derived2, StrictlyUpper, Dynamic>
}
};
template<typename Derived1, typename Derived2>
struct ei_part_assignment_impl<Derived1, Derived2, StrictlyLower, Dynamic>
struct ei_part_assignment_impl<Derived1, Derived2, StrictlyLowerTriangular, Dynamic>
{
inline static void run(Derived1 &dst, const Derived2 &src)
{
@@ -285,8 +285,8 @@ void Part<MatrixType, Mode>::lazyAssign(const Other& other)
/** \returns a lvalue pseudo-expression allowing to perform special operations on \c *this.
*
* The \a Mode parameter can have the following values: \c Upper, \c StrictlyUpper, \c Lower,
* \c StrictlyLower, \c SelfAdjoint.
* 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
*
@@ -305,44 +305,44 @@ inline Part<Derived, Mode> MatrixBase<Derived>::part()
/** \returns true if *this is approximately equal to an upper triangular matrix,
* within the precision given by \a prec.
*
* \sa isLower(), extract(), part(), marked()
* \sa isLowerTriangular(), extract(), part(), marked()
*/
template<typename Derived>
bool MatrixBase<Derived>::isUpper(RealScalar prec) const
bool MatrixBase<Derived>::isUpperTriangular(RealScalar prec) const
{
if(cols() != rows()) return false;
RealScalar maxAbsOnUpperPart = static_cast<RealScalar>(-1);
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 > maxAbsOnUpperPart) maxAbsOnUpperPart = absValue;
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), maxAbsOnUpperPart, prec)) return false;
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 isUpper(), extract(), part(), marked()
* \sa isUpperTriangular(), extract(), part(), marked()
*/
template<typename Derived>
bool MatrixBase<Derived>::isLower(RealScalar prec) const
bool MatrixBase<Derived>::isLowerTriangular(RealScalar prec) const
{
if(cols() != rows()) return false;
RealScalar maxAbsOnLowerPart = static_cast<RealScalar>(-1);
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 > maxAbsOnLowerPart) maxAbsOnLowerPart = absValue;
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), maxAbsOnLowerPart, prec)) return false;
if(!ei_isMuchSmallerThan(coeff(i, j), maxAbsOnLowerTriangularPart, prec)) return false;
return true;
}

74
Eigen/src/Core/SolveTriangular.h
View File

@@ -30,9 +30,9 @@ template<typename XprType, unsigned int Mode> struct ei_is_part<Part<XprType,Mod
template<typename Lhs, typename Rhs,
int TriangularPart = (int(Lhs::Flags) & LowerTriangularBit)
? Lower
? LowerTriangular
: (int(Lhs::Flags) & UpperTriangularBit)
? Upper
? UpperTriangular
: -1,
int StorageOrder = ei_is_part<Lhs>::value ? -1 // this is to solve ambiguous specializations
: int(Lhs::Flags) & (RowMajorBit|SparseBit)
@@ -56,14 +56,14 @@ struct ei_solve_triangular_selector<Lhs,Rhs,UpLo,RowMajor|IsDense>
typedef typename Rhs::Scalar Scalar;
static void run(const Lhs& lhs, Rhs& other)
{
const bool IsLower = (UpLo==Lower);
const bool IsLowerTriangular = (UpLo==LowerTriangular);
const int size = lhs.cols();
/* We perform the inverse product per block of 4 rows such that we perfectly match
* our optimized matrix * vector product. blockyStart represents the number of rows
* we have process first using the non-block version.
*/
int blockyStart = (std::max(size-5,0)/4)*4;
if (IsLower)
if (IsLowerTriangular)
blockyStart = size - blockyStart;
else
blockyStart -= 1;
@@ -72,15 +72,15 @@ struct ei_solve_triangular_selector<Lhs,Rhs,UpLo,RowMajor|IsDense>
// process first rows using the non block version
if(!(Lhs::Flags & UnitDiagBit))
{
if (IsLower)
if (IsLowerTriangular)
other.coeffRef(0,c) = other.coeff(0,c)/lhs.coeff(0, 0);
else
other.coeffRef(size-1,c) = other.coeff(size-1, c)/lhs.coeff(size-1, size-1);
}
for(int i=(IsLower ? 1 : size-2); IsLower ? i<blockyStart : i>blockyStart; i += (IsLower ? 1 : -1) )
for(int i=(IsLowerTriangular ? 1 : size-2); IsLowerTriangular ? i<blockyStart : i>blockyStart; i += (IsLowerTriangular ? 1 : -1) )
{
Scalar tmp = other.coeff(i,c)
- (IsLower ? ((lhs.row(i).start(i)) * other.col(c).start(i)).coeff(0,0)
- (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)
other.coeffRef(i,c) = tmp;
@@ -89,15 +89,15 @@ struct ei_solve_triangular_selector<Lhs,Rhs,UpLo,RowMajor|IsDense>
}
// now let's process the remaining rows 4 at once
for(int i=blockyStart; IsLower ? i<size : i>0; )
for(int i=blockyStart; IsLowerTriangular ? i<size : i>0; )
{
int startBlock = i;
int endBlock = startBlock + (IsLower ? 4 : -4);
int endBlock = startBlock + (IsLowerTriangular ? 4 : -4);
/* Process the i cols times 4 rows block, and keep the result in a temporary vector */
// FIXME use fixed size block but take care to small fixed size matrices...
Matrix<Scalar,Dynamic,1> btmp(4);
if (IsLower)
if (IsLowerTriangular)
btmp = lhs.block(startBlock,0,4,i) * other.col(c).start(i);
else
btmp = lhs.block(i-3,i+1,4,size-1-i) * other.col(c).end(size-1-i);
@@ -106,21 +106,21 @@ struct ei_solve_triangular_selector<Lhs,Rhs,UpLo,RowMajor|IsDense>
* btmp stores the diagonal coefficients used to update the remaining part of the result.
*/
{
Scalar tmp = other.coeff(startBlock,c)-btmp.coeff(IsLower?0:3);
Scalar tmp = other.coeff(startBlock,c)-btmp.coeff(IsLowerTriangular?0:3);
if (Lhs::Flags & UnitDiagBit)
other.coeffRef(i,c) = tmp;
else
other.coeffRef(i,c) = tmp/lhs.coeff(i,i);
}
i += IsLower ? 1 : -1;
for (;IsLower ? i<endBlock : i>endBlock; i += IsLower ? 1 : -1)
i += IsLowerTriangular ? 1 : -1;
for (;IsLowerTriangular ? i<endBlock : i>endBlock; i += IsLowerTriangular ? 1 : -1)
{
int remainingSize = IsLower ? i-startBlock : startBlock-i;
int remainingSize = IsLowerTriangular ? i-startBlock : startBlock-i;
Scalar tmp = other.coeff(i,c)
- btmp.coeff(IsLower ? remainingSize : 3-remainingSize)
- ( lhs.row(i).segment(IsLower ? startBlock : i+1, remainingSize)
* other.col(c).segment(IsLower ? startBlock : i+1, remainingSize)).coeff(0,0);
- btmp.coeff(IsLowerTriangular ? remainingSize : 3-remainingSize)
- ( 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)
other.coeffRef(i,c) = tmp;
@@ -133,10 +133,10 @@ struct ei_solve_triangular_selector<Lhs,Rhs,UpLo,RowMajor|IsDense>
};
// Implements the following configurations:
// - inv(Lower, ColMajor) * Column vector
// - inv(Lower,UnitDiag,ColMajor) * Column vector
// - inv(Upper, ColMajor) * Column vector
// - inv(Upper,UnitDiag,ColMajor) * Column vector
// - inv(LowerTriangular, ColMajor) * Column vector
// - 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>
{
@@ -146,7 +146,7 @@ struct ei_solve_triangular_selector<Lhs,Rhs,UpLo,ColMajor|IsDense>
static void run(const Lhs& lhs, Rhs& other)
{
static const bool IsLower = (UpLo==Lower);
static const bool IsLowerTriangular = (UpLo==LowerTriangular);
const int size = lhs.cols();
for(int c=0 ; c<other.cols() ; ++c)
{
@@ -155,27 +155,27 @@ struct ei_solve_triangular_selector<Lhs,Rhs,UpLo,ColMajor|IsDense>
* we can process using the block version.
*/
int blockyEnd = (std::max(size-5,0)/4)*4;
if (!IsLower)
if (!IsLowerTriangular)
blockyEnd = size-1 - blockyEnd;
for(int i=IsLower ? 0 : size-1; IsLower ? i<blockyEnd : i>blockyEnd;)
for(int i=IsLowerTriangular ? 0 : size-1; IsLowerTriangular ? i<blockyEnd : i>blockyEnd;)
{
/* Let's process the 4x4 sub-matrix as usual.
* btmp stores the diagonal coefficients used to update the remaining part of the result.
*/
int startBlock = i;
int endBlock = startBlock + (IsLower ? 4 : -4);
int endBlock = startBlock + (IsLowerTriangular ? 4 : -4);
Matrix<Scalar,4,1> btmp;
for (;IsLower ? i<endBlock : i>endBlock;
i += IsLower ? 1 : -1)
for (;IsLowerTriangular ? i<endBlock : i>endBlock;
i += IsLowerTriangular ? 1 : -1)
{
if(!(Lhs::Flags & UnitDiagBit))
other.coeffRef(i,c) /= lhs.coeff(i,i);
int remainingSize = IsLower ? endBlock-i-1 : i-endBlock-1;
int remainingSize = IsLowerTriangular ? endBlock-i-1 : i-endBlock-1;
if (remainingSize>0)
other.col(c).segment((IsLower ? i : endBlock) + 1, remainingSize) -=
other.col(c).segment((IsLowerTriangular ? i : endBlock) + 1, remainingSize) -=
other.coeffRef(i,c)
* Block<Lhs,Dynamic,1>(lhs, (IsLower ? i : endBlock) + 1, i, remainingSize, 1);
btmp.coeffRef(IsLower ? i-startBlock : remainingSize) = -other.coeffRef(i,c);
* Block<Lhs,Dynamic,1>(lhs, (IsLowerTriangular ? i : endBlock) + 1, i, remainingSize, 1);
btmp.coeffRef(IsLowerTriangular ? i-startBlock : remainingSize) = -other.coeffRef(i,c);
}
/* Now we can efficiently update the remaining part of the result as a matrix * vector product.
@@ -187,11 +187,11 @@ struct ei_solve_triangular_selector<Lhs,Rhs,UpLo,ColMajor|IsDense>
// FIXME this is cool but what about conjugate/adjoint expressions ? do we want to evaluate them ?
// this is a more general problem though.
ei_cache_friendly_product_colmajor_times_vector(
IsLower ? size-endBlock : endBlock+1,
&(lhs.const_cast_derived().coeffRef(IsLower ? endBlock : 0, IsLower ? startBlock : endBlock+1)),
IsLowerTriangular ? size-endBlock : endBlock+1,
&(lhs.const_cast_derived().coeffRef(IsLowerTriangular ? endBlock : 0, IsLowerTriangular ? startBlock : endBlock+1)),
lhs.stride(),
btmp, &(other.coeffRef(IsLower ? endBlock : 0, c)));
// if (IsLower)
btmp, &(other.coeffRef(IsLowerTriangular ? endBlock : 0, c)));
// if (IsLowerTriangular)
// other.col(c).end(size-endBlock) += (lhs.block(endBlock, startBlock, size-endBlock, endBlock-startBlock)
// * other.col(c).block(startBlock,endBlock-startBlock)).lazy();
// else
@@ -201,7 +201,7 @@ struct ei_solve_triangular_selector<Lhs,Rhs,UpLo,ColMajor|IsDense>
/* Now we have to process the remaining part as usual */
int i;
for(i=blockyEnd; IsLower ? i<size-1 : i>0; i += (IsLower ? 1 : -1) )
for(i=blockyEnd; IsLowerTriangular ? i<size-1 : i>0; i += (IsLowerTriangular ? 1 : -1) )
{
if(!(Lhs::Flags & UnitDiagBit))
other.coeffRef(i,c) /= lhs.coeff(i,i);
@@ -209,7 +209,7 @@ struct ei_solve_triangular_selector<Lhs,Rhs,UpLo,ColMajor|IsDense>
/* NOTE we cannot use lhs.col(i).end(size-i-1) because Part::coeffRef gets called by .col() to
* get the address of the start of the row
*/
if(IsLower)
if(IsLowerTriangular)
other.col(c).end(size-i-1) -= other.coeffRef(i,c) * Block<Lhs,Dynamic,1>(lhs, i+1,i, size-i-1,1);
else
other.col(c).start(i) -= other.coeffRef(i,c) * Block<Lhs,Dynamic,1>(lhs, 0,i, i, 1);

4
Eigen/src/Core/Transpose.h
View File

@@ -167,7 +167,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<StrictlyUpper>().swap(m.transpose());
m.template part<StrictlyUpperTriangular>().swap(m.transpose());
}
};
@@ -175,7 +175,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<StrictlyUpper>().swap(m.transpose());
m.template part<StrictlyUpperTriangular>().swap(m.transpose());
else
m.set(m.transpose().eval());
}

14
Eigen/src/Core/util/Constants.h
View File

@@ -185,16 +185,16 @@ const unsigned int HereditaryBits = RowMajorBit
| SparseBit;
// Possible values for the Mode parameter of part() and of extract()
const unsigned int Upper = UpperTriangularBit;
const unsigned int StrictlyUpper = UpperTriangularBit | ZeroDiagBit;
const unsigned int Lower = LowerTriangularBit;
const unsigned int StrictlyLower = LowerTriangularBit | ZeroDiagBit;
const unsigned int UpperTriangular = UpperTriangularBit;
const unsigned int StrictlyUpperTriangular = UpperTriangularBit | ZeroDiagBit;
const unsigned int LowerTriangular = LowerTriangularBit;
const unsigned int StrictlyLowerTriangular = LowerTriangularBit | ZeroDiagBit;
const unsigned int SelfAdjoint = SelfAdjointBit;
// additional possible values for the Mode parameter of extract()
const unsigned int UnitUpper = UpperTriangularBit | UnitDiagBit;
const unsigned int UnitLower = LowerTriangularBit | UnitDiagBit;
const unsigned int Diagonal = Upper | Lower;
const unsigned int UnitUpperTriangular = UpperTriangularBit | UnitDiagBit;
const unsigned int UnitLowerTriangular = LowerTriangularBit | UnitDiagBit;
const unsigned int Diagonal = UpperTriangular | LowerTriangular;
enum { Aligned, Unaligned };
enum { ForceAligned, AsRequested };