1f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org// This file is part of Eigen, a lightweight C++ template library 2f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org// for linear algebra. 3f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org// 4f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org// Copyright (C) 2012 Desire Nuentsa <desire.nuentsa_wakam@inria.fr> 5f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org// 6f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org// This Source Code Form is subject to the terms of the Mozilla 7f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org// Public License v. 2.0. If a copy of the MPL was not distributed 8f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org 10f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org#ifndef EIGEN_SUITESPARSEQRSUPPORT_H 11f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org#define EIGEN_SUITESPARSEQRSUPPORT_H 12f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org 13f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.orgnamespace Eigen { 14f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org 15f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org template<typename MatrixType> class SPQR; 16f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org template<typename SPQRType> struct SPQRMatrixQReturnType; 17f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org template<typename SPQRType> struct SPQRMatrixQTransposeReturnType; 18f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org template <typename SPQRType, typename Derived> struct SPQR_QProduct; 19f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org namespace internal { 20f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org template <typename SPQRType> struct traits<SPQRMatrixQReturnType<SPQRType> > 21f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org { 22f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org typedef typename SPQRType::MatrixType ReturnType; 23f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org }; 24f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org template <typename SPQRType> struct traits<SPQRMatrixQTransposeReturnType<SPQRType> > 25f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org { 26f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org typedef typename SPQRType::MatrixType ReturnType; 27f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org }; 28f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org template <typename SPQRType, typename Derived> struct traits<SPQR_QProduct<SPQRType, Derived> > 29f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org { 30f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org typedef typename Derived::PlainObject ReturnType; 31f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org }; 32f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org } // End namespace internal 33f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org 34f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org/** 35f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org * \ingroup SPQRSupport_Module 36f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org * \class SPQR 37f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org * \brief Sparse QR factorization based on SuiteSparseQR library 38f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org * 39f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org * This class is used to perform a multithreaded and multifrontal rank-revealing QR decomposition 40f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org * of sparse matrices. The result is then used to solve linear leasts_square systems. 41f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org * Clearly, a QR factorization is returned such that A*P = Q*R where : 42f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org * 43f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org * P is the column permutation. Use colsPermutation() to get it. 44f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org * 45f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org * Q is the orthogonal matrix represented as Householder reflectors. 46f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org * Use matrixQ() to get an expression and matrixQ().transpose() to get the transpose. 47f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org * You can then apply it to a vector. 48f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org * 49f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org * R is the sparse triangular factor. Use matrixQR() to get it as SparseMatrix. 50f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org * NOTE : The Index type of R is always UF_long. You can get it with SPQR::Index 51f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org * 52f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org * \tparam _MatrixType The type of the sparse matrix A, must be a column-major SparseMatrix<> 53f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org * NOTE 54f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org * 55f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org */ 56f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.orgtemplate<typename _MatrixType> 57f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.orgclass SPQR 58f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org{ 59f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org public: 60f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org typedef typename _MatrixType::Scalar Scalar; 61f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org typedef typename _MatrixType::RealScalar RealScalar; 62f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org typedef UF_long Index ; 63f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org typedef SparseMatrix<Scalar, ColMajor, Index> MatrixType; 64f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org typedef PermutationMatrix<Dynamic, Dynamic> PermutationType; 65f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org public: 66f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org SPQR() 67f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org : m_isInitialized(false), 68f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org m_ordering(SPQR_ORDERING_DEFAULT), 69f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org m_allow_tol(SPQR_DEFAULT_TOL), 70f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org m_tolerance (NumTraits<Scalar>::epsilon()) 71f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org { 72f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org cholmod_l_start(&m_cc); 73f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org } 74f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org 75f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org SPQR(const _MatrixType& matrix) 76f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org : m_isInitialized(false), 77f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org m_ordering(SPQR_ORDERING_DEFAULT), 78f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org m_allow_tol(SPQR_DEFAULT_TOL), 79f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org m_tolerance (NumTraits<Scalar>::epsilon()) 80f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org { 81f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org cholmod_l_start(&m_cc); 82f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org compute(matrix); 83f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org } 84f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org 85f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org ~SPQR() 86f2ba7591b1407a7ee9209f842c50696914dc2dedkbr@chromium.org { 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 236template <typename SPQRType, typename Derived> 237struct 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}; 264template<typename SPQRType> 265struct 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 285template<typename SPQRType> 286struct 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 296namespace internal { 297 298template<typename _MatrixType, typename Rhs> 299struct 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 315