From 50625834e64616d798e8a28a66f0eef638dc2801 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Sat, 12 Jan 2013 09:40:31 +0100 Subject: SparseQR: clean a bit the documentation, fix rows/cols methods, remove rowsQ methods and rename matrixQR to matrixR. --- Eigen/src/SparseQR/SparseQR.h | 187 ++++++++++++++++++++++-------------------- 1 file changed, 96 insertions(+), 91 deletions(-) (limited to 'Eigen/src/SparseQR') diff --git a/Eigen/src/SparseQR/SparseQR.h b/Eigen/src/SparseQR/SparseQR.h index e1016ac90..d78d5436d 100644 --- a/Eigen/src/SparseQR/SparseQR.h +++ b/Eigen/src/SparseQR/SparseQR.h @@ -34,30 +34,30 @@ namespace internal { } // End namespace internal /** - * \ingroup SparseQRSupport_Module - * \class SparseQR - * \brief Sparse left-looking QR factorization - * - * This class is used to perform a left-looking QR decomposition - * of sparse matrices. The result is then used to solve linear leasts_square systems. - * Clearly, a QR factorization is returned such that A*P = Q*R where : - * - * P is the column permutation. Use colsPermutation() to get it. - * - * Q is the orthogonal matrix represented as Householder reflectors. - * Use matrixQ() to get an expression and matrixQ().transpose() to get the transpose. - * You can then apply it to a vector. - * - * R is the sparse triangular factor. Use matrixQR() to get it as SparseMatrix. - * - * NOTE This is not a rank-revealing QR decomposition. - * - * \tparam _MatrixType The type of the sparse matrix A, must be a column-major SparseMatrix<> - * \tparam _OrderingType The fill-reducing ordering method. See \link OrderingMethods_Module - * OrderingMethods_Module \endlink for the list of built-in and external ordering methods. - * - * - */ + * \ingroup SparseQR_Module + * \class SparseQR + * \brief Sparse left-looking QR factorization + * + * This class is used to perform a left-looking QR decomposition + * of sparse matrices. The result is then used to solve linear leasts_square systems. + * Clearly, a QR factorization is returned such that A*P = Q*R where : + * + * P is the column permutation. Use colsPermutation() to get it. + * + * Q is the orthogonal matrix represented as Householder reflectors. + * Use matrixQ() to get an expression and matrixQ().transpose() to get the transpose. + * You can then apply it to a vector. + * + * R is the sparse triangular factor. Use matrixR() to get it as SparseMatrix. + * + * \note This is not a rank-revealing QR decomposition. + * + * \tparam _MatrixType The type of the sparse matrix A, must be a column-major SparseMatrix<> + * \tparam _OrderingType The fill-reducing ordering method. See the \link OrderingMethods_Module + * OrderingMethods \endlink module for the list of built-in and external ordering methods. + * + * + */ template class SparseQR { @@ -87,57 +87,56 @@ class SparseQR void analyzePattern(const MatrixType& mat); void factorize(const MatrixType& mat); - /** - * Get the number of rows of the triangular matrix. - */ - inline Index rows() const { return m_R.rows(); } + /** \returns the number of rows of the represented matrix. + */ + inline Index rows() const { return m_pmat.rows(); } - /** - * Get the number of columns of the triangular matrix. - */ - inline Index cols() const { return m_R.cols();} + /** \returns the number of columns of the represented matrix. + */ + inline Index cols() const { return m_pmat.cols();} - /** - * This is the number of rows in the input matrix and the Q matrix - */ - inline Index rowsQ() const {return m_pmat.rows(); } + /** \returns a const reference to the \b sparse upper triangular matrix R of the QR factorization. + */ + const MatrixType& matrixR() const { return m_R; } - /// Get the sparse triangular matrix R. It is a sparse matrix - MatrixType matrixQR() const - { - return m_R; - } - /// Get an expression of the matrix Q + /** \returns an expression of the matrix Q as products of sparse Householder reflectors. + * You can do the following to get an actual SparseMatrix representation of Q: + * \code + * SparseMatrix Q = SparseQR >(A).matrixQ(); + * \endcode + */ SparseQRMatrixQReturnType matrixQ() const - { - return SparseQRMatrixQReturnType(*this); - } + { return SparseQRMatrixQReturnType(*this); } - /// Get the permutation that was applied to columns of A - PermutationType colsPermutation() const + /** \returns a const reference to the fill-in reducing permutation that was applied to the columns of A + */ + const PermutationType& colsPermutation() const { eigen_assert(m_isInitialized && "Decomposition is not initialized."); return m_perm_c; } + + /** \internal */ template bool _solve(const MatrixBase &B, MatrixBase &dest) const { eigen_assert(m_isInitialized && "The factorization should be called first, use compute()"); - eigen_assert(this->rowsQ() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix"); - Index rank = this->matrixQR().cols(); + eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix"); + Index rank = this->matrixR().cols(); // Compute Q^T * b; dest = this->matrixQ().transpose() * B; // Solve with the triangular matrix R Dest y; - y = this->matrixQR().template triangularView().solve(dest.derived().topRows(rank)); + y = this->matrixR().template triangularView().solve(dest.derived().topRows(rank)); + // Apply the column permutation - if (m_perm_c.size()) - dest.topRows(rank) = colsPermutation().inverse() * y; - else - dest = y; + if (m_perm_c.size()) dest.topRows(rank) = colsPermutation().inverse() * y; + else dest = y; + m_info = Success; return true; } + /** \returns the solution X of \f$ A X = B \f$ using the current decomposition of A. * * \sa compute() @@ -146,13 +145,14 @@ class SparseQR inline const internal::solve_retval solve(const MatrixBase& B) const { eigen_assert(m_isInitialized && "The factorization should be called first, use compute()"); - eigen_assert(this->rowsQ() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix"); + eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix"); return internal::solve_retval(*this, B.derived()); } + /** \brief Reports whether previous computation was successful. * * \returns \c Success if computation was succesful, - * \c NumericalIssue if the QR factorization reports a problem + * \c NumericalIssue if the QR factorization reports a numerical problem * \c InvalidInput if the input matrix is invalid * * \sa iparm() @@ -162,30 +162,31 @@ class SparseQR eigen_assert(m_isInitialized && "Decomposition is not initialized."); return m_info; } + protected: bool m_isInitialized; bool m_analysisIsok; bool m_factorizationIsok; mutable ComputationInfo m_info; - QRMatrixType m_pmat; // Temporary matrix - QRMatrixType m_R; // The triangular factor matrix - QRMatrixType m_Q; // The orthogonal reflectors - ScalarVector m_hcoeffs; // The Householder coefficients - PermutationType m_perm_c; // Column permutation - PermutationType m_perm_r; // Column permutation - IndexVector m_etree; // Column elimination tree - IndexVector m_firstRowElt; // First element in each row - IndexVector m_found_diag_elem; // Existence of diagonal elements + QRMatrixType m_pmat; // Temporary matrix + QRMatrixType m_R; // The triangular factor matrix + QRMatrixType m_Q; // The orthogonal reflectors + ScalarVector m_hcoeffs; // The Householder coefficients + PermutationType m_perm_c; // Column permutation + PermutationType m_perm_r; // Column permutation + IndexVector m_etree; // Column elimination tree + IndexVector m_firstRowElt; // First element in each row + IndexVector m_found_diag_elem; // Existence of diagonal elements template friend struct SparseQR_QProduct; }; -/** - * \brief Preprocessing step of a QR factorization - * - * In this step, the fill-reducing permutation is computed and applied to the columns of A - * and the column elimination tree is computed as well. - * \note In this step it is assumed that there is no empty row in the matrix \p mat - */ + +/** \brief Preprocessing step of a QR factorization + * + * In this step, the fill-reducing permutation is computed and applied to the columns of A + * and the column elimination tree is computed as well. Only the sparcity pattern of \a mat is exploited. + * \note In this step it is assumed that there is no empty row in the matrix \a mat + */ template void SparseQR::analyzePattern(const MatrixType& mat) { @@ -195,18 +196,17 @@ void SparseQR::analyzePattern(const MatrixType& mat) Index n = mat.cols(); Index m = mat.rows(); // Permute the input matrix... only the column pointers are permuted + // FIXME: directly send "m_perm.inverse() * mat" to coletree -> need an InnerIterator to the sparse-permutation-product expression. m_pmat = mat; m_pmat.uncompress(); for (int i = 0; i < n; i++) { Index p = m_perm_c.size() ? m_perm_c.indices()(i) : i; -// m_pmat.outerIndexPtr()[i] = mat.outerIndexPtr()[p]; -// m_pmat.innerNonZeroPtr()[i] = mat.outerIndexPtr()[p+1] - mat.outerIndexPtr()[p]; m_pmat.outerIndexPtr()[p] = mat.outerIndexPtr()[i]; m_pmat.innerNonZeroPtr()[p] = mat.outerIndexPtr()[i+1] - mat.outerIndexPtr()[i]; } // Compute the column elimination tree of the permuted matrix - internal::coletree(m_pmat, m_etree, m_firstRowElt/*, m_found_diag_elem*/); + internal::coletree(m_pmat, m_etree, m_firstRowElt); m_R.resize(n, n); m_Q.resize(m, m); @@ -217,22 +217,24 @@ void SparseQR::analyzePattern(const MatrixType& mat) m_analysisIsok = true; } -/** - * \brief Perform the numerical QR factorization of the input matrix - * - * \param mat The sparse column-major matrix - */ +/** \brief Perform the numerical QR factorization of the input matrix + * + * The function SparseQR::analyzePattern(const MatrixType&) must have been called beforehand with + * a matrix having the same sparcity pattern than \a mat. + * + * \param mat The sparse column-major matrix + */ template void SparseQR::factorize(const MatrixType& mat) { eigen_assert(m_analysisIsok && "analyzePattern() should be called before this step"); Index m = mat.rows(); Index n = mat.cols(); - IndexVector mark(m); mark.setConstant(-1);// Record the visited nodes - IndexVector Ridx(n),Qidx(m); // Store temporarily the row indexes for the current column of R and Q - Index nzcolR, nzcolQ; // Number of nonzero for the current column of R and Q + IndexVector mark(m); mark.setConstant(-1); // Record the visited nodes + IndexVector Ridx(n), Qidx(m); // Store temporarily the row indexes for the current column of R and Q + Index nzcolR, nzcolQ; // Number of nonzero for the current column of R and Q Index pcol; - ScalarVector tval(m); tval.setZero(); // Temporary vector + ScalarVector tval(m); tval.setZero(); // Temporary vector IndexVector iperm(m); bool found_diag; if (m_perm_c.size()) @@ -290,7 +292,7 @@ void SparseQR::factorize(const MatrixType& mat) } } - //Browse all the indexes of R(:,col) in reverse order + // Browse all the indexes of R(:,col) in reverse order for (Index i = nzcolR-1; i >= 0; i--) { Index curIdx = Ridx(i); @@ -311,7 +313,7 @@ void SparseQR::factorize(const MatrixType& mat) m_R.insertBackByOuterInnerUnordered(col, curIdx) = tval(curIdx); tval(curIdx) = Scalar(0.); - //Detect fill-in for the current column of Q + // Detect fill-in for the current column of Q if(m_etree(curIdx) == col) { for (typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq) @@ -410,7 +412,7 @@ struct SparseQR_QProduct : ReturnByValue void evalTo(DesType& res) const { - Index m = m_qr.rowsQ(); + Index m = m_qr.rows(); Index n = m_qr.cols(); if (m_transpose) { @@ -447,8 +449,8 @@ struct SparseQR_QProduct : ReturnByValue -struct SparseQRMatrixQReturnType{ - +struct SparseQRMatrixQReturnType +{ SparseQRMatrixQReturnType(const SparseQRType& qr) : m_qr(qr) {} template SparseQR_QProduct operator*(const MatrixBase& other) @@ -468,7 +470,8 @@ struct SparseQRMatrixQReturnType{ }; template -struct SparseQRMatrixQTransposeReturnType{ +struct SparseQRMatrixQTransposeReturnType +{ SparseQRMatrixQTransposeReturnType(const SparseQRType& qr) : m_qr(qr) {} template SparseQR_QProduct operator*(const MatrixBase& other) @@ -477,5 +480,7 @@ struct SparseQRMatrixQTransposeReturnType{ } const SparseQRType& m_qr; }; + } // end namespace Eigen -#endif \ No newline at end of file + +#endif -- cgit v1.2.3