GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-06-04 18:54:00 +08:00
Code Issues Packages Projects Releases Wiki Activity

remove the m_originalMatrix member. Instead, image() now takes the original matrix as parameter. It was the only method to use it anyway. Introduce m_isInitialized.

Browse Source
This commit is contained in:
Benoit Jacob 2009-10-18 15:21:19 -04:00
parent d71c7f42d3
commit 47eeb40380
2 changed files with 37 additions and 37 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

68
Eigen/src/LU/LU.h
View File

@ -107,7 +107,7 @@ template<typename MatrixType> class LU
*/ */
inline const MatrixType& matrixLU() const inline const MatrixType& matrixLU() const
{ {
ei_assert(m_originalMatrix != 0 && "LU is not initialized."); ei_assert(m_isInitialized && "LU is not initialized.");
return m_lu; return m_lu;
} }
@ -120,7 +120,7 @@ template<typename MatrixType> class LU
*/ */
inline int nonzeroPivots() const inline int nonzeroPivots() const
{ {
ei_assert(m_originalMatrix != 0 && "LU is not initialized."); ei_assert(m_isInitialized && "LU is not initialized.");
return m_nonzero_pivots; return m_nonzero_pivots;
} }
@ -129,14 +129,6 @@ template<typename MatrixType> class LU
*/ */
RealScalar maxPivot() const { return m_maxpivot; } RealScalar maxPivot() const { return m_maxpivot; }
/** \returns a pointer to the matrix of which *this is the LU decomposition.
*/
inline const MatrixType* originalMatrix() const
{
ei_assert(m_originalMatrix != 0 && "LU is not initialized.");
return m_originalMatrix;
}
/** \returns a vector of integers, whose size is the number of rows of the matrix being decomposed, /** \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, * 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 LU.
@ -145,7 +137,7 @@ template<typename MatrixType> class LU
*/ */
inline const IntColVectorType& permutationP() const inline const IntColVectorType& permutationP() const
{ {
ei_assert(m_originalMatrix != 0 && "LU is not initialized."); ei_assert(m_isInitialized && "LU is not initialized.");
return m_p; return m_p;
} }
@ -157,7 +149,7 @@ template<typename MatrixType> class LU
*/ */
inline const IntRowVectorType& permutationQ() const inline const IntRowVectorType& permutationQ() const
{ {
ei_assert(m_originalMatrix != 0 && "LU is not initialized."); ei_assert(m_isInitialized && "LU is not initialized.");
return m_q; return m_q;
} }
@ -177,13 +169,18 @@ template<typename MatrixType> class LU
*/ */
inline const ei_lu_kernel_impl<MatrixType> kernel() const inline const ei_lu_kernel_impl<MatrixType> kernel() const
{ {
ei_assert(m_originalMatrix != 0 && "LU is not initialized."); ei_assert(m_isInitialized && "LU is not initialized.");
return ei_lu_kernel_impl<MatrixType>(*this); return ei_lu_kernel_impl<MatrixType>(*this);
} }
/** \returns the image of the matrix, also called its column-space. The columns of the returned matrix /** \returns the image of the matrix, also called its column-space. The columns of the returned matrix
* will form a basis of the kernel. * will form a basis of the kernel.
* *
* \param originalMatrix the original matrix, of which *this is the LU decomposition.
* The reason why it is needed to pass it here, is that this allows
* a large optimization, as otherwise this method would need to reconstruct it
* from the LU decomposition.
*
* \note If the image has dimension zero, then the returned matrix is a column-vector filled with zeros. * \note If the image has dimension zero, then the returned matrix is a column-vector filled with zeros.
* *
* \note This method has to determine which pivots should be considered nonzero. * \note This method has to determine which pivots should be considered nonzero.
@ -195,10 +192,12 @@ template<typename MatrixType> class LU
* *
* \sa kernel() * \sa kernel()
*/ */
inline const ei_lu_image_impl<MatrixType> image() const template<typename OriginalMatrixType>
inline const ei_lu_image_impl<MatrixType>
image(const MatrixBase<OriginalMatrixType>& originalMatrix) const
{ {
ei_assert(m_originalMatrix != 0 && "LU is not initialized."); ei_assert(m_isInitialized && "LU is not initialized.");
return ei_lu_image_impl<MatrixType>(*this); return ei_lu_image_impl<MatrixType>(*this, originalMatrix.derived());
} }
/** This method returns a solution x to the equation Ax=b, where A is the matrix of which /** This method returns a solution x to the equation Ax=b, where A is the matrix of which
@ -224,7 +223,7 @@ template<typename MatrixType> class LU
inline const ei_lu_solve_impl<MatrixType, Rhs> inline const ei_lu_solve_impl<MatrixType, Rhs>
solve(const MatrixBase<Rhs>& b) const solve(const MatrixBase<Rhs>& b) const
{ {
ei_assert(m_originalMatrix != 0 && "LU is not initialized."); ei_assert(m_isInitialized && "LU is not initialized.");
return ei_lu_solve_impl<MatrixType, Rhs>(*this, b.derived()); return ei_lu_solve_impl<MatrixType, Rhs>(*this, b.derived());
} }
@ -286,7 +285,7 @@ template<typename MatrixType> class LU
*/ */
RealScalar threshold() const RealScalar threshold() const
{ {
ei_assert(m_originalMatrix != 0 || m_usePrescribedThreshold); ei_assert(m_isInitialized || m_usePrescribedThreshold);
return m_usePrescribedThreshold ? m_prescribedThreshold return m_usePrescribedThreshold ? m_prescribedThreshold
// this formula comes from experimenting (see "LU precision tuning" thread on the list) // this formula comes from experimenting (see "LU precision tuning" thread on the list)
// and turns out to be identical to Higham's formula used already in LDLt. // and turns out to be identical to Higham's formula used already in LDLt.
@ -301,7 +300,7 @@ template<typename MatrixType> class LU
*/ */
inline int rank() const inline int rank() const
{ {
ei_assert(m_originalMatrix != 0 && "LU is not initialized."); ei_assert(m_isInitialized && "LU is not initialized.");
RealScalar premultiplied_threshold = ei_abs(m_maxpivot) * threshold(); RealScalar premultiplied_threshold = ei_abs(m_maxpivot) * threshold();
int result = 0; int result = 0;
for(int i = 0; i < m_nonzero_pivots; ++i) for(int i = 0; i < m_nonzero_pivots; ++i)
@ -317,7 +316,7 @@ template<typename MatrixType> class LU
*/ */
inline int dimensionOfKernel() const inline int dimensionOfKernel() const
{ {
ei_assert(m_originalMatrix != 0 && "LU is not initialized."); ei_assert(m_isInitialized && "LU is not initialized.");
return m_lu.cols() - rank(); return m_lu.cols() - rank();
} }
@ -330,7 +329,7 @@ template<typename MatrixType> class LU
*/ */
inline bool isInjective() const inline bool isInjective() const
{ {
ei_assert(m_originalMatrix != 0 && "LU is not initialized."); ei_assert(m_isInitialized && "LU is not initialized.");
return rank() == m_lu.cols(); return rank() == m_lu.cols();
} }
@ -343,7 +342,7 @@ template<typename MatrixType> class LU
*/ */
inline bool isSurjective() const inline bool isSurjective() const
{ {
ei_assert(m_originalMatrix != 0 && "LU is not initialized."); ei_assert(m_isInitialized && "LU is not initialized.");
return rank() == m_lu.rows(); return rank() == m_lu.rows();
} }
@ -355,7 +354,7 @@ template<typename MatrixType> class LU
*/ */
inline bool isInvertible() const inline bool isInvertible() const
{ {
ei_assert(m_originalMatrix != 0 && "LU is not initialized."); ei_assert(m_isInitialized && "LU is not initialized.");
return isInjective() && (m_lu.rows() == m_lu.cols()); return isInjective() && (m_lu.rows() == m_lu.cols());
} }
@ -368,14 +367,14 @@ template<typename MatrixType> class LU
*/ */
inline const ei_lu_solve_impl<MatrixType,NestByValue<typename MatrixType::IdentityReturnType> > inverse() const inline const ei_lu_solve_impl<MatrixType,NestByValue<typename MatrixType::IdentityReturnType> > inverse() const
{ {
ei_assert(m_originalMatrix != 0 && "LU is not initialized."); ei_assert(m_isInitialized && "LU is not initialized.");
ei_assert(m_lu.rows() == m_lu.cols() && "You can't take the inverse of a non-square matrix!"); ei_assert(m_lu.rows() == m_lu.cols() && "You can't take the inverse of a non-square matrix!");
return ei_lu_solve_impl<MatrixType,NestByValue<typename MatrixType::IdentityReturnType> > return ei_lu_solve_impl<MatrixType,NestByValue<typename MatrixType::IdentityReturnType> >
(*this, MatrixType::Identity(m_lu.rows(), m_lu.cols()).nestByValue()); (*this, MatrixType::Identity(m_lu.rows(), m_lu.cols()).nestByValue());
} }
protected: protected:
const MatrixType* m_originalMatrix; bool m_isInitialized;
MatrixType m_lu; MatrixType m_lu;
IntColVectorType m_p; IntColVectorType m_p;
IntRowVectorType m_q; IntRowVectorType m_q;
@ -388,13 +387,13 @@ template<typename MatrixType> class LU
template<typename MatrixType> template<typename MatrixType>
LU<MatrixType>::LU() LU<MatrixType>::LU()
: m_originalMatrix(0), m_usePrescribedThreshold(false) : m_isInitialized(false), m_usePrescribedThreshold(false)
{ {
} }
template<typename MatrixType> template<typename MatrixType>
LU<MatrixType>::LU(const MatrixType& matrix) LU<MatrixType>::LU(const MatrixType& matrix)
: m_originalMatrix(0), m_usePrescribedThreshold(false) : m_isInitialized(false), m_usePrescribedThreshold(false)
{ {
compute(matrix); compute(matrix);
} }
@ -402,7 +401,7 @@ LU<MatrixType>::LU(const MatrixType& matrix)
template<typename MatrixType> template<typename MatrixType>
LU<MatrixType>& LU<MatrixType>::compute(const MatrixType& matrix) LU<MatrixType>& LU<MatrixType>::compute(const MatrixType& matrix)
{ {
m_originalMatrix = &matrix; m_isInitialized = true;
m_lu = matrix; m_lu = matrix;
m_p.resize(matrix.rows()); m_p.resize(matrix.rows());
m_q.resize(matrix.cols()); m_q.resize(matrix.cols());
@ -475,7 +474,7 @@ LU<MatrixType>& LU<MatrixType>::compute(const MatrixType& matrix)
template<typename MatrixType> template<typename MatrixType>
typename ei_traits<MatrixType>::Scalar LU<MatrixType>::determinant() const typename ei_traits<MatrixType>::Scalar LU<MatrixType>::determinant() const
{ {
ei_assert(m_originalMatrix != 0 && "LU is not initialized."); 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!"); ei_assert(m_lu.rows() == m_lu.cols() && "You can't take the determinant of a non-square matrix!");
return Scalar(m_det_pq) * m_lu.diagonal().prod(); return Scalar(m_det_pq) * m_lu.diagonal().prod();
} }
@ -605,9 +604,10 @@ struct ei_lu_image_impl : public ReturnByValue<ei_lu_image_impl<MatrixType> >
typedef typename MatrixType::RealScalar RealScalar; typedef typename MatrixType::RealScalar RealScalar;
const LUType& m_lu; const LUType& m_lu;
int m_rank; int m_rank;
const MatrixType& m_originalMatrix;
ei_lu_image_impl(const LUType& lu) ei_lu_image_impl(const LUType& lu, const MatrixType& originalMatrix)
: m_lu(lu), m_rank(lu.rank()) {} : m_lu(lu), m_rank(lu.rank()), m_originalMatrix(originalMatrix) {}
inline int rows() const { return m_lu.matrixLU().cols(); } inline int rows() const { return m_lu.matrixLU().cols(); }
inline int cols() const { return m_rank; } inline int cols() const { return m_rank; }
@ -619,7 +619,7 @@ struct ei_lu_image_impl : public ReturnByValue<ei_lu_image_impl<MatrixType> >
// The Image is just {0}, so it doesn't have a basis properly speaking, but let's // The Image is just {0}, so it doesn't have a basis properly speaking, but let's
// avoid crashing/asserting as that depends on floating point calculations. Let's // avoid crashing/asserting as that depends on floating point calculations. Let's
// just return a single column vector filled with zeros. // just return a single column vector filled with zeros.
dst.resize(m_lu.originalMatrix()->rows(), 1); dst.resize(m_originalMatrix.rows(), 1);
dst.setZero(); dst.setZero();
return; return;
} }
@ -632,9 +632,9 @@ struct ei_lu_image_impl : public ReturnByValue<ei_lu_image_impl<MatrixType> >
pivots.coeffRef(p++) = i; pivots.coeffRef(p++) = i;
ei_assert(p == m_rank && "You hit a bug in Eigen! Please report (backtrace and matrix)!"); ei_assert(p == m_rank && "You hit a bug in Eigen! Please report (backtrace and matrix)!");
dst.resize(m_lu.originalMatrix()->rows(), m_rank); dst.resize(m_originalMatrix.rows(), m_rank);
for(int i = 0; i < m_rank; ++i) for(int i = 0; i < m_rank; ++i)
dst.col(i) = m_lu.originalMatrix()->col(m_lu.permutationQ().coeff(pivots.coeff(i))); dst.col(i) = m_originalMatrix.col(m_lu.permutationQ().coeff(pivots.coeff(i)));
} }
}; };

6
test/lu.cpp
View File

@ -38,7 +38,7 @@ template<typename MatrixType> void lu_non_invertible()
LU<MatrixType> lu(m1); LU<MatrixType> lu(m1);
typename ei_lu_kernel_impl<MatrixType>::ReturnMatrixType m1kernel = lu.kernel(); typename ei_lu_kernel_impl<MatrixType>::ReturnMatrixType m1kernel = lu.kernel();
typename ei_lu_image_impl <MatrixType>::ReturnMatrixType m1image = lu.image(); typename ei_lu_image_impl <MatrixType>::ReturnMatrixType m1image = lu.image(m1);
// std::cerr << rank << " " << lu.rank() << std::endl; // std::cerr << rank << " " << lu.rank() << std::endl;
VERIFY(rank == lu.rank()); VERIFY(rank == lu.rank());
@ -88,7 +88,7 @@ template<typename MatrixType> void lu_invertible()
VERIFY(lu.isInjective()); VERIFY(lu.isInjective());
VERIFY(lu.isSurjective()); VERIFY(lu.isSurjective());
VERIFY(lu.isInvertible()); VERIFY(lu.isInvertible());
VERIFY(lu.image().lu().isInvertible()); VERIFY(lu.image(m1).lu().isInvertible());
m3 = MatrixType::Random(size,size); m3 = MatrixType::Random(size,size);
m2 = lu.solve(m3); m2 = lu.solve(m3);
VERIFY_IS_APPROX(m3, m1*m2); VERIFY_IS_APPROX(m3, m1*m2);
@ -104,7 +104,7 @@ template<typename MatrixType> void lu_verify_assert()
VERIFY_RAISES_ASSERT(lu.permutationP()) VERIFY_RAISES_ASSERT(lu.permutationP())
VERIFY_RAISES_ASSERT(lu.permutationQ()) VERIFY_RAISES_ASSERT(lu.permutationQ())
VERIFY_RAISES_ASSERT(lu.kernel()) VERIFY_RAISES_ASSERT(lu.kernel())
VERIFY_RAISES_ASSERT(lu.image()) VERIFY_RAISES_ASSERT(lu.image(tmp))
VERIFY_RAISES_ASSERT(lu.solve(tmp)) VERIFY_RAISES_ASSERT(lu.solve(tmp))
VERIFY_RAISES_ASSERT(lu.determinant()) VERIFY_RAISES_ASSERT(lu.determinant())
VERIFY_RAISES_ASSERT(lu.rank()) VERIFY_RAISES_ASSERT(lu.rank())