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

big huge changes, so i dont remember everything.

Browse Source 
* renaming, e.g. LU ---> FullPivLU
* split tests framework: more robust, e.g. dont generate empty tests if a number is skipped
* make all remaining tests use that splitting, as needed.
* Fix 4x4 inversion (see stable branch)
* Transform::inverse() and geo_transform test : adapt to new inverse() API, it was also trying to instantiate inverse() for 3x4 matrices.
* CMakeLists: more robust regexp to parse the version number
* misc fixes in unit tests
This commit is contained in:
Benoit Jacob
2009-10-28 18:19:29 -04:00
parent 1f1c04cac1
commit 2840ac7e94
99 changed files with 816 additions and 710 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
Eigen/src/LU/Determinant.h
View File

@@ -53,7 +53,7 @@ template<typename Derived,
{
static inline typename ei_traits<Derived>::Scalar run(const Derived& m)
{
return m.partialLu().determinant();
return m.partialPivLu().determinant();
}
};

62
Eigen/src/LU/LU.h → Eigen/src/LU/FullPivLU.h
View File

@@ -31,7 +31,7 @@ template<typename MatrixType> struct ei_lu_image_impl;
/** \ingroup LU_Module
*
* \class LU
* \class FullPivLU
*
* \brief LU decomposition of a matrix with complete pivoting, and related features
*
@@ -54,12 +54,12 @@ template<typename MatrixType> struct ei_lu_image_impl;
* permutationP(), permutationQ().
*
* As an exemple, here is how the original matrix can be retrieved:
* \include class_LU.cpp
* Output: \verbinclude class_LU.out
* \include class_FullPivLU.cpp
* Output: \verbinclude class_FullPivLU.out
*
* \sa MatrixBase::lu(), MatrixBase::determinant(), MatrixBase::inverse()
* \sa MatrixBase::fullPivLu(), MatrixBase::determinant(), MatrixBase::inverse()
*/
template<typename MatrixType> class LU
template<typename MatrixType> class FullPivLU
{
public:
@@ -81,14 +81,14 @@ template<typename MatrixType> class LU
* The default constructor is useful in cases in which the user intends to
* perform decompositions via LU::compute(const MatrixType&).
*/
LU();
FullPivLU();
/** Constructor.
*
* \param matrix the matrix of which to compute the LU decomposition.
* It is required to be nonzero.
*/
LU(const MatrixType& matrix);
FullPivLU(const MatrixType& matrix);
/** Computes the LU decomposition of the given matrix.
*
@@ -97,11 +97,11 @@ template<typename MatrixType> class LU
*
* \returns a reference to *this
*/
LU& compute(const MatrixType& matrix);
FullPivLU& compute(const MatrixType& matrix);
/** \returns the LU decomposition matrix: the upper-triangular part is U, the
* unit-lower-triangular part is L (at least for square matrices; in the non-square
* case, special care is needed, see the documentation of class LU).
* case, special care is needed, see the documentation of class FullPivLU).
*
* \sa matrixL(), matrixU()
*/
@@ -131,7 +131,7 @@ template<typename MatrixType> class LU
/** \returns a vector of integers, whose size is the number of rows of the matrix being decomposed,
* representing the P permutation i.e. the permutation of the rows. For its precise meaning,
* see the examples given in the documentation of class LU.
* see the examples given in the documentation of class FullPivLU.
*
* \sa permutationQ()
*/
@@ -143,7 +143,7 @@ template<typename MatrixType> class LU
/** \returns a vector of integers, whose size is the number of columns of the matrix being
* decomposed, representing the Q permutation i.e. the permutation of the columns.
* For its precise meaning, see the examples given in the documentation of class LU.
* For its precise meaning, see the examples given in the documentation of class FullPivLU.
*
* \sa permutationP()
*/
@@ -162,8 +162,8 @@ template<typename MatrixType> class LU
* For that, it uses the threshold value that you can control by calling
* setThreshold(const RealScalar&).
*
* Example: \include LU_kernel.cpp
* Output: \verbinclude LU_kernel.out
* Example: \include FullPivLU_kernel.cpp
* Output: \verbinclude FullPivLU_kernel.out
*
* \sa image()
*/
@@ -187,8 +187,8 @@ template<typename MatrixType> class LU
* For that, it uses the threshold value that you can control by calling
* setThreshold(const RealScalar&).
*
* Example: \include LU_image.cpp
* Output: \verbinclude LU_image.out
* Example: \include FullPivLU_image.cpp
* Output: \verbinclude FullPivLU_image.out
*
* \sa kernel()
*/
@@ -214,8 +214,8 @@ template<typename MatrixType> class LU
* \note_about_arbitrary_choice_of_solution
* \note_about_using_kernel_to_study_multiple_solutions
*
* Example: \include LU_solve.cpp
* Output: \verbinclude LU_solve.out
* Example: \include FullPivLU_solve.cpp
* Output: \verbinclude FullPivLU_solve.out
*
* \sa TriangularView::solve(), kernel(), inverse()
*/
@@ -260,7 +260,7 @@ template<typename MatrixType> class LU
*
* If you want to come back to the default behavior, call setThreshold(Default_t)
*/
LU& setThreshold(const RealScalar& threshold)
FullPivLU& setThreshold(const RealScalar& threshold)
{
m_usePrescribedThreshold = true;
m_prescribedThreshold = threshold;
@@ -274,7 +274,7 @@ template<typename MatrixType> class LU
*
* See the documentation of setThreshold(const RealScalar&).
*/
LU& setThreshold(Default_t)
FullPivLU& setThreshold(Default_t)
{
m_usePrescribedThreshold = false;
}
@@ -383,20 +383,20 @@ template<typename MatrixType> class LU
};
template<typename MatrixType>
LU<MatrixType>::LU()
FullPivLU<MatrixType>::FullPivLU()
: m_isInitialized(false), m_usePrescribedThreshold(false)
{
}
template<typename MatrixType>
LU<MatrixType>::LU(const MatrixType& matrix)
FullPivLU<MatrixType>::FullPivLU(const MatrixType& matrix)
: m_isInitialized(false), m_usePrescribedThreshold(false)
{
compute(matrix);
}
template<typename MatrixType>
LU<MatrixType>& LU<MatrixType>::compute(const MatrixType& matrix)
FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const MatrixType& matrix)
{
m_isInitialized = true;
m_lu = matrix;
@@ -483,7 +483,7 @@ LU<MatrixType>& LU<MatrixType>::compute(const MatrixType& matrix)
}
template<typename MatrixType>
typename ei_traits<MatrixType>::Scalar LU<MatrixType>::determinant() const
typename ei_traits<MatrixType>::Scalar FullPivLU<MatrixType>::determinant() const
{
ei_assert(m_isInitialized && "LU is not initialized.");
ei_assert(m_lu.rows() == m_lu.cols() && "You can't take the determinant of a non-square matrix!");
@@ -511,7 +511,7 @@ struct ei_traits<ei_lu_kernel_impl<MatrixType> >
template<typename MatrixType>
struct ei_lu_kernel_impl : public ReturnByValue<ei_lu_kernel_impl<MatrixType> >
{
typedef LU<MatrixType> LUType;
typedef FullPivLU<MatrixType> LUType;
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
const LUType& m_lu;
@@ -615,7 +615,7 @@ struct ei_traits<ei_lu_image_impl<MatrixType> >
template<typename MatrixType>
struct ei_lu_image_impl : public ReturnByValue<ei_lu_image_impl<MatrixType> >
{
typedef LU<MatrixType> LUType;
typedef FullPivLU<MatrixType> LUType;
typedef typename MatrixType::RealScalar RealScalar;
const LUType& m_lu;
int m_rank, m_cols;
@@ -670,7 +670,7 @@ template<typename MatrixType, typename Rhs>
struct ei_lu_solve_impl : public ReturnByValue<ei_lu_solve_impl<MatrixType, Rhs> >
{
typedef typename ei_cleantype<typename Rhs::Nested>::type RhsNested;
typedef LU<MatrixType> LUType;
typedef FullPivLU<MatrixType> LUType;
const LUType& m_lu;
const typename Rhs::Nested m_rhs;
@@ -739,15 +739,15 @@ struct ei_lu_solve_impl : public ReturnByValue<ei_lu_solve_impl<MatrixType, Rhs>
/** \lu_module
*
* \return the LU decomposition of \c *this.
* \return the full-pivoting LU decomposition of \c *this.
*
* \sa class LU
* \sa class FullPivLU
*/
template<typename Derived>
inline const LU<typename MatrixBase<Derived>::PlainMatrixType>
MatrixBase<Derived>::lu() const
inline const FullPivLU<typename MatrixBase<Derived>::PlainMatrixType>
MatrixBase<Derived>::fullPivLu() const
{
return LU<PlainMatrixType>(eval());
return FullPivLU<PlainMatrixType>(eval());
}
#endif // EIGEN_LU_H

47
Eigen/src/LU/Inverse.h
View File

@@ -34,7 +34,7 @@ struct ei_compute_inverse
{
static inline void run(const MatrixType& matrix, ResultType& result)
{
result = matrix.partialLu().inverse();
result = matrix.partialPivLu().inverse();
}
};
@@ -232,22 +232,31 @@ struct ei_compute_inverse<MatrixType, ResultType, 4>
typename MatrixType::PlainMatrixType matrix(_matrix);
// let's extract from the 2 first colums a 2x2 block whose determinant is as big as possible.
int good_row0=0, good_row1=1;
RealScalar good_absdet(-1);
// this double for loop shouldn't be too costly: only 6 iterations
for(int row0=0; row0<4; ++row0) {
for(int row1=row0+1; row1<4; ++row1)
{
RealScalar absdet = ei_abs(matrix.coeff(row0,0)*matrix.coeff(row1,1)
- matrix.coeff(row0,1)*matrix.coeff(row1,0));
if(absdet > good_absdet)
{
good_absdet = absdet;
good_row0 = row0;
good_row1 = row1;
}
}
}
int good_row0, good_row1, good_i;
Matrix<RealScalar,6,1> absdet;
// any 2x2 block with determinant above this threshold will be considered good enough
RealScalar d = (matrix.col(0).squaredNorm()+matrix.col(1).squaredNorm()) * RealScalar(1e-2);
#define ei_inv_size4_helper_macro(i,row0,row1) \
absdet[i] = ei_abs(matrix.coeff(row0,0)*matrix.coeff(row1,1) \
- matrix.coeff(row0,1)*matrix.coeff(row1,0)); \
if(absdet[i] > d) { good_row0=row0; good_row1=row1; goto good; }
ei_inv_size4_helper_macro(0,0,1)
ei_inv_size4_helper_macro(1,0,2)
ei_inv_size4_helper_macro(2,0,3)
ei_inv_size4_helper_macro(3,1,2)
ei_inv_size4_helper_macro(4,1,3)
ei_inv_size4_helper_macro(5,2,3)
// no 2x2 block has determinant bigger than the threshold. So just take the one that
// has the biggest determinant
absdet.maxCoeff(&good_i);
good_row0 = good_i <= 2 ? 0 : good_i <= 4 ? 1 : 2;
good_row1 = good_i <= 2 ? good_i+1 : good_i <= 4 ? good_i-1 : 3;
// now good_row0 and good_row1 are correctly set
good:
// do row permutations to move this 2x2 block to the top
matrix.row(0).swap(matrix.row(good_row0));
matrix.row(1).swap(matrix.row(good_row1));
@@ -318,12 +327,12 @@ struct ei_inverse_impl : public ReturnByValue<ei_inverse_impl<MatrixType> >
* \returns the matrix inverse of this matrix.
*
* For small fixed sizes up to 4x4, this method uses ad-hoc methods (cofactors up to 3x3, Euler's trick for 4x4).
* In the general case, this method uses class PartialLU.
* In the general case, this method uses class PartialPivLU.
*
* \note This matrix must be invertible, otherwise the result is undefined. If you need an
* invertibility check, do the following:
* \li for fixed sizes up to 4x4, use computeInverseAndDetWithCheck().
* \li for the general case, use class LU.
* \li for the general case, use class FullPivLU.
*
* Example: \include MatrixBase_inverse.cpp
* Output: \verbinclude MatrixBase_inverse.out

75
Eigen/src/LU/PartialLU.h → Eigen/src/LU/PartialPivLU.h
View File

@@ -30,7 +30,7 @@ template<typename MatrixType, typename Rhs> struct ei_partiallu_solve_impl;
/** \ingroup LU_Module
*
* \class PartialLU
* \class PartialPivLU
*
* \brief LU decomposition of a matrix with partial pivoting, and related features
*
@@ -46,10 +46,10 @@ template<typename MatrixType, typename Rhs> struct ei_partiallu_solve_impl;
* matrix is invertible: it is your task to check that you only use this decomposition on invertible matrices.
*
* The guaranteed safe alternative, working for all matrices, is the full pivoting LU decomposition, provided
* by class LU.
* by class FullPivLU.
*
* This is \b not a rank-revealing LU decomposition. Many features are intentionally absent from this class,
* such as rank computation. If you need these features, use class LU.
* such as rank computation. If you need these features, use class FullPivLU.
*
* This LU decomposition is suitable to invert invertible matrices. It is what MatrixBase::inverse() uses
* in the general case.
@@ -57,9 +57,9 @@ template<typename MatrixType, typename Rhs> struct ei_partiallu_solve_impl;
*
* The data of the LU decomposition can be directly accessed through the methods matrixLU(), permutationP().
*
* \sa MatrixBase::partialLu(), MatrixBase::determinant(), MatrixBase::inverse(), MatrixBase::computeInverse(), class LU
* \sa MatrixBase::partialPivLu(), MatrixBase::determinant(), MatrixBase::inverse(), MatrixBase::computeInverse(), class FullPivLU
*/
template<typename MatrixType> class PartialLU
template<typename MatrixType> class PartialPivLU
{
public:
@@ -79,40 +79,40 @@ template<typename MatrixType> class PartialLU
* \brief Default Constructor.
*
* The default constructor is useful in cases in which the user intends to
* perform decompositions via PartialLU::compute(const MatrixType&).
* perform decompositions via PartialPivLU::compute(const MatrixType&).
*/
PartialLU();
PartialPivLU();
/** Constructor.
*
* \param matrix the matrix of which to compute the LU decomposition.
*
* \warning The matrix should have full rank (e.g. if it's square, it should be invertible).
* If you need to deal with non-full rank, use class LU instead.
* If you need to deal with non-full rank, use class FullPivLU instead.
*/
PartialLU(const MatrixType& matrix);
PartialPivLU(const MatrixType& matrix);
PartialLU& compute(const MatrixType& matrix);
PartialPivLU& compute(const MatrixType& matrix);
/** \returns the LU decomposition matrix: the upper-triangular part is U, the
* unit-lower-triangular part is L (at least for square matrices; in the non-square
* case, special care is needed, see the documentation of class LU).
* case, special care is needed, see the documentation of class FullPivLU).
*
* \sa matrixL(), matrixU()
*/
inline const MatrixType& matrixLU() const
{
ei_assert(m_isInitialized && "PartialLU is not initialized.");
ei_assert(m_isInitialized && "PartialPivLU is not initialized.");
return m_lu;
}
/** \returns a vector of integers, whose size is the number of rows of the matrix being decomposed,
* representing the P permutation i.e. the permutation of the rows. For its precise meaning,
* see the examples given in the documentation of class LU.
* see the examples given in the documentation of class FullPivLU.
*/
inline const IntColVectorType& permutationP() const
{
ei_assert(m_isInitialized && "PartialLU is not initialized.");
ei_assert(m_isInitialized && "PartialPivLU is not initialized.");
return m_p;
}
@@ -125,10 +125,10 @@ template<typename MatrixType> class PartialLU
*
* \returns the solution.
*
* Example: \include PartialLU_solve.cpp
* Output: \verbinclude PartialLU_solve.out
* Example: \include PartialPivLU_solve.cpp
* Output: \verbinclude PartialPivLU_solve.out
*
* Since this PartialLU class assumes anyway that the matrix A is invertible, the solution
* Since this PartialPivLU class assumes anyway that the matrix A is invertible, the solution
* theoretically exists and is unique regardless of b.
*
* \note_about_checking_solutions
@@ -146,7 +146,7 @@ template<typename MatrixType> class PartialLU
/** \returns the inverse of the matrix of which *this is the LU decomposition.
*
* \warning The matrix being decomposed here is assumed to be invertible. If you need to check for
* invertibility, use class LU instead.
* invertibility, use class FullPivLU instead.
*
* \sa MatrixBase::inverse(), LU::inverse()
*/
@@ -180,7 +180,7 @@ template<typename MatrixType> class PartialLU
};
template<typename MatrixType>
PartialLU<MatrixType>::PartialLU()
PartialPivLU<MatrixType>::PartialPivLU()
: m_lu(),
m_p(),
m_det_p(0),
@@ -189,7 +189,7 @@ PartialLU<MatrixType>::PartialLU()
}
template<typename MatrixType>
PartialLU<MatrixType>::PartialLU(const MatrixType& matrix)
PartialPivLU<MatrixType>::PartialPivLU(const MatrixType& matrix)
: m_lu(),
m_p(),
m_det_p(0),
@@ -378,12 +378,12 @@ void ei_partial_lu_inplace(MatrixType& lu, IntVector& row_transpositions, int& n
}
template<typename MatrixType>
PartialLU<MatrixType>& PartialLU<MatrixType>::compute(const MatrixType& matrix)
PartialPivLU<MatrixType>& PartialPivLU<MatrixType>::compute(const MatrixType& matrix)
{
m_lu = matrix;
m_p.resize(matrix.rows());
ei_assert(matrix.rows() == matrix.cols() && "PartialLU is only for square (and moreover invertible) matrices");
ei_assert(matrix.rows() == matrix.cols() && "PartialPivLU is only for square (and moreover invertible) matrices");
const int size = matrix.rows();
IntColVectorType rows_transpositions(size);
@@ -401,9 +401,9 @@ PartialLU<MatrixType>& PartialLU<MatrixType>::compute(const MatrixType& matrix)
}
template<typename MatrixType>
typename ei_traits<MatrixType>::Scalar PartialLU<MatrixType>::determinant() const
typename ei_traits<MatrixType>::Scalar PartialPivLU<MatrixType>::determinant() const
{
ei_assert(m_isInitialized && "PartialLU is not initialized.");
ei_assert(m_isInitialized && "PartialPivLU is not initialized.");
return Scalar(m_det_p) * m_lu.diagonal().prod();
}
@@ -424,7 +424,7 @@ template<typename MatrixType, typename Rhs>
struct ei_partiallu_solve_impl : public ReturnByValue<ei_partiallu_solve_impl<MatrixType, Rhs> >
{
typedef typename ei_cleantype<typename Rhs::Nested>::type RhsNested;
typedef PartialLU<MatrixType> LUType;
typedef PartialPivLU<MatrixType> LUType;
const LUType& m_lu;
const typename Rhs::Nested m_rhs;
@@ -464,15 +464,30 @@ struct ei_partiallu_solve_impl : public ReturnByValue<ei_partiallu_solve_impl<Ma
/** \lu_module
*
* \return the LU decomposition of \c *this.
* \return the partial-pivoting LU decomposition of \c *this.
*
* \sa class LU
* \sa class PartialPivLU
*/
template<typename Derived>
inline const PartialLU<typename MatrixBase<Derived>::PlainMatrixType>
MatrixBase<Derived>::partialLu() const
inline const PartialPivLU<typename MatrixBase<Derived>::PlainMatrixType>
MatrixBase<Derived>::partialPivLu() const
{
return PartialLU<PlainMatrixType>(eval());
return PartialPivLU<PlainMatrixType>(eval());
}
/** \lu_module
*
* Synonym of partialPivLu().
*
* \return the partial-pivoting LU decomposition of \c *this.
*
* \sa class PartialPivLU
*/
template<typename Derived>
inline const PartialPivLU<typename MatrixBase<Derived>::PlainMatrixType>
MatrixBase<Derived>::lu() const
{
return PartialPivLU<PlainMatrixType>(eval());
}
#endif // EIGEN_PARTIALLU_H