$treeview $search $mathjax
Eigen
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2012 Desire Nuentsa <desire.nuentsa_wakam@inria.fr> 00005 // 00006 // This Source Code Form is subject to the terms of the Mozilla 00007 // Public License v. 2.0. If a copy of the MPL was not distributed 00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00009 00010 #ifndef EIGEN_SUITESPARSEQRSUPPORT_H 00011 #define EIGEN_SUITESPARSEQRSUPPORT_H 00012 00013 namespace Eigen { 00014 00015 template<typename MatrixType> class SPQR; 00016 template<typename SPQRType> struct SPQRMatrixQReturnType; 00017 template<typename SPQRType> struct SPQRMatrixQTransposeReturnType; 00018 template <typename SPQRType, typename Derived> struct SPQR_QProduct; 00019 namespace internal { 00020 template <typename SPQRType> struct traits<SPQRMatrixQReturnType<SPQRType> > 00021 { 00022 typedef typename SPQRType::MatrixType ReturnType; 00023 }; 00024 template <typename SPQRType> struct traits<SPQRMatrixQTransposeReturnType<SPQRType> > 00025 { 00026 typedef typename SPQRType::MatrixType ReturnType; 00027 }; 00028 template <typename SPQRType, typename Derived> struct traits<SPQR_QProduct<SPQRType, Derived> > 00029 { 00030 typedef typename Derived::PlainObject ReturnType; 00031 }; 00032 } // End namespace internal 00033 00056 template<typename _MatrixType> 00057 class SPQR 00058 { 00059 public: 00060 typedef typename _MatrixType::Scalar Scalar; 00061 typedef typename _MatrixType::RealScalar RealScalar; 00062 typedef UF_long Index ; 00063 typedef SparseMatrix<Scalar, ColMajor, Index> MatrixType; 00064 typedef PermutationMatrix<Dynamic, Dynamic> PermutationType; 00065 public: 00066 SPQR() 00067 : m_isInitialized(false), m_ordering(SPQR_ORDERING_DEFAULT), m_allow_tol(SPQR_DEFAULT_TOL), m_tolerance (NumTraits<Scalar>::epsilon()), m_useDefaultThreshold(true) 00068 { 00069 cholmod_l_start(&m_cc); 00070 } 00071 00072 SPQR(const _MatrixType& matrix) 00073 : m_isInitialized(false), m_ordering(SPQR_ORDERING_DEFAULT), m_allow_tol(SPQR_DEFAULT_TOL), m_tolerance (NumTraits<Scalar>::epsilon()), m_useDefaultThreshold(true) 00074 { 00075 cholmod_l_start(&m_cc); 00076 compute(matrix); 00077 } 00078 00079 ~SPQR() 00080 { 00081 SPQR_free(); 00082 cholmod_l_finish(&m_cc); 00083 } 00084 void SPQR_free() 00085 { 00086 cholmod_l_free_sparse(&m_H, &m_cc); 00087 cholmod_l_free_sparse(&m_cR, &m_cc); 00088 cholmod_l_free_dense(&m_HTau, &m_cc); 00089 std::free(m_E); 00090 std::free(m_HPinv); 00091 } 00092 00093 void compute(const _MatrixType& matrix) 00094 { 00095 if(m_isInitialized) SPQR_free(); 00096 00097 MatrixType mat(matrix); 00098 00099 /* Compute the default threshold as in MatLab, see: 00100 * Tim Davis, "Algorithm 915, SuiteSparseQR: Multifrontal Multithreaded Rank-Revealing 00101 * Sparse QR Factorization, ACM Trans. on Math. Soft. 38(1), 2011, Page 8:3 00102 */ 00103 RealScalar pivotThreshold = m_tolerance; 00104 if(m_useDefaultThreshold) 00105 { 00106 using std::max; 00107 RealScalar max2Norm = 0.0; 00108 for (int j = 0; j < mat.cols(); j++) max2Norm = (max)(max2Norm, mat.col(j).norm()); 00109 if(max2Norm==RealScalar(0)) 00110 max2Norm = RealScalar(1); 00111 pivotThreshold = 20 * (mat.rows() + mat.cols()) * max2Norm * NumTraits<RealScalar>::epsilon(); 00112 } 00113 00114 cholmod_sparse A; 00115 A = viewAsCholmod(mat); 00116 Index col = matrix.cols(); 00117 m_rank = SuiteSparseQR<Scalar>(m_ordering, pivotThreshold, col, &A, 00118 &m_cR, &m_E, &m_H, &m_HPinv, &m_HTau, &m_cc); 00119 00120 if (!m_cR) 00121 { 00122 m_info = NumericalIssue; 00123 m_isInitialized = false; 00124 return; 00125 } 00126 m_info = Success; 00127 m_isInitialized = true; 00128 m_isRUpToDate = false; 00129 } 00133 inline Index rows() const {return m_cR->nrow; } 00134 00138 inline Index cols() const { return m_cR->ncol; } 00139 00144 template<typename Rhs> 00145 inline const internal::solve_retval<SPQR, Rhs> solve(const MatrixBase<Rhs>& B) const 00146 { 00147 eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()"); 00148 eigen_assert(this->rows()==B.rows() 00149 && "SPQR::solve(): invalid number of rows of the right hand side matrix B"); 00150 return internal::solve_retval<SPQR, Rhs>(*this, B.derived()); 00151 } 00152 00153 template<typename Rhs, typename Dest> 00154 void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const 00155 { 00156 eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()"); 00157 eigen_assert(b.cols()==1 && "This method is for vectors only"); 00158 00159 //Compute Q^T * b 00160 typename Dest::PlainObject y, y2; 00161 y = matrixQ().transpose() * b; 00162 00163 // Solves with the triangular matrix R 00164 Index rk = this->rank(); 00165 y2 = y; 00166 y.resize((std::max)(cols(),Index(y.rows())),y.cols()); 00167 y.topRows(rk) = this->matrixR().topLeftCorner(rk, rk).template triangularView<Upper>().solve(y2.topRows(rk)); 00168 00169 // Apply the column permutation 00170 // colsPermutation() performs a copy of the permutation, 00171 // so let's apply it manually: 00172 for(Index i = 0; i < rk; ++i) dest.row(m_E[i]) = y.row(i); 00173 for(Index i = rk; i < cols(); ++i) dest.row(m_E[i]).setZero(); 00174 00175 // y.bottomRows(y.rows()-rk).setZero(); 00176 // dest = colsPermutation() * y.topRows(cols()); 00177 00178 m_info = Success; 00179 } 00180 00183 const MatrixType matrixR() const 00184 { 00185 eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()"); 00186 if(!m_isRUpToDate) { 00187 m_R = viewAsEigen<Scalar,ColMajor, typename MatrixType::Index>(*m_cR); 00188 m_isRUpToDate = true; 00189 } 00190 return m_R; 00191 } 00193 SPQRMatrixQReturnType<SPQR> matrixQ() const 00194 { 00195 return SPQRMatrixQReturnType<SPQR>(*this); 00196 } 00198 PermutationType colsPermutation() const 00199 { 00200 eigen_assert(m_isInitialized && "Decomposition is not initialized."); 00201 Index n = m_cR->ncol; 00202 PermutationType colsPerm(n); 00203 for(Index j = 0; j <n; j++) colsPerm.indices()(j) = m_E[j]; 00204 return colsPerm; 00205 00206 } 00211 Index rank() const 00212 { 00213 eigen_assert(m_isInitialized && "Decomposition is not initialized."); 00214 return m_cc.SPQR_istat[4]; 00215 } 00217 void setSPQROrdering(int ord) { m_ordering = ord;} 00219 void setPivotThreshold(const RealScalar& tol) 00220 { 00221 m_useDefaultThreshold = false; 00222 m_tolerance = tol; 00223 } 00224 00226 cholmod_common *cholmodCommon() const { return &m_cc; } 00227 00228 00234 ComputationInfo info() const 00235 { 00236 eigen_assert(m_isInitialized && "Decomposition is not initialized."); 00237 return m_info; 00238 } 00239 protected: 00240 bool m_isInitialized; 00241 bool m_analysisIsOk; 00242 bool m_factorizationIsOk; 00243 mutable bool m_isRUpToDate; 00244 mutable ComputationInfo m_info; 00245 int m_ordering; // Ordering method to use, see SPQR's manual 00246 int m_allow_tol; // Allow to use some tolerance during numerical factorization. 00247 RealScalar m_tolerance; // treat columns with 2-norm below this tolerance as zero 00248 mutable cholmod_sparse *m_cR; // The sparse R factor in cholmod format 00249 mutable MatrixType m_R; // The sparse matrix R in Eigen format 00250 mutable Index *m_E; // The permutation applied to columns 00251 mutable cholmod_sparse *m_H; //The householder vectors 00252 mutable Index *m_HPinv; // The row permutation of H 00253 mutable cholmod_dense *m_HTau; // The Householder coefficients 00254 mutable Index m_rank; // The rank of the matrix 00255 mutable cholmod_common m_cc; // Workspace and parameters 00256 bool m_useDefaultThreshold; // Use default threshold 00257 template<typename ,typename > friend struct SPQR_QProduct; 00258 }; 00259 00260 template <typename SPQRType, typename Derived> 00261 struct SPQR_QProduct : ReturnByValue<SPQR_QProduct<SPQRType,Derived> > 00262 { 00263 typedef typename SPQRType::Scalar Scalar; 00264 typedef typename SPQRType::Index Index; 00265 //Define the constructor to get reference to argument types 00266 SPQR_QProduct(const SPQRType& spqr, const Derived& other, bool transpose) : m_spqr(spqr),m_other(other),m_transpose(transpose) {} 00267 00268 inline Index rows() const { return m_transpose ? m_spqr.rows() : m_spqr.cols(); } 00269 inline Index cols() const { return m_other.cols(); } 00270 // Assign to a vector 00271 template<typename ResType> 00272 void evalTo(ResType& res) const 00273 { 00274 cholmod_dense y_cd; 00275 cholmod_dense *x_cd; 00276 int method = m_transpose ? SPQR_QTX : SPQR_QX; 00277 cholmod_common *cc = m_spqr.cholmodCommon(); 00278 y_cd = viewAsCholmod(m_other.const_cast_derived()); 00279 x_cd = SuiteSparseQR_qmult<Scalar>(method, m_spqr.m_H, m_spqr.m_HTau, m_spqr.m_HPinv, &y_cd, cc); 00280 res = Matrix<Scalar,ResType::RowsAtCompileTime,ResType::ColsAtCompileTime>::Map(reinterpret_cast<Scalar*>(x_cd->x), x_cd->nrow, x_cd->ncol); 00281 cholmod_l_free_dense(&x_cd, cc); 00282 } 00283 const SPQRType& m_spqr; 00284 const Derived& m_other; 00285 bool m_transpose; 00286 00287 }; 00288 template<typename SPQRType> 00289 struct SPQRMatrixQReturnType{ 00290 00291 SPQRMatrixQReturnType(const SPQRType& spqr) : m_spqr(spqr) {} 00292 template<typename Derived> 00293 SPQR_QProduct<SPQRType, Derived> operator*(const MatrixBase<Derived>& other) 00294 { 00295 return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(),false); 00296 } 00297 SPQRMatrixQTransposeReturnType<SPQRType> adjoint() const 00298 { 00299 return SPQRMatrixQTransposeReturnType<SPQRType>(m_spqr); 00300 } 00301 // To use for operations with the transpose of Q 00302 SPQRMatrixQTransposeReturnType<SPQRType> transpose() const 00303 { 00304 return SPQRMatrixQTransposeReturnType<SPQRType>(m_spqr); 00305 } 00306 const SPQRType& m_spqr; 00307 }; 00308 00309 template<typename SPQRType> 00310 struct SPQRMatrixQTransposeReturnType{ 00311 SPQRMatrixQTransposeReturnType(const SPQRType& spqr) : m_spqr(spqr) {} 00312 template<typename Derived> 00313 SPQR_QProduct<SPQRType,Derived> operator*(const MatrixBase<Derived>& other) 00314 { 00315 return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(), true); 00316 } 00317 const SPQRType& m_spqr; 00318 }; 00319 00320 namespace internal { 00321 00322 template<typename _MatrixType, typename Rhs> 00323 struct solve_retval<SPQR<_MatrixType>, Rhs> 00324 : solve_retval_base<SPQR<_MatrixType>, Rhs> 00325 { 00326 typedef SPQR<_MatrixType> Dec; 00327 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) 00328 00329 template<typename Dest> void evalTo(Dest& dst) const 00330 { 00331 dec()._solve(rhs(),dst); 00332 } 00333 }; 00334 00335 } // end namespace internal 00336 00337 }// End namespace Eigen 00338 #endif