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