diff options
Diffstat (limited to 'eigen/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h')
| -rw-r--r-- | eigen/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h | 123 |
1 files changed, 49 insertions, 74 deletions
diff --git a/eigen/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h b/eigen/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h index 3613810..953d57c 100644 --- a/eigen/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h +++ b/eigen/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h @@ -2,6 +2,7 @@ // for linear algebra. // // Copyright (C) 2012 Desire Nuentsa <desire.nuentsa_wakam@inria.fr> +// Copyright (C) 2014 Gael Guennebaud <gael.guennebaud@inria.fr> // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed @@ -32,45 +33,54 @@ namespace Eigen { } // End namespace internal /** - * \ingroup SPQRSupport_Module - * \class SPQR - * \brief Sparse QR factorization based on SuiteSparseQR library - * - * This class is used to perform a multithreaded and multifrontal rank-revealing 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 : The Index type of R is always SuiteSparse_long. You can get it with SPQR::Index - * - * \tparam _MatrixType The type of the sparse matrix A, must be a column-major SparseMatrix<> - * NOTE - * - */ + * \ingroup SPQRSupport_Module + * \class SPQR + * \brief Sparse QR factorization based on SuiteSparseQR library + * + * This class is used to perform a multithreaded and multifrontal rank-revealing 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 : The Index type of R is always SuiteSparse_long. You can get it with SPQR::Index + * + * \tparam _MatrixType The type of the sparse matrix A, must be a column-major SparseMatrix<> + * + * \implsparsesolverconcept + * + * + */ template<typename _MatrixType> -class SPQR +class SPQR : public SparseSolverBase<SPQR<_MatrixType> > { + protected: + typedef SparseSolverBase<SPQR<_MatrixType> > Base; + using Base::m_isInitialized; public: typedef typename _MatrixType::Scalar Scalar; typedef typename _MatrixType::RealScalar RealScalar; - typedef SuiteSparse_long Index ; - typedef SparseMatrix<Scalar, ColMajor, Index> MatrixType; - typedef PermutationMatrix<Dynamic, Dynamic> PermutationType; + typedef SuiteSparse_long StorageIndex ; + typedef SparseMatrix<Scalar, ColMajor, StorageIndex> MatrixType; + typedef Map<PermutationMatrix<Dynamic, Dynamic, StorageIndex> > PermutationType; + enum { + ColsAtCompileTime = Dynamic, + MaxColsAtCompileTime = Dynamic + }; public: SPQR() - : m_isInitialized(false), m_ordering(SPQR_ORDERING_DEFAULT), m_allow_tol(SPQR_DEFAULT_TOL), m_tolerance (NumTraits<Scalar>::epsilon()), m_useDefaultThreshold(true) + : m_ordering(SPQR_ORDERING_DEFAULT), m_allow_tol(SPQR_DEFAULT_TOL), m_tolerance (NumTraits<Scalar>::epsilon()), m_useDefaultThreshold(true) { cholmod_l_start(&m_cc); } - SPQR(const _MatrixType& matrix) - : m_isInitialized(false), m_ordering(SPQR_ORDERING_DEFAULT), m_allow_tol(SPQR_DEFAULT_TOL), m_tolerance (NumTraits<Scalar>::epsilon()), m_useDefaultThreshold(true) + explicit SPQR(const _MatrixType& matrix) + : m_ordering(SPQR_ORDERING_DEFAULT), m_allow_tol(SPQR_DEFAULT_TOL), m_tolerance (NumTraits<Scalar>::epsilon()), m_useDefaultThreshold(true) { cholmod_l_start(&m_cc); compute(matrix); @@ -103,23 +113,22 @@ class SPQR RealScalar pivotThreshold = m_tolerance; if(m_useDefaultThreshold) { - using std::max; RealScalar max2Norm = 0.0; - for (int j = 0; j < mat.cols(); j++) max2Norm = (max)(max2Norm, mat.col(j).norm()); + for (int j = 0; j < mat.cols(); j++) max2Norm = numext::maxi(max2Norm, mat.col(j).norm()); if(max2Norm==RealScalar(0)) max2Norm = RealScalar(1); pivotThreshold = 20 * (mat.rows() + mat.cols()) * max2Norm * NumTraits<RealScalar>::epsilon(); } - cholmod_sparse A; A = viewAsCholmod(mat); + m_rows = matrix.rows(); Index col = matrix.cols(); m_rank = SuiteSparseQR<Scalar>(m_ordering, pivotThreshold, col, &A, &m_cR, &m_E, &m_H, &m_HPinv, &m_HTau, &m_cc); if (!m_cR) { - m_info = NumericalIssue; + m_info = NumericalIssue; m_isInitialized = false; return; } @@ -130,28 +139,15 @@ class SPQR /** * Get the number of rows of the input matrix and the Q matrix */ - inline Index rows() const {return m_cR->nrow; } + inline Index rows() const {return m_rows; } /** * Get the number of columns of the input matrix. */ inline Index cols() const { return m_cR->ncol; } - - /** \returns the solution X of \f$ A X = B \f$ using the current decomposition of A. - * - * \sa compute() - */ - template<typename Rhs> - inline const internal::solve_retval<SPQR, Rhs> solve(const MatrixBase<Rhs>& B) const - { - eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()"); - eigen_assert(this->rows()==B.rows() - && "SPQR::solve(): invalid number of rows of the right hand side matrix B"); - return internal::solve_retval<SPQR, Rhs>(*this, B.derived()); - } template<typename Rhs, typename Dest> - void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const + void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const { eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()"); eigen_assert(b.cols()==1 && "This method is for vectors only"); @@ -184,7 +180,7 @@ class SPQR { eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()"); if(!m_isRUpToDate) { - m_R = viewAsEigen<Scalar,ColMajor, typename MatrixType::Index>(*m_cR); + m_R = viewAsEigen<Scalar,ColMajor, typename MatrixType::StorageIndex>(*m_cR); m_isRUpToDate = true; } return m_R; @@ -198,11 +194,7 @@ class SPQR PermutationType colsPermutation() const { eigen_assert(m_isInitialized && "Decomposition is not initialized."); - Index n = m_cR->ncol; - PermutationType colsPerm(n); - for(Index j = 0; j <n; j++) colsPerm.indices()(j) = m_E[j]; - return colsPerm; - + return PermutationType(m_E, m_cR->ncol); } /** * Gets the rank of the matrix. @@ -237,7 +229,6 @@ class SPQR return m_info; } protected: - bool m_isInitialized; bool m_analysisIsOk; bool m_factorizationIsOk; mutable bool m_isRUpToDate; @@ -247,13 +238,14 @@ class SPQR RealScalar m_tolerance; // treat columns with 2-norm below this tolerance as zero mutable cholmod_sparse *m_cR; // The sparse R factor in cholmod format mutable MatrixType m_R; // The sparse matrix R in Eigen format - mutable Index *m_E; // The permutation applied to columns + mutable StorageIndex *m_E; // The permutation applied to columns mutable cholmod_sparse *m_H; //The householder vectors - mutable Index *m_HPinv; // The row permutation of H + mutable StorageIndex *m_HPinv; // The row permutation of H mutable cholmod_dense *m_HTau; // The Householder coefficients mutable Index m_rank; // The rank of the matrix mutable cholmod_common m_cc; // Workspace and parameters bool m_useDefaultThreshold; // Use default threshold + Index m_rows; template<typename ,typename > friend struct SPQR_QProduct; }; @@ -261,7 +253,7 @@ template <typename SPQRType, typename Derived> struct SPQR_QProduct : ReturnByValue<SPQR_QProduct<SPQRType,Derived> > { typedef typename SPQRType::Scalar Scalar; - typedef typename SPQRType::Index Index; + typedef typename SPQRType::StorageIndex StorageIndex; //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) {} @@ -317,22 +309,5 @@ struct SPQRMatrixQTransposeReturnType{ const SPQRType& m_spqr; }; -namespace internal { - -template<typename _MatrixType, typename Rhs> -struct solve_retval<SPQR<_MatrixType>, Rhs> - : solve_retval_base<SPQR<_MatrixType>, Rhs> -{ - typedef SPQR<_MatrixType> Dec; - EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) - - template<typename Dest> void evalTo(Dest& dst) const - { - dec()._solve(rhs(),dst); - } -}; - -} // end namespace internal - }// End namespace Eigen #endif |