GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-11 03:09:01 +08:00
Code Issues Packages Projects Releases Wiki Activity

SparseQR: clean a bit the documentation, fix rows/cols methods, remove rowsQ methods and rename matrixQR to matrixR.

Browse Source
This commit is contained in:
Gael Guennebaud 2013-01-12 09:40:31 +01:00
parent 581e389c54
commit 50625834e6
3 changed files with 100 additions and 93 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
Eigen/SparseQR
View File

@ -6,12 +6,12 @@
#include "src/Core/util/DisableStupidWarnings.h" #include "src/Core/util/DisableStupidWarnings.h"
/** \defgroup SparseQR_Module SparseQR module /** \defgroup SparseQR_Module SparseQR module
* * \brief Provides QR decomposition for sparse matrices
* *
* This module provides a simplicial version of the left-looking Sparse QR decomposition. * This module provides a simplicial version of the left-looking Sparse QR decomposition.
* The columns of the input matrix should be reordered to limit the fill-in during the * The columns of the input matrix should be reordered to limit the fill-in during the
* decomposition. Built-in methods (COLAMD, AMD) or external methods (METIS) can be used to this end. * decomposition. Built-in methods (COLAMD, AMD) or external methods (METIS) can be used to this end.
* See \link OrderingMethods_Module OrderingMethods_Module \endlink for the list * See the \link OrderingMethods_Module OrderingMethods\endlink module for the list
* of built-in and external ordering methods. * of built-in and external ordering methods.
* *
* \code * \code

187
Eigen/src/SparseQR/SparseQR.h
View File

@ -34,30 +34,30 @@ namespace internal {
} // End namespace internal } // End namespace internal
/** /**
* \ingroup SparseQRSupport_Module * \ingroup SparseQR_Module
* \class SparseQR * \class SparseQR
* \brief Sparse left-looking QR factorization * \brief Sparse left-looking QR factorization
* *
* This class is used to perform a left-looking QR decomposition * 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. * 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 : * Clearly, a QR factorization is returned such that A*P = Q*R where :
* *
* P is the column permutation. Use colsPermutation() to get it. * P is the column permutation. Use colsPermutation() to get it.
* *
* Q is the orthogonal matrix represented as Householder reflectors. * Q is the orthogonal matrix represented as Householder reflectors.
* Use matrixQ() to get an expression and matrixQ().transpose() to get the transpose. * Use matrixQ() to get an expression and matrixQ().transpose() to get the transpose.
* You can then apply it to a vector. * You can then apply it to a vector.
* *
* R is the sparse triangular factor. Use matrixQR() to get it as SparseMatrix. * R is the sparse triangular factor. Use matrixR() to get it as SparseMatrix.
* *
* NOTE This is not a rank-revealing QR decomposition. * \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 _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 * \tparam _OrderingType The fill-reducing ordering method. See the \link OrderingMethods_Module
* OrderingMethods_Module \endlink for the list of built-in and external ordering methods. * OrderingMethods \endlink module for the list of built-in and external ordering methods.
* *
* *
*/ */
template<typename _MatrixType, typename _OrderingType> template<typename _MatrixType, typename _OrderingType>
class SparseQR class SparseQR
{ {
@ -87,57 +87,56 @@ class SparseQR
void analyzePattern(const MatrixType& mat); void analyzePattern(const MatrixType& mat);
void factorize(const MatrixType& mat); void factorize(const MatrixType& mat);
/** /** \returns the number of rows of the represented matrix.
* Get the number of rows of the triangular matrix. */
*/ inline Index rows() const { return m_pmat.rows(); }
inline Index rows() const { return m_R.rows(); }
/** /** \returns the number of columns of the represented matrix.
* Get the number of columns of the triangular matrix. */
*/ inline Index cols() const { return m_pmat.cols();}
inline Index cols() const { return m_R.cols();}
/** /** \returns a const reference to the \b sparse upper triangular matrix R of the QR factorization.
* This is the number of rows in the input matrix and the Q matrix */
*/ const MatrixType& matrixR() const { return m_R; }
inline Index rowsQ() const {return m_pmat.rows(); }
/// Get the sparse triangular matrix R. It is a sparse matrix /** \returns an expression of the matrix Q as products of sparse Householder reflectors.
MatrixType matrixQR() const * You can do the following to get an actual SparseMatrix representation of Q:
{ * \code
return m_R; * SparseMatrix<double> Q = SparseQR<SparseMatrix<double> >(A).matrixQ();
} * \endcode
/// Get an expression of the matrix Q */
SparseQRMatrixQReturnType<SparseQR> matrixQ() const SparseQRMatrixQReturnType<SparseQR> matrixQ() const
{ { return SparseQRMatrixQReturnType<SparseQR>(*this); }
return SparseQRMatrixQReturnType<SparseQR>(*this);
}
/// Get the permutation that was applied to columns of A /** \returns a const reference to the fill-in reducing permutation that was applied to the columns of A
PermutationType colsPermutation() const */
const PermutationType& colsPermutation() const
{ {
eigen_assert(m_isInitialized && "Decomposition is not initialized."); eigen_assert(m_isInitialized && "Decomposition is not initialized.");
return m_perm_c; return m_perm_c;
} }
/** \internal */
template<typename Rhs, typename Dest> template<typename Rhs, typename Dest>
bool _solve(const MatrixBase<Rhs> &B, MatrixBase<Dest> &dest) const bool _solve(const MatrixBase<Rhs> &B, MatrixBase<Dest> &dest) const
{ {
eigen_assert(m_isInitialized && "The factorization should be called first, use compute()"); 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");
Index rank = this->matrixQR().cols(); Index rank = this->matrixR().cols();
// Compute Q^T * b; // Compute Q^T * b;
dest = this->matrixQ().transpose() * B; dest = this->matrixQ().transpose() * B;
// Solve with the triangular matrix R // Solve with the triangular matrix R
Dest y; Dest y;
y = this->matrixQR().template triangularView<Upper>().solve(dest.derived().topRows(rank)); y = this->matrixR().template triangularView<Upper>().solve(dest.derived().topRows(rank));
// Apply the column permutation // Apply the column permutation
if (m_perm_c.size()) if (m_perm_c.size()) dest.topRows(rank) = colsPermutation().inverse() * y;
dest.topRows(rank) = colsPermutation().inverse() * y; else dest = y;
else
dest = y;
m_info = Success; m_info = Success;
return true; return true;
} }
/** \returns the solution X of \f$ A X = B \f$ using the current decomposition of A. /** \returns the solution X of \f$ A X = B \f$ using the current decomposition of A.
* *
* \sa compute() * \sa compute()
@ -146,13 +145,14 @@ class SparseQR
inline const internal::solve_retval<SparseQR, Rhs> solve(const MatrixBase<Rhs>& B) const inline const internal::solve_retval<SparseQR, Rhs> solve(const MatrixBase<Rhs>& B) const
{ {
eigen_assert(m_isInitialized && "The factorization should be called first, use compute()"); 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<SparseQR, Rhs>(*this, B.derived()); return internal::solve_retval<SparseQR, Rhs>(*this, B.derived());
} }
/** \brief Reports whether previous computation was successful. /** \brief Reports whether previous computation was successful.
* *
* \returns \c Success if computation was succesful, * \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 * \c InvalidInput if the input matrix is invalid
* *
* \sa iparm() * \sa iparm()
@ -162,30 +162,31 @@ class SparseQR
eigen_assert(m_isInitialized && "Decomposition is not initialized."); eigen_assert(m_isInitialized && "Decomposition is not initialized.");
return m_info; return m_info;
} }
protected: protected:
bool m_isInitialized; bool m_isInitialized;
bool m_analysisIsok; bool m_analysisIsok;
bool m_factorizationIsok; bool m_factorizationIsok;
mutable ComputationInfo m_info; mutable ComputationInfo m_info;
QRMatrixType m_pmat; // Temporary matrix QRMatrixType m_pmat; // Temporary matrix
QRMatrixType m_R; // The triangular factor matrix QRMatrixType m_R; // The triangular factor matrix
QRMatrixType m_Q; // The orthogonal reflectors QRMatrixType m_Q; // The orthogonal reflectors
ScalarVector m_hcoeffs; // The Householder coefficients ScalarVector m_hcoeffs; // The Householder coefficients
PermutationType m_perm_c; // Column permutation PermutationType m_perm_c; // Column permutation
PermutationType m_perm_r; // Column permutation PermutationType m_perm_r; // Column permutation
IndexVector m_etree; // Column elimination tree IndexVector m_etree; // Column elimination tree
IndexVector m_firstRowElt; // First element in each row IndexVector m_firstRowElt; // First element in each row
IndexVector m_found_diag_elem; // Existence of diagonal elements IndexVector m_found_diag_elem; // Existence of diagonal elements
template <typename, typename > friend struct SparseQR_QProduct; template <typename, typename > friend struct SparseQR_QProduct;
}; };
/**
* \brief Preprocessing step of a QR factorization /** \brief Preprocessing step of a QR factorization
* *
* In this step, the fill-reducing permutation is computed and applied to the columns of A * 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. * 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 \p mat * \note In this step it is assumed that there is no empty row in the matrix \a mat
*/ */
template <typename MatrixType, typename OrderingType> template <typename MatrixType, typename OrderingType>
void SparseQR<MatrixType,OrderingType>::analyzePattern(const MatrixType& mat) void SparseQR<MatrixType,OrderingType>::analyzePattern(const MatrixType& mat)
{ {
@ -195,18 +196,17 @@ void SparseQR<MatrixType,OrderingType>::analyzePattern(const MatrixType& mat)
Index n = mat.cols(); Index n = mat.cols();
Index m = mat.rows(); Index m = mat.rows();
// Permute the input matrix... only the column pointers are permuted // 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 = mat;
m_pmat.uncompress(); m_pmat.uncompress();
for (int i = 0; i < n; i++) for (int i = 0; i < n; i++)
{ {
Index p = m_perm_c.size() ? m_perm_c.indices()(i) : 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.outerIndexPtr()[p] = mat.outerIndexPtr()[i];
m_pmat.innerNonZeroPtr()[p] = mat.outerIndexPtr()[i+1] - mat.outerIndexPtr()[i]; m_pmat.innerNonZeroPtr()[p] = mat.outerIndexPtr()[i+1] - mat.outerIndexPtr()[i];
} }
// Compute the column elimination tree of the permuted matrix // 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_R.resize(n, n);
m_Q.resize(m, m); m_Q.resize(m, m);
@ -217,22 +217,24 @@ void SparseQR<MatrixType,OrderingType>::analyzePattern(const MatrixType& mat)
m_analysisIsok = true; m_analysisIsok = true;
} }
/** /** \brief Perform the numerical QR factorization of the input matrix
* \brief Perform the numerical QR factorization of the input matrix *
* * The function SparseQR::analyzePattern(const MatrixType&) must have been called beforehand with
* \param mat The sparse column-major matrix * a matrix having the same sparcity pattern than \a mat.
*/ *
* \param mat The sparse column-major matrix
*/
template <typename MatrixType, typename OrderingType> template <typename MatrixType, typename OrderingType>
void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat) void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat)
{ {
eigen_assert(m_analysisIsok && "analyzePattern() should be called before this step"); eigen_assert(m_analysisIsok && "analyzePattern() should be called before this step");
Index m = mat.rows(); Index m = mat.rows();
Index n = mat.cols(); Index n = mat.cols();
IndexVector mark(m); mark.setConstant(-1);// Record the visited nodes 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 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 nzcolR, nzcolQ; // Number of nonzero for the current column of R and Q
Index pcol; Index pcol;
ScalarVector tval(m); tval.setZero(); // Temporary vector ScalarVector tval(m); tval.setZero(); // Temporary vector
IndexVector iperm(m); IndexVector iperm(m);
bool found_diag; bool found_diag;
if (m_perm_c.size()) if (m_perm_c.size())
@ -290,7 +292,7 @@ void SparseQR<MatrixType,OrderingType>::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--) for (Index i = nzcolR-1; i >= 0; i--)
{ {
Index curIdx = Ridx(i); Index curIdx = Ridx(i);
@ -311,7 +313,7 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat)
m_R.insertBackByOuterInnerUnordered(col, curIdx) = tval(curIdx); m_R.insertBackByOuterInnerUnordered(col, curIdx) = tval(curIdx);
tval(curIdx) = Scalar(0.); 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) if(m_etree(curIdx) == col)
{ {
for (typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq) for (typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq)
@ -410,7 +412,7 @@ struct SparseQR_QProduct : ReturnByValue<SparseQR_QProduct<SparseQRType, Derived
template<typename DesType> template<typename DesType>
void evalTo(DesType& res) const void evalTo(DesType& res) const
{ {
Index m = m_qr.rowsQ(); Index m = m_qr.rows();
Index n = m_qr.cols(); Index n = m_qr.cols();
if (m_transpose) if (m_transpose)
{ {
@ -447,8 +449,8 @@ struct SparseQR_QProduct : ReturnByValue<SparseQR_QProduct<SparseQRType, Derived
}; };
template<typename SparseQRType> template<typename SparseQRType>
struct SparseQRMatrixQReturnType{ struct SparseQRMatrixQReturnType
{
SparseQRMatrixQReturnType(const SparseQRType& qr) : m_qr(qr) {} SparseQRMatrixQReturnType(const SparseQRType& qr) : m_qr(qr) {}
template<typename Derived> template<typename Derived>
SparseQR_QProduct<SparseQRType, Derived> operator*(const MatrixBase<Derived>& other) SparseQR_QProduct<SparseQRType, Derived> operator*(const MatrixBase<Derived>& other)
@ -468,7 +470,8 @@ struct SparseQRMatrixQReturnType{
}; };
template<typename SparseQRType> template<typename SparseQRType>
struct SparseQRMatrixQTransposeReturnType{ struct SparseQRMatrixQTransposeReturnType
{
SparseQRMatrixQTransposeReturnType(const SparseQRType& qr) : m_qr(qr) {} SparseQRMatrixQTransposeReturnType(const SparseQRType& qr) : m_qr(qr) {}
template<typename Derived> template<typename Derived>
SparseQR_QProduct<SparseQRType,Derived> operator*(const MatrixBase<Derived>& other) SparseQR_QProduct<SparseQRType,Derived> operator*(const MatrixBase<Derived>& other)
@ -477,5 +480,7 @@ struct SparseQRMatrixQTransposeReturnType{
} }
const SparseQRType& m_qr; const SparseQRType& m_qr;
}; };
} // end namespace Eigen } // end namespace Eigen
#endif
#endif

2
doc/Manual.dox
View File

@ -131,6 +131,8 @@ namespace Eigen {
\ingroup Sparse_Reference */ \ingroup Sparse_Reference */
/** \addtogroup SparseLU_Module /** \addtogroup SparseLU_Module
\ingroup Sparse_Reference */ \ingroup Sparse_Reference */
/** \addtogroup SparseQR_Module
\ingroup Sparse_Reference */
/** \addtogroup IterativeLinearSolvers_Module /** \addtogroup IterativeLinearSolvers_Module
\ingroup Sparse_Reference */ \ingroup Sparse_Reference */
/** \addtogroup Support_modules /** \addtogroup Support_modules