diff options
-rw-r--r-- | Eigen/src/QR/CompleteOrthogonalDecomposition.h | 41 | ||||
-rw-r--r-- | test/qr_colpivoting.cpp | 3 |
2 files changed, 26 insertions, 18 deletions
diff --git a/Eigen/src/QR/CompleteOrthogonalDecomposition.h b/Eigen/src/QR/CompleteOrthogonalDecomposition.h index bee5bf47e..9bc768b7c 100644 --- a/Eigen/src/QR/CompleteOrthogonalDecomposition.h +++ b/Eigen/src/QR/CompleteOrthogonalDecomposition.h @@ -70,6 +70,7 @@ class CompleteOrthogonalDecomposition { MatrixType, typename internal::remove_all< typename HCoeffsType::ConjugateReturnType>::type> HouseholderSequenceType; + typedef typename MatrixType::PlainObject PlainObject; private: typedef typename PermutationType::Index PermIndexType; @@ -246,26 +247,15 @@ class CompleteOrthogonalDecomposition { */ inline bool isInvertible() const { return m_cpqr.isInvertible(); } - /** \returns the inverse of the matrix of which *this is the complete + /** \returns the pseudo-inverse of the matrix of which *this is the complete * orthogonal decomposition. - * - * \note If this matrix is not invertible, the returned matrix has undefined - * coefficients. Use isInvertible() to first determine whether this matrix is - * invertible. + * \warning: Do not compute \c this->inverse()*rhs to solve a linear systems. + * It is more efficient and numerically stable to call \c this->solve(rhs). */ - - // TODO(rmlarsen): Add method for pseudo-inverse. - // inline const - // internal::solve_retval<CompleteOrthogonalDecomposition, typename - // MatrixType::IdentityReturnType> - // inverse() const - // { - // eigen_assert(m_isInitialized && "CompleteOrthogonalDecomposition is not - // initialized."); - // return internal::solve_retval<CompleteOrthogonalDecomposition,typename - // MatrixType::IdentityReturnType> - // (*this, MatrixType::Identity(m_cpqr.rows(), m_cpqr.cols())); - // } + inline const Inverse<CompleteOrthogonalDecomposition> inverse() const + { + return Inverse<CompleteOrthogonalDecomposition>(*this); + } inline Index rows() const { return m_cpqr.rows(); } inline Index cols() const { return m_cpqr.cols(); } @@ -513,6 +503,21 @@ void CompleteOrthogonalDecomposition<_MatrixType>::_solve_impl( } #endif +namespace internal { + +template<typename DstXprType, typename MatrixType, typename Scalar> +struct Assignment<DstXprType, Inverse<CompleteOrthogonalDecomposition<MatrixType> >, internal::assign_op<Scalar>, Dense2Dense, Scalar> +{ + typedef CompleteOrthogonalDecomposition<MatrixType> CodType; + typedef Inverse<CodType> SrcXprType; + static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar> &) + { + dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.rows())); + } +}; + +} // end namespace internal + /** \returns the matrix Q as a sequence of householder transformations */ template <typename MatrixType> typename CompleteOrthogonalDecomposition<MatrixType>::HouseholderSequenceType diff --git a/test/qr_colpivoting.cpp b/test/qr_colpivoting.cpp index 16f80d8b5..d8672ce33 100644 --- a/test/qr_colpivoting.cpp +++ b/test/qr_colpivoting.cpp @@ -57,6 +57,9 @@ void cod() { JacobiSVD<MatrixType> svd(matrix, ComputeThinU | ComputeThinV); MatrixType svd_solution = svd.solve(rhs); VERIFY_IS_APPROX(cod_solution, svd_solution); + + MatrixType pinv = cod.inverse(); + VERIFY_IS_APPROX(cod_solution, pinv * rhs); } template <typename MatrixType, int Cols2> |