From 1f398dfc826c2427375e6aa6a15bfe23383d76f7 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 1 Sep 2014 18:31:54 +0200 Subject: [PATCH] Factorize *SVD::solve to SVDBase --- Eigen/src/SVD/JacobiSVD.h | 71 --------------------------------- Eigen/src/SVD/SVDBase.h | 70 ++++++++++++++++++++++++++++++++ unsupported/test/CMakeLists.txt | 2 +- 3 files changed, 71 insertions(+), 72 deletions(-) diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h index 3d778516a..2b9d03c2b 100644 --- a/Eigen/src/SVD/JacobiSVD.h +++ b/Eigen/src/SVD/JacobiSVD.h @@ -592,47 +592,12 @@ template class JacobiSVD return compute(matrix, m_computationOptions); } - /** \returns a (least squares) solution of \f$ A x = b \f$ using the current SVD decomposition of A. - * - * \param b the right-hand-side of the equation to solve. - * - * \note Solving requires both U and V to be computed. Thin U and V are enough, there is no need for full U or V. - * - * \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$. - */ -#ifdef EIGEN_TEST_EVALUATORS - template - inline const Solve - solve(const MatrixBase& 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(*this, b.derived()); - } -#else - template - inline const internal::solve_retval - solve(const MatrixBase& 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 internal::solve_retval(*this, b.derived()); - } -#endif - using Base::computeU; using Base::computeV; using Base::rows; using Base::cols; using Base::rank; - #ifndef EIGEN_PARSED_BY_DOXYGEN - template - EIGEN_DEVICE_FUNC - void _solve_impl(const RhsType &rhs, DstType &dst) const; - #endif - private: void allocate(Index rows, Index cols, unsigned int computationOptions); @@ -817,42 +782,6 @@ JacobiSVD::compute(const MatrixType& matrix, unsig return *this; } -#ifndef EIGEN_PARSED_BY_DOXYGEN -template -template -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 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 { -#ifndef EIGEN_TEST_EVALUATORS -template -struct solve_retval, Rhs> - : solve_retval_base, Rhs> -{ - typedef JacobiSVD<_MatrixType, QRPreconditioner> JacobiSVDType; - EIGEN_MAKE_SOLVE_HELPERS(JacobiSVDType,Rhs) - - template void evalTo(Dest& dst) const - { - dec()._solve_impl(rhs(), dst); - } -}; -#endif -} // end namespace internal - #ifndef __CUDACC__ /** \svd_module * diff --git a/Eigen/src/SVD/SVDBase.h b/Eigen/src/SVD/SVDBase.h index 61b01fb8a..a4bb97f4b 100644 --- a/Eigen/src/SVD/SVDBase.h +++ b/Eigen/src/SVD/SVDBase.h @@ -190,6 +190,41 @@ public: inline Index rows() const { return m_rows; } inline Index cols() const { return m_cols; } + + /** \returns a (least squares) solution of \f$ A x = b \f$ using the current SVD decomposition of A. + * + * \param b the right-hand-side of the equation to solve. + * + * \note Solving requires both U and V to be computed. Thin U and V are enough, there is no need for full U or V. + * + * \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$. + */ +#ifdef EIGEN_TEST_EVALUATORS + template + inline const Solve + solve(const MatrixBase& b) const + { + eigen_assert(m_isInitialized && "SVD is not initialized."); + eigen_assert(computeU() && computeV() && "SVD::solve() requires both unitaries U and V to be computed (thin unitaries suffice)."); + return Solve(derived(), b.derived()); + } +#else + template + inline const internal::solve_retval + solve(const MatrixBase& b) const + { + eigen_assert(m_isInitialized && "SVD is not initialized."); + eigen_assert(computeU() && computeV() && "SVD::solve() requires both unitaries U and V to be computed (thin unitaries suffice)."); + return internal::solve_retval(*this, b.derived()); + } +#endif + + #ifndef EIGEN_PARSED_BY_DOXYGEN + template + EIGEN_DEVICE_FUNC + void _solve_impl(const RhsType &rhs, DstType &dst) const; + #endif protected: // return true if already allocated @@ -220,6 +255,41 @@ protected: }; +#ifndef EIGEN_PARSED_BY_DOXYGEN +template +template +void SVDBase::_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 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 { +#ifndef EIGEN_TEST_EVALUATORS +template +struct solve_retval, Rhs> + : solve_retval_base, Rhs> +{ + typedef SVDBase SVDType; + EIGEN_MAKE_SOLVE_HELPERS(SVDType,Rhs) + + template void evalTo(Dest& dst) const + { + dec().derived()._solve_impl(rhs(), dst); + } +}; +#endif +} // end namespace internal template bool SVDBase::allocate(Index rows, Index cols, unsigned int computationOptions) diff --git a/unsupported/test/CMakeLists.txt b/unsupported/test/CMakeLists.txt index 5e9e74050..970a05bbd 100644 --- a/unsupported/test/CMakeLists.txt +++ b/unsupported/test/CMakeLists.txt @@ -100,8 +100,8 @@ ei_add_test(splines) ei_add_test(gmres) ei_add_test(minres) ei_add_test(levenberg_marquardt) -if(EIGEN_TEST_NO_EVALUATORS) ei_add_test(bdcsvd) +if(EIGEN_TEST_NO_EVALUATORS) ei_add_test(kronecker_product) endif()