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

Migrate JacobiSVD to Solver

Browse Source
This commit is contained in:
Gael Guennebaud 2014-03-11 13:43:46 +01:00
parent 082f7ddc37
commit 7eefdb948c
2 changed files with 41 additions and 13 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -702,8 +702,8 @@ struct image_retval<FullPivLU<_MatrixType> >
#ifndef EIGEN_PARSED_BY_DOXYGEN #ifndef EIGEN_PARSED_BY_DOXYGEN
template<typename _MatrixType> template<typename _MatrixType>
template<typename RhsType, typename DstType> template<typename RhsType, typename DstType>
void FullPivLU<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const { void FullPivLU<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const
{
/* The decomposition PAQ = LU can be rewritten as A = P^{-1} L U Q^{-1}. /* The decomposition PAQ = LU can be rewritten as A = P^{-1} L U Q^{-1}.
* So we proceed as follows: * So we proceed as follows:
* Step 1: compute c = P * rhs. * Step 1: compute c = P * rhs.

50
Eigen/src/SVD/JacobiSVD.h
View File

@ -653,6 +653,16 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
* \note SVD solving is implicitly least-squares. Thus, this method serves both purposes of exact solving and least-squares solving. * \note SVD solving is implicitly least-squares. Thus, this method serves both purposes of exact solving and least-squares solving.
* In other words, the returned solution is guaranteed to minimize the Euclidean norm \f$ \Vert A x - b \Vert \f$. * In other words, the returned solution is guaranteed to minimize the Euclidean norm \f$ \Vert A x - b \Vert \f$.
*/ */
#ifdef EIGEN_TEST_EVALUATORS
template<typename Rhs>
inline const Solve<JacobiSVD, Rhs>
solve(const MatrixBase<Rhs>& b) const
{
eigen_assert(m_isInitialized && "JacobiSVD is not initialized.");
eigen_assert(computeU() && computeV() && "JacobiSVD::solve() requires both unitaries U and V to be computed (thin unitaries suffice).");
return Solve<JacobiSVD, Rhs>(*this, b.derived());
}
#else
template<typename Rhs> template<typename Rhs>
inline const internal::solve_retval<JacobiSVD, Rhs> inline const internal::solve_retval<JacobiSVD, Rhs>
solve(const MatrixBase<Rhs>& b) const solve(const MatrixBase<Rhs>& b) const
@ -661,6 +671,7 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
eigen_assert(computeU() && computeV() && "JacobiSVD::solve() requires both unitaries U and V to be computed (thin unitaries suffice)."); eigen_assert(computeU() && computeV() && "JacobiSVD::solve() requires both unitaries U and V to be computed (thin unitaries suffice).");
return internal::solve_retval<JacobiSVD, Rhs>(*this, b.derived()); return internal::solve_retval<JacobiSVD, Rhs>(*this, b.derived());
} }
#endif
/** \returns the number of singular values that are not exactly 0 */ /** \returns the number of singular values that are not exactly 0 */
Index nonzeroSingularValues() const Index nonzeroSingularValues() const
@ -734,6 +745,12 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
inline Index rows() const { return m_rows; } inline Index rows() const { return m_rows; }
inline Index cols() const { return m_cols; } inline Index cols() const { return m_cols; }
#ifndef EIGEN_PARSED_BY_DOXYGEN
template<typename RhsType, typename DstType>
EIGEN_DEVICE_FUNC
void _solve_impl(const RhsType &rhs, DstType &dst) const;
#endif
private: private:
void allocate(Index rows, Index cols, unsigned int computationOptions); void allocate(Index rows, Index cols, unsigned int computationOptions);
@ -912,7 +929,27 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig
return *this; return *this;
} }
#ifndef EIGEN_PARSED_BY_DOXYGEN
template<typename _MatrixType, int QRPreconditioner>
template<typename RhsType, typename DstType>
void JacobiSVD<_MatrixType,QRPreconditioner>::_solve_impl(const RhsType &rhs, DstType &dst) const
{
eigen_assert(rhs.rows() == rows());
// A = U S V^*
// So A^{-1} = V S^{-1} U^*
Matrix<Scalar, Dynamic, RhsType::ColsAtCompileTime, 0, _MatrixType::MaxRowsAtCompileTime, RhsType::MaxColsAtCompileTime> tmp;
Index l_rank = rank();
tmp.noalias() = m_matrixU.leftCols(l_rank).adjoint() * rhs;
tmp = m_singularValues.head(l_rank).asDiagonal().inverse() * tmp;
dst = m_matrixV.leftCols(l_rank) * tmp;
}
#endif
namespace internal { namespace internal {
#ifndef EIGEN_TEST_EVALUATORS
template<typename _MatrixType, int QRPreconditioner, typename Rhs> template<typename _MatrixType, int QRPreconditioner, typename Rhs>
struct solve_retval<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs> struct solve_retval<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs>
: solve_retval_base<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs> : solve_retval_base<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs>
@ -922,19 +959,10 @@ struct solve_retval<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs>
template<typename Dest> void evalTo(Dest& dst) const template<typename Dest> void evalTo(Dest& dst) const
{ {
eigen_assert(rhs().rows() == dec().rows()); dec()._solve_impl(rhs(), dst);
// A = U S V^*
// So A^{-1} = V S^{-1} U^*
Matrix<Scalar, Dynamic, Rhs::ColsAtCompileTime, 0, _MatrixType::MaxRowsAtCompileTime, Rhs::MaxColsAtCompileTime> tmp;
Index rank = dec().rank();
tmp.noalias() = dec().matrixU().leftCols(rank).adjoint() * rhs();
tmp = dec().singularValues().head(rank).asDiagonal().inverse() * tmp;
dst = dec().matrixV().leftCols(rank) * tmp;
} }
}; };
#endif
} // end namespace internal } // end namespace internal
#ifndef __CUDACC__ #ifndef __CUDACC__