Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 1 | // This file is part of Eigen, a lightweight C++ template library |
| 2 | // for linear algebra. |
| 3 | // |
| 4 | // Copyright (C) 2012 Desire Nuentsa <desire.nuentsa_wakam@inria.fr> |
| 5 | // |
| 6 | // This Source Code Form is subject to the terms of the Mozilla |
| 7 | // Public License v. 2.0. If a copy of the MPL was not distributed |
| 8 | // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. |
| 9 | |
| 10 | #ifndef EIGEN_SUITESPARSEQRSUPPORT_H |
| 11 | #define EIGEN_SUITESPARSEQRSUPPORT_H |
| 12 | |
| 13 | namespace Eigen { |
| 14 | |
| 15 | template<typename MatrixType> class SPQR; |
| 16 | template<typename SPQRType> struct SPQRMatrixQReturnType; |
| 17 | template<typename SPQRType> struct SPQRMatrixQTransposeReturnType; |
| 18 | template <typename SPQRType, typename Derived> struct SPQR_QProduct; |
| 19 | namespace internal { |
| 20 | template <typename SPQRType> struct traits<SPQRMatrixQReturnType<SPQRType> > |
| 21 | { |
| 22 | typedef typename SPQRType::MatrixType ReturnType; |
| 23 | }; |
| 24 | template <typename SPQRType> struct traits<SPQRMatrixQTransposeReturnType<SPQRType> > |
| 25 | { |
| 26 | typedef typename SPQRType::MatrixType ReturnType; |
| 27 | }; |
| 28 | template <typename SPQRType, typename Derived> struct traits<SPQR_QProduct<SPQRType, Derived> > |
| 29 | { |
| 30 | typedef typename Derived::PlainObject ReturnType; |
| 31 | }; |
| 32 | } // End namespace internal |
| 33 | |
| 34 | /** |
| 35 | * \ingroup SPQRSupport_Module |
| 36 | * \class SPQR |
| 37 | * \brief Sparse QR factorization based on SuiteSparseQR library |
| 38 | * |
| 39 | * This class is used to perform a multithreaded and multifrontal rank-revealing QR decomposition |
| 40 | * of sparse matrices. The result is then used to solve linear leasts_square systems. |
| 41 | * Clearly, a QR factorization is returned such that A*P = Q*R where : |
| 42 | * |
| 43 | * P is the column permutation. Use colsPermutation() to get it. |
| 44 | * |
| 45 | * Q is the orthogonal matrix represented as Householder reflectors. |
| 46 | * Use matrixQ() to get an expression and matrixQ().transpose() to get the transpose. |
| 47 | * You can then apply it to a vector. |
| 48 | * |
| 49 | * R is the sparse triangular factor. Use matrixQR() to get it as SparseMatrix. |
| 50 | * NOTE : The Index type of R is always UF_long. You can get it with SPQR::Index |
| 51 | * |
| 52 | * \tparam _MatrixType The type of the sparse matrix A, must be a column-major SparseMatrix<> |
| 53 | * NOTE |
| 54 | * |
| 55 | */ |
| 56 | template<typename _MatrixType> |
| 57 | class SPQR |
| 58 | { |
| 59 | public: |
| 60 | typedef typename _MatrixType::Scalar Scalar; |
| 61 | typedef typename _MatrixType::RealScalar RealScalar; |
| 62 | typedef UF_long Index ; |
| 63 | typedef SparseMatrix<Scalar, ColMajor, Index> MatrixType; |
| 64 | typedef PermutationMatrix<Dynamic, Dynamic> PermutationType; |
| 65 | public: |
| 66 | SPQR() |
| 67 | : m_isInitialized(false), |
| 68 | m_ordering(SPQR_ORDERING_DEFAULT), |
| 69 | m_allow_tol(SPQR_DEFAULT_TOL), |
| 70 | m_tolerance (NumTraits<Scalar>::epsilon()) |
| 71 | { |
| 72 | cholmod_l_start(&m_cc); |
| 73 | } |
| 74 | |
| 75 | SPQR(const _MatrixType& matrix) |
| 76 | : m_isInitialized(false), |
| 77 | m_ordering(SPQR_ORDERING_DEFAULT), |
| 78 | m_allow_tol(SPQR_DEFAULT_TOL), |
| 79 | m_tolerance (NumTraits<Scalar>::epsilon()) |
| 80 | { |
| 81 | cholmod_l_start(&m_cc); |
| 82 | compute(matrix); |
| 83 | } |
| 84 | |
| 85 | ~SPQR() |
| 86 | { |
| 87 | SPQR_free(); |
| 88 | cholmod_l_finish(&m_cc); |
| 89 | } |
| 90 | void SPQR_free() |
| 91 | { |
| 92 | cholmod_l_free_sparse(&m_H, &m_cc); |
| 93 | cholmod_l_free_sparse(&m_cR, &m_cc); |
| 94 | cholmod_l_free_dense(&m_HTau, &m_cc); |
| 95 | std::free(m_E); |
| 96 | std::free(m_HPinv); |
| 97 | } |
| 98 | |
| 99 | void compute(const _MatrixType& matrix) |
| 100 | { |
| 101 | if(m_isInitialized) SPQR_free(); |
| 102 | |
| 103 | MatrixType mat(matrix); |
| 104 | cholmod_sparse A; |
| 105 | A = viewAsCholmod(mat); |
| 106 | Index col = matrix.cols(); |
| 107 | m_rank = SuiteSparseQR<Scalar>(m_ordering, m_tolerance, col, &A, |
| 108 | &m_cR, &m_E, &m_H, &m_HPinv, &m_HTau, &m_cc); |
| 109 | |
| 110 | if (!m_cR) |
| 111 | { |
| 112 | m_info = NumericalIssue; |
| 113 | m_isInitialized = false; |
| 114 | return; |
| 115 | } |
| 116 | m_info = Success; |
| 117 | m_isInitialized = true; |
| 118 | m_isRUpToDate = false; |
| 119 | } |
| 120 | /** |
| 121 | * Get the number of rows of the input matrix and the Q matrix |
| 122 | */ |
| 123 | inline Index rows() const {return m_H->nrow; } |
| 124 | |
| 125 | /** |
| 126 | * Get the number of columns of the input matrix. |
| 127 | */ |
| 128 | inline Index cols() const { return m_cR->ncol; } |
| 129 | |
| 130 | /** \returns the solution X of \f$ A X = B \f$ using the current decomposition of A. |
| 131 | * |
| 132 | * \sa compute() |
| 133 | */ |
| 134 | template<typename Rhs> |
| 135 | inline const internal::solve_retval<SPQR, Rhs> solve(const MatrixBase<Rhs>& B) const |
| 136 | { |
| 137 | eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()"); |
| 138 | eigen_assert(this->rows()==B.rows() |
| 139 | && "SPQR::solve(): invalid number of rows of the right hand side matrix B"); |
| 140 | return internal::solve_retval<SPQR, Rhs>(*this, B.derived()); |
| 141 | } |
| 142 | |
| 143 | template<typename Rhs, typename Dest> |
| 144 | void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const |
| 145 | { |
| 146 | eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()"); |
| 147 | eigen_assert(b.cols()==1 && "This method is for vectors only"); |
| 148 | |
| 149 | //Compute Q^T * b |
| 150 | typename Dest::PlainObject y; |
| 151 | y = matrixQ().transpose() * b; |
| 152 | // Solves with the triangular matrix R |
| 153 | Index rk = this->rank(); |
| 154 | y.topRows(rk) = this->matrixR().topLeftCorner(rk, rk).template triangularView<Upper>().solve(y.topRows(rk)); |
| 155 | y.bottomRows(cols()-rk).setZero(); |
| 156 | // Apply the column permutation |
| 157 | dest.topRows(cols()) = colsPermutation() * y.topRows(cols()); |
| 158 | |
| 159 | m_info = Success; |
| 160 | } |
| 161 | |
| 162 | /** \returns the sparse triangular factor R. It is a sparse matrix |
| 163 | */ |
| 164 | const MatrixType matrixR() const |
| 165 | { |
| 166 | eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()"); |
| 167 | if(!m_isRUpToDate) { |
| 168 | m_R = viewAsEigen<Scalar,ColMajor, typename MatrixType::Index>(*m_cR); |
| 169 | m_isRUpToDate = true; |
| 170 | } |
| 171 | return m_R; |
| 172 | } |
| 173 | /// Get an expression of the matrix Q |
| 174 | SPQRMatrixQReturnType<SPQR> matrixQ() const |
| 175 | { |
| 176 | return SPQRMatrixQReturnType<SPQR>(*this); |
| 177 | } |
| 178 | /// Get the permutation that was applied to columns of A |
| 179 | PermutationType colsPermutation() const |
| 180 | { |
| 181 | eigen_assert(m_isInitialized && "Decomposition is not initialized."); |
| 182 | Index n = m_cR->ncol; |
| 183 | PermutationType colsPerm(n); |
| 184 | for(Index j = 0; j <n; j++) colsPerm.indices()(j) = m_E[j]; |
| 185 | return colsPerm; |
| 186 | |
| 187 | } |
| 188 | /** |
| 189 | * Gets the rank of the matrix. |
| 190 | * It should be equal to matrixQR().cols if the matrix is full-rank |
| 191 | */ |
| 192 | Index rank() const |
| 193 | { |
| 194 | eigen_assert(m_isInitialized && "Decomposition is not initialized."); |
| 195 | return m_cc.SPQR_istat[4]; |
| 196 | } |
| 197 | /// Set the fill-reducing ordering method to be used |
| 198 | void setSPQROrdering(int ord) { m_ordering = ord;} |
| 199 | /// Set the tolerance tol to treat columns with 2-norm < =tol as zero |
| 200 | void setPivotThreshold(const RealScalar& tol) { m_tolerance = tol; } |
| 201 | |
| 202 | /** \returns a pointer to the SPQR workspace */ |
| 203 | cholmod_common *cholmodCommon() const { return &m_cc; } |
| 204 | |
| 205 | |
| 206 | /** \brief Reports whether previous computation was successful. |
| 207 | * |
| 208 | * \returns \c Success if computation was succesful, |
| 209 | * \c NumericalIssue if the sparse QR can not be computed |
| 210 | */ |
| 211 | ComputationInfo info() const |
| 212 | { |
| 213 | eigen_assert(m_isInitialized && "Decomposition is not initialized."); |
| 214 | return m_info; |
| 215 | } |
| 216 | protected: |
| 217 | bool m_isInitialized; |
| 218 | bool m_analysisIsOk; |
| 219 | bool m_factorizationIsOk; |
| 220 | mutable bool m_isRUpToDate; |
| 221 | mutable ComputationInfo m_info; |
| 222 | int m_ordering; // Ordering method to use, see SPQR's manual |
| 223 | int m_allow_tol; // Allow to use some tolerance during numerical factorization. |
| 224 | RealScalar m_tolerance; // treat columns with 2-norm below this tolerance as zero |
| 225 | mutable cholmod_sparse *m_cR; // The sparse R factor in cholmod format |
| 226 | mutable MatrixType m_R; // The sparse matrix R in Eigen format |
| 227 | mutable Index *m_E; // The permutation applied to columns |
| 228 | mutable cholmod_sparse *m_H; //The householder vectors |
| 229 | mutable Index *m_HPinv; // The row permutation of H |
| 230 | mutable cholmod_dense *m_HTau; // The Householder coefficients |
| 231 | mutable Index m_rank; // The rank of the matrix |
| 232 | mutable cholmod_common m_cc; // Workspace and parameters |
| 233 | template<typename ,typename > friend struct SPQR_QProduct; |
| 234 | }; |
| 235 | |
| 236 | template <typename SPQRType, typename Derived> |
| 237 | struct SPQR_QProduct : ReturnByValue<SPQR_QProduct<SPQRType,Derived> > |
| 238 | { |
| 239 | typedef typename SPQRType::Scalar Scalar; |
| 240 | typedef typename SPQRType::Index Index; |
| 241 | //Define the constructor to get reference to argument types |
| 242 | SPQR_QProduct(const SPQRType& spqr, const Derived& other, bool transpose) : m_spqr(spqr),m_other(other),m_transpose(transpose) {} |
| 243 | |
| 244 | inline Index rows() const { return m_transpose ? m_spqr.rows() : m_spqr.cols(); } |
| 245 | inline Index cols() const { return m_other.cols(); } |
| 246 | // Assign to a vector |
| 247 | template<typename ResType> |
| 248 | void evalTo(ResType& res) const |
| 249 | { |
| 250 | cholmod_dense y_cd; |
| 251 | cholmod_dense *x_cd; |
| 252 | int method = m_transpose ? SPQR_QTX : SPQR_QX; |
| 253 | cholmod_common *cc = m_spqr.cholmodCommon(); |
| 254 | y_cd = viewAsCholmod(m_other.const_cast_derived()); |
| 255 | x_cd = SuiteSparseQR_qmult<Scalar>(method, m_spqr.m_H, m_spqr.m_HTau, m_spqr.m_HPinv, &y_cd, cc); |
| 256 | res = Matrix<Scalar,ResType::RowsAtCompileTime,ResType::ColsAtCompileTime>::Map(reinterpret_cast<Scalar*>(x_cd->x), x_cd->nrow, x_cd->ncol); |
| 257 | cholmod_l_free_dense(&x_cd, cc); |
| 258 | } |
| 259 | const SPQRType& m_spqr; |
| 260 | const Derived& m_other; |
| 261 | bool m_transpose; |
| 262 | |
| 263 | }; |
| 264 | template<typename SPQRType> |
| 265 | struct SPQRMatrixQReturnType{ |
| 266 | |
| 267 | SPQRMatrixQReturnType(const SPQRType& spqr) : m_spqr(spqr) {} |
| 268 | template<typename Derived> |
| 269 | SPQR_QProduct<SPQRType, Derived> operator*(const MatrixBase<Derived>& other) |
| 270 | { |
| 271 | return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(),false); |
| 272 | } |
| 273 | SPQRMatrixQTransposeReturnType<SPQRType> adjoint() const |
| 274 | { |
| 275 | return SPQRMatrixQTransposeReturnType<SPQRType>(m_spqr); |
| 276 | } |
| 277 | // To use for operations with the transpose of Q |
| 278 | SPQRMatrixQTransposeReturnType<SPQRType> transpose() const |
| 279 | { |
| 280 | return SPQRMatrixQTransposeReturnType<SPQRType>(m_spqr); |
| 281 | } |
| 282 | const SPQRType& m_spqr; |
| 283 | }; |
| 284 | |
| 285 | template<typename SPQRType> |
| 286 | struct SPQRMatrixQTransposeReturnType{ |
| 287 | SPQRMatrixQTransposeReturnType(const SPQRType& spqr) : m_spqr(spqr) {} |
| 288 | template<typename Derived> |
| 289 | SPQR_QProduct<SPQRType,Derived> operator*(const MatrixBase<Derived>& other) |
| 290 | { |
| 291 | return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(), true); |
| 292 | } |
| 293 | const SPQRType& m_spqr; |
| 294 | }; |
| 295 | |
| 296 | namespace internal { |
| 297 | |
| 298 | template<typename _MatrixType, typename Rhs> |
| 299 | struct solve_retval<SPQR<_MatrixType>, Rhs> |
| 300 | : solve_retval_base<SPQR<_MatrixType>, Rhs> |
| 301 | { |
| 302 | typedef SPQR<_MatrixType> Dec; |
| 303 | EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) |
| 304 | |
| 305 | template<typename Dest> void evalTo(Dest& dst) const |
| 306 | { |
| 307 | dec()._solve(rhs(),dst); |
| 308 | } |
| 309 | }; |
| 310 | |
| 311 | } // end namespace internal |
| 312 | |
| 313 | }// End namespace Eigen |
| 314 | #endif |