GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-09-17 20:03:17 +08:00
Code Issues Packages Projects Releases Wiki Activity

add reconstructedMatrix() to LLT, and LUs

Browse Source 
=> they show that some improvements have still to be done
   for permutations, tr*tr, trapezoidal matrices
This commit is contained in:
Gael Guennebaud 2010-02-24 19:16:10 +01:00
parent a7e4c0f825
commit 7c98c04412
6 changed files with 87 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
Eigen/src/Cholesky/LDLT.h
View File

@ -155,7 +155,7 @@ template<typename _MatrixType> class LDLT
return m_matrix;
}
const MatrixType reconstructedMatrix() const;
MatrixType reconstructedMatrix() const;
inline int rows() const { return m_matrix.rows(); }
inline int cols() const { return m_matrix.cols(); }
@ -324,7 +324,7 @@ bool LDLT<MatrixType>::solveInPlace(MatrixBase<Derived> &bAndX) const
* i.e., it returns the product: P^T L D L^* P.
* This function is provided for debug purpose. */
template<typename MatrixType>
const MatrixType LDLT<MatrixType>::reconstructedMatrix() const
MatrixType LDLT<MatrixType>::reconstructedMatrix() const
{
ei_assert(m_isInitialized && "LDLT is not initialized.");
const int size = m_matrix.rows();

12
Eigen/src/Cholesky/LLT.h
View File

@ -133,6 +133,8 @@ template<typename _MatrixType, int _UpLo> class LLT
return m_matrix;
}
MatrixType reconstructedMatrix() const;
inline int rows() const { return m_matrix.rows(); }
inline int cols() const { return m_matrix.cols(); }
@ -295,6 +297,16 @@ bool LLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const
return true;
}
/** \returns the matrix represented by the decomposition,
* i.e., it returns the product: L L^*.
* This function is provided for debug purpose. */
template<typename MatrixType, int _UpLo>
MatrixType LLT<MatrixType,_UpLo>::reconstructedMatrix() const
{
ei_assert(m_isInitialized && "LLT is not initialized.");
return matrixL() * matrixL().adjoint().toDenseMatrix();
}
/** \cholesky_module
* \returns the LLT decomposition of \c *this
*/

29
Eigen/src/LU/FullPivLU.h
View File

@ -361,6 +361,8 @@ template<typename _MatrixType> class FullPivLU
(*this, MatrixType::Identity(m_lu.rows(), m_lu.cols()));
}
MatrixType reconstructedMatrix() const;
inline int rows() const { return m_lu.rows(); }
inline int cols() const { return m_lu.cols(); }
@ -487,6 +489,33 @@ typename ei_traits<MatrixType>::Scalar FullPivLU<MatrixType>::determinant() cons
return Scalar(m_det_pq) * Scalar(m_lu.diagonal().prod());
}
/** \returns the matrix represented by the decomposition,
* i.e., it returns the product: P^{-1} L U Q^{-1}.
* This function is provided for debug purpose. */
template<typename MatrixType>
MatrixType FullPivLU<MatrixType>::reconstructedMatrix() const
{
ei_assert(m_isInitialized && "LU is not initialized.");
const int smalldim = std::min(m_lu.rows(), m_lu.cols());
// LU
MatrixType res(m_lu.rows(),m_lu.cols());
// FIXME the .toDenseMatrix() should not be needed...
res = m_lu.corner(TopLeft,m_lu.rows(),smalldim)
.template triangularView<UnitLower>().toDenseMatrix()
* m_lu.corner(TopLeft,smalldim,m_lu.cols())
.template triangularView<Upper>().toDenseMatrix();
// P^{-1}(LU)
// FIXME implement inplace permutation
res = (m_p.inverse() * res).eval();
// (P^{-1}LU)Q^{-1}
// FIXME implement inplace permutation
res = (res * m_q.inverse()).eval();
return res;
}
/********* Implementation of kernel() **************************************************/
template<typename _MatrixType>

20
Eigen/src/LU/PartialPivLU.h
View File

@ -165,6 +165,8 @@ template<typename _MatrixType> class PartialPivLU
*/
typename ei_traits<MatrixType>::Scalar determinant() const;
MatrixType reconstructedMatrix() const;
inline int rows() const { return m_lu.rows(); }
inline int cols() const { return m_lu.cols(); }
@ -400,6 +402,24 @@ typename ei_traits<MatrixType>::Scalar PartialPivLU<MatrixType>::determinant() c
return Scalar(m_det_p) * m_lu.diagonal().prod();
}
/** \returns the matrix represented by the decomposition,
* i.e., it returns the product: P^{-1} L U.
* This function is provided for debug purpose. */
template<typename MatrixType>
MatrixType PartialPivLU<MatrixType>::reconstructedMatrix() const
{
ei_assert(m_isInitialized && "LU is not initialized.");
// LU
MatrixType res = m_lu.template triangularView<UnitLower>().toDenseMatrix()
* m_lu.template triangularView<Upper>();
// P^{-1}(LU)
// FIXME implement inplace permutation
res = (m_p.inverse() * res).eval();
return res;
}
/***** Implementation of solve() *****************************************************/
template<typename _MatrixType, typename Rhs>

7
test/cholesky.cpp
View File

@ -95,7 +95,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
{
LLT<SquareMatrixType,Lower> chollo(symmLo);
VERIFY_IS_APPROX(symm, chollo.matrixL().toDenseMatrix() * chollo.matrixL().adjoint().toDenseMatrix());
VERIFY_IS_APPROX(symm, chollo.reconstructedMatrix());
vecX = chollo.solve(vecB);
VERIFY_IS_APPROX(symm * vecX, vecB);
matX = chollo.solve(matB);
@ -103,7 +103,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
// test the upper mode
LLT<SquareMatrixType,Upper> cholup(symmUp);
VERIFY_IS_APPROX(symm, cholup.matrixL().toDenseMatrix() * cholup.matrixL().adjoint().toDenseMatrix());
VERIFY_IS_APPROX(symm, cholup.reconstructedMatrix());
vecX = cholup.solve(vecB);
VERIFY_IS_APPROX(symm * vecX, vecB);
matX = cholup.solve(matB);
@ -119,8 +119,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
{
LDLT<SquareMatrixType> ldlt(symm);
// TODO(keir): This doesn't make sense now that LDLT pivots.
//VERIFY_IS_APPROX(symm, ldlt.matrixL() * ldlt.vectorD().asDiagonal() * ldlt.matrixL().adjoint());
VERIFY_IS_APPROX(symm, ldlt.reconstructedMatrix());
vecX = ldlt.solve(vecB);
VERIFY_IS_APPROX(symm * vecX, vecB);
matX = ldlt.solve(matB);

21
test/lu.cpp
View File

@ -91,6 +91,7 @@ template<typename MatrixType> void lu_non_invertible()
KernelMatrixType m1kernel = lu.kernel();
ImageMatrixType m1image = lu.image(m1);
VERIFY_IS_APPROX(m1, lu.reconstructedMatrix());
VERIFY(rank == lu.rank());
VERIFY(cols - lu.rank() == lu.dimensionOfKernel());
VERIFY(!lu.isInjective());
@ -125,6 +126,7 @@ template<typename MatrixType> void lu_invertible()
lu.compute(m1);
} while(!lu.isInvertible());
VERIFY_IS_APPROX(m1, lu.reconstructedMatrix());
VERIFY(0 == lu.dimensionOfKernel());
VERIFY(lu.kernel().cols() == 1); // the kernel() should consist of a single (zero) column vector
VERIFY(size == lu.rank());
@ -138,6 +140,23 @@ template<typename MatrixType> void lu_invertible()
VERIFY_IS_APPROX(m2, lu.inverse()*m3);
}
template<typename MatrixType> void lu_partial_piv()
{
/* this test covers the following files:
PartialPivLU.h
*/
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
int rows = ei_random<int>(1,4);
int cols = rows;
MatrixType m1(cols, rows);
m1.setRandom();
PartialPivLU<MatrixType> plu(m1);
VERIFY_IS_APPROX(m1, plu.reconstructedMatrix());
}
template<typename MatrixType> void lu_verify_assert()
{
MatrixType tmp;
@ -180,6 +199,7 @@ void test_lu()
CALL_SUBTEST_4( lu_non_invertible<MatrixXd>() );
CALL_SUBTEST_4( lu_invertible<MatrixXd>() );
CALL_SUBTEST_4( lu_partial_piv<MatrixXd>() );
CALL_SUBTEST_4( lu_verify_assert<MatrixXd>() );
CALL_SUBTEST_5( lu_non_invertible<MatrixXcf>() );
@ -188,6 +208,7 @@ void test_lu()
CALL_SUBTEST_6( lu_non_invertible<MatrixXcd>() );
CALL_SUBTEST_6( lu_invertible<MatrixXcd>() );
CALL_SUBTEST_6( lu_partial_piv<MatrixXcd>() );
CALL_SUBTEST_6( lu_verify_assert<MatrixXcd>() );
CALL_SUBTEST_7(( lu_non_invertible<Matrix<float,Dynamic,16> >() ));