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