diff --git a/Eigen/src/QR/ColPivHouseholderQR.h b/Eigen/src/QR/ColPivHouseholderQR.h index ea63e6f8e..01f6ac055 100644 --- a/Eigen/src/QR/ColPivHouseholderQR.h +++ b/Eigen/src/QR/ColPivHouseholderQR.h @@ -127,6 +127,10 @@ template class ColPivHouseholderQR } HouseholderSequenceType householderQ(void) const; + HouseholderSequenceType matrixQ(void) const + { + return householderQ(); + } /** \returns a reference to the matrix where the Householder QR decomposition is stored */ diff --git a/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h b/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h index 2647d22f0..a3880c9f8 100644 --- a/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h +++ b/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h @@ -71,8 +71,12 @@ class SPQR cholmod_l_start(&m_cc); } - SPQR(const _MatrixType& matrix) : SPQR() + SPQR(const _MatrixType& matrix) + : m_ordering(SPQR_ORDERING_DEFAULT), + m_allow_tol(SPQR_DEFAULT_TOL), + m_tolerance (NumTraits::epsilon()) { + cholmod_l_start(&m_cc); compute(matrix); } @@ -102,6 +106,32 @@ class SPQR m_info = Success; m_isInitialized = true; } + /** + * Get the number of rows of the triangular matrix. + */ + inline Index rows() const { return m_cR->nrow; } + + /** + * Get the number of columns of the triangular matrix. + */ + inline Index cols() const { return m_cR->ncol; } + /** + * This is the number of rows in the input matrix and the Q matrix + */ + inline Index rowsQ() const {return m_HTau->nrow; } + /** \returns the solution X of \f$ A X = B \f$ using the current decomposition of A. + * + * \sa compute() + */ + template + inline const internal::solve_retval solve(const MatrixBase& B) const + { + eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()"); + eigen_assert(rows()==B.rows() + && "SPQR::solve(): invalid number of rows of the right hand side matrix B"); + return internal::solve_retval(*this, B.derived()); + } + template void _solve(const MatrixBase &b, MatrixBase &dest) const { @@ -109,8 +139,6 @@ class SPQR eigen_assert(b.cols()==1 && "This method is for vectors only"); //Compute Q^T * b - // NOTE : We may have called directly the corresponding routines in SPQR codes. - // This version is used to test directly the corresponding part of the code dest = matrixQ().transpose() * b; // Solves with the triangular matrix R @@ -195,9 +223,12 @@ template struct SPQR_QProduct : ReturnByValue > { typedef typename SPQRType::Scalar Scalar; + typedef typename SPQRType::Index Index; //Define the constructor to get reference to argument types SPQR_QProduct(const SPQRType& spqr, const Derived& other, bool transpose) : m_spqr(spqr),m_other(other),m_transpose(transpose) {} + inline Index rows() const { return m_transpose ? m_spqr.rowsQ() : m_spqr.cols(); } + inline Index cols() const { return m_other.cols(); } // Assign to a vector template void evalTo(ResType& res) const @@ -225,6 +256,10 @@ struct SPQRMatrixQReturnType{ { return SPQR_QProduct(m_spqr,other.derived(),false); } + SPQRMatrixQTransposeReturnType adjoint() const + { + return SPQRMatrixQTransposeReturnType(m_spqr); + } // To use for operations with the transpose of Q SPQRMatrixQTransposeReturnType transpose() const { @@ -243,5 +278,23 @@ struct SPQRMatrixQTransposeReturnType{ } const SPQRType& m_spqr; }; + +namespace internal { + +template +struct solve_retval, Rhs> + : solve_retval_base, Rhs> +{ + typedef SPQR<_MatrixType> Dec; + EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) + + template void evalTo(Dest& dst) const + { + dec()._solve(rhs(),dst); + } +}; + +} // end namespace internal + }// End namespace Eigen #endif \ No newline at end of file diff --git a/test/spqr_support.cpp b/test/spqr_support.cpp index 9ff4f9cba..fbfd4c58b 100644 --- a/test/spqr_support.cpp +++ b/test/spqr_support.cpp @@ -44,7 +44,7 @@ template void test_spqr_scalar() exit(0); return; } - solver._solve(b,x); + x = solver.solve(b); if (solver.info() != Success) { std::cerr << "sparse QR factorization failed\n";