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 SuiteSparse_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 SuiteSparse_long Index ;
|
---|
63 | typedef SparseMatrix<Scalar, ColMajor, Index> MatrixType;
|
---|
64 | typedef PermutationMatrix<Dynamic, Dynamic> PermutationType;
|
---|
65 | public:
|
---|
66 | SPQR()
|
---|
67 | : m_isInitialized(false), m_ordering(SPQR_ORDERING_DEFAULT), m_allow_tol(SPQR_DEFAULT_TOL), m_tolerance (NumTraits<Scalar>::epsilon()), m_useDefaultThreshold(true)
|
---|
68 | {
|
---|
69 | cholmod_l_start(&m_cc);
|
---|
70 | }
|
---|
71 |
|
---|
72 | SPQR(const _MatrixType& matrix)
|
---|
73 | : m_isInitialized(false), m_ordering(SPQR_ORDERING_DEFAULT), m_allow_tol(SPQR_DEFAULT_TOL), m_tolerance (NumTraits<Scalar>::epsilon()), m_useDefaultThreshold(true)
|
---|
74 | {
|
---|
75 | cholmod_l_start(&m_cc);
|
---|
76 | compute(matrix);
|
---|
77 | }
|
---|
78 |
|
---|
79 | ~SPQR()
|
---|
80 | {
|
---|
81 | SPQR_free();
|
---|
82 | cholmod_l_finish(&m_cc);
|
---|
83 | }
|
---|
84 | void SPQR_free()
|
---|
85 | {
|
---|
86 | cholmod_l_free_sparse(&m_H, &m_cc);
|
---|
87 | cholmod_l_free_sparse(&m_cR, &m_cc);
|
---|
88 | cholmod_l_free_dense(&m_HTau, &m_cc);
|
---|
89 | std::free(m_E);
|
---|
90 | std::free(m_HPinv);
|
---|
91 | }
|
---|
92 |
|
---|
93 | void compute(const _MatrixType& matrix)
|
---|
94 | {
|
---|
95 | if(m_isInitialized) SPQR_free();
|
---|
96 |
|
---|
97 | MatrixType mat(matrix);
|
---|
98 |
|
---|
99 | /* Compute the default threshold as in MatLab, see:
|
---|
100 | * Tim Davis, "Algorithm 915, SuiteSparseQR: Multifrontal Multithreaded Rank-Revealing
|
---|
101 | * Sparse QR Factorization, ACM Trans. on Math. Soft. 38(1), 2011, Page 8:3
|
---|
102 | */
|
---|
103 | RealScalar pivotThreshold = m_tolerance;
|
---|
104 | if(m_useDefaultThreshold)
|
---|
105 | {
|
---|
106 | using std::max;
|
---|
107 | RealScalar max2Norm = 0.0;
|
---|
108 | for (int j = 0; j < mat.cols(); j++) max2Norm = (max)(max2Norm, mat.col(j).norm());
|
---|
109 | if(max2Norm==RealScalar(0))
|
---|
110 | max2Norm = RealScalar(1);
|
---|
111 | pivotThreshold = 20 * (mat.rows() + mat.cols()) * max2Norm * NumTraits<RealScalar>::epsilon();
|
---|
112 | }
|
---|
113 |
|
---|
114 | cholmod_sparse A;
|
---|
115 | A = viewAsCholmod(mat);
|
---|
116 | Index col = matrix.cols();
|
---|
117 | m_rank = SuiteSparseQR<Scalar>(m_ordering, pivotThreshold, col, &A,
|
---|
118 | &m_cR, &m_E, &m_H, &m_HPinv, &m_HTau, &m_cc);
|
---|
119 |
|
---|
120 | if (!m_cR)
|
---|
121 | {
|
---|
122 | m_info = NumericalIssue;
|
---|
123 | m_isInitialized = false;
|
---|
124 | return;
|
---|
125 | }
|
---|
126 | m_info = Success;
|
---|
127 | m_isInitialized = true;
|
---|
128 | m_isRUpToDate = false;
|
---|
129 | }
|
---|
130 | /**
|
---|
131 | * Get the number of rows of the input matrix and the Q matrix
|
---|
132 | */
|
---|
133 | inline Index rows() const {return m_cR->nrow; }
|
---|
134 |
|
---|
135 | /**
|
---|
136 | * Get the number of columns of the input matrix.
|
---|
137 | */
|
---|
138 | inline Index cols() const { return m_cR->ncol; }
|
---|
139 |
|
---|
140 | /** \returns the solution X of \f$ A X = B \f$ using the current decomposition of A.
|
---|
141 | *
|
---|
142 | * \sa compute()
|
---|
143 | */
|
---|
144 | template<typename Rhs>
|
---|
145 | inline const internal::solve_retval<SPQR, Rhs> solve(const MatrixBase<Rhs>& B) const
|
---|
146 | {
|
---|
147 | eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()");
|
---|
148 | eigen_assert(this->rows()==B.rows()
|
---|
149 | && "SPQR::solve(): invalid number of rows of the right hand side matrix B");
|
---|
150 | return internal::solve_retval<SPQR, Rhs>(*this, B.derived());
|
---|
151 | }
|
---|
152 |
|
---|
153 | template<typename Rhs, typename Dest>
|
---|
154 | void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
|
---|
155 | {
|
---|
156 | eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()");
|
---|
157 | eigen_assert(b.cols()==1 && "This method is for vectors only");
|
---|
158 |
|
---|
159 | //Compute Q^T * b
|
---|
160 | typename Dest::PlainObject y, y2;
|
---|
161 | y = matrixQ().transpose() * b;
|
---|
162 |
|
---|
163 | // Solves with the triangular matrix R
|
---|
164 | Index rk = this->rank();
|
---|
165 | y2 = y;
|
---|
166 | y.resize((std::max)(cols(),Index(y.rows())),y.cols());
|
---|
167 | y.topRows(rk) = this->matrixR().topLeftCorner(rk, rk).template triangularView<Upper>().solve(y2.topRows(rk));
|
---|
168 |
|
---|
169 | // Apply the column permutation
|
---|
170 | // colsPermutation() performs a copy of the permutation,
|
---|
171 | // so let's apply it manually:
|
---|
172 | for(Index i = 0; i < rk; ++i) dest.row(m_E[i]) = y.row(i);
|
---|
173 | for(Index i = rk; i < cols(); ++i) dest.row(m_E[i]).setZero();
|
---|
174 |
|
---|
175 | // y.bottomRows(y.rows()-rk).setZero();
|
---|
176 | // dest = colsPermutation() * y.topRows(cols());
|
---|
177 |
|
---|
178 | m_info = Success;
|
---|
179 | }
|
---|
180 |
|
---|
181 | /** \returns the sparse triangular factor R. It is a sparse matrix
|
---|
182 | */
|
---|
183 | const MatrixType matrixR() const
|
---|
184 | {
|
---|
185 | eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()");
|
---|
186 | if(!m_isRUpToDate) {
|
---|
187 | m_R = viewAsEigen<Scalar,ColMajor, typename MatrixType::Index>(*m_cR);
|
---|
188 | m_isRUpToDate = true;
|
---|
189 | }
|
---|
190 | return m_R;
|
---|
191 | }
|
---|
192 | /// Get an expression of the matrix Q
|
---|
193 | SPQRMatrixQReturnType<SPQR> matrixQ() const
|
---|
194 | {
|
---|
195 | return SPQRMatrixQReturnType<SPQR>(*this);
|
---|
196 | }
|
---|
197 | /// Get the permutation that was applied to columns of A
|
---|
198 | PermutationType colsPermutation() const
|
---|
199 | {
|
---|
200 | eigen_assert(m_isInitialized && "Decomposition is not initialized.");
|
---|
201 | Index n = m_cR->ncol;
|
---|
202 | PermutationType colsPerm(n);
|
---|
203 | for(Index j = 0; j <n; j++) colsPerm.indices()(j) = m_E[j];
|
---|
204 | return colsPerm;
|
---|
205 |
|
---|
206 | }
|
---|
207 | /**
|
---|
208 | * Gets the rank of the matrix.
|
---|
209 | * It should be equal to matrixQR().cols if the matrix is full-rank
|
---|
210 | */
|
---|
211 | Index rank() const
|
---|
212 | {
|
---|
213 | eigen_assert(m_isInitialized && "Decomposition is not initialized.");
|
---|
214 | return m_cc.SPQR_istat[4];
|
---|
215 | }
|
---|
216 | /// Set the fill-reducing ordering method to be used
|
---|
217 | void setSPQROrdering(int ord) { m_ordering = ord;}
|
---|
218 | /// Set the tolerance tol to treat columns with 2-norm < =tol as zero
|
---|
219 | void setPivotThreshold(const RealScalar& tol)
|
---|
220 | {
|
---|
221 | m_useDefaultThreshold = false;
|
---|
222 | m_tolerance = tol;
|
---|
223 | }
|
---|
224 |
|
---|
225 | /** \returns a pointer to the SPQR workspace */
|
---|
226 | cholmod_common *cholmodCommon() const { return &m_cc; }
|
---|
227 |
|
---|
228 |
|
---|
229 | /** \brief Reports whether previous computation was successful.
|
---|
230 | *
|
---|
231 | * \returns \c Success if computation was succesful,
|
---|
232 | * \c NumericalIssue if the sparse QR can not be computed
|
---|
233 | */
|
---|
234 | ComputationInfo info() const
|
---|
235 | {
|
---|
236 | eigen_assert(m_isInitialized && "Decomposition is not initialized.");
|
---|
237 | return m_info;
|
---|
238 | }
|
---|
239 | protected:
|
---|
240 | bool m_isInitialized;
|
---|
241 | bool m_analysisIsOk;
|
---|
242 | bool m_factorizationIsOk;
|
---|
243 | mutable bool m_isRUpToDate;
|
---|
244 | mutable ComputationInfo m_info;
|
---|
245 | int m_ordering; // Ordering method to use, see SPQR's manual
|
---|
246 | int m_allow_tol; // Allow to use some tolerance during numerical factorization.
|
---|
247 | RealScalar m_tolerance; // treat columns with 2-norm below this tolerance as zero
|
---|
248 | mutable cholmod_sparse *m_cR; // The sparse R factor in cholmod format
|
---|
249 | mutable MatrixType m_R; // The sparse matrix R in Eigen format
|
---|
250 | mutable Index *m_E; // The permutation applied to columns
|
---|
251 | mutable cholmod_sparse *m_H; //The householder vectors
|
---|
252 | mutable Index *m_HPinv; // The row permutation of H
|
---|
253 | mutable cholmod_dense *m_HTau; // The Householder coefficients
|
---|
254 | mutable Index m_rank; // The rank of the matrix
|
---|
255 | mutable cholmod_common m_cc; // Workspace and parameters
|
---|
256 | bool m_useDefaultThreshold; // Use default threshold
|
---|
257 | template<typename ,typename > friend struct SPQR_QProduct;
|
---|
258 | };
|
---|
259 |
|
---|
260 | template <typename SPQRType, typename Derived>
|
---|
261 | struct SPQR_QProduct : ReturnByValue<SPQR_QProduct<SPQRType,Derived> >
|
---|
262 | {
|
---|
263 | typedef typename SPQRType::Scalar Scalar;
|
---|
264 | typedef typename SPQRType::Index Index;
|
---|
265 | //Define the constructor to get reference to argument types
|
---|
266 | SPQR_QProduct(const SPQRType& spqr, const Derived& other, bool transpose) : m_spqr(spqr),m_other(other),m_transpose(transpose) {}
|
---|
267 |
|
---|
268 | inline Index rows() const { return m_transpose ? m_spqr.rows() : m_spqr.cols(); }
|
---|
269 | inline Index cols() const { return m_other.cols(); }
|
---|
270 | // Assign to a vector
|
---|
271 | template<typename ResType>
|
---|
272 | void evalTo(ResType& res) const
|
---|
273 | {
|
---|
274 | cholmod_dense y_cd;
|
---|
275 | cholmod_dense *x_cd;
|
---|
276 | int method = m_transpose ? SPQR_QTX : SPQR_QX;
|
---|
277 | cholmod_common *cc = m_spqr.cholmodCommon();
|
---|
278 | y_cd = viewAsCholmod(m_other.const_cast_derived());
|
---|
279 | x_cd = SuiteSparseQR_qmult<Scalar>(method, m_spqr.m_H, m_spqr.m_HTau, m_spqr.m_HPinv, &y_cd, cc);
|
---|
280 | res = Matrix<Scalar,ResType::RowsAtCompileTime,ResType::ColsAtCompileTime>::Map(reinterpret_cast<Scalar*>(x_cd->x), x_cd->nrow, x_cd->ncol);
|
---|
281 | cholmod_l_free_dense(&x_cd, cc);
|
---|
282 | }
|
---|
283 | const SPQRType& m_spqr;
|
---|
284 | const Derived& m_other;
|
---|
285 | bool m_transpose;
|
---|
286 |
|
---|
287 | };
|
---|
288 | template<typename SPQRType>
|
---|
289 | struct SPQRMatrixQReturnType{
|
---|
290 |
|
---|
291 | SPQRMatrixQReturnType(const SPQRType& spqr) : m_spqr(spqr) {}
|
---|
292 | template<typename Derived>
|
---|
293 | SPQR_QProduct<SPQRType, Derived> operator*(const MatrixBase<Derived>& other)
|
---|
294 | {
|
---|
295 | return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(),false);
|
---|
296 | }
|
---|
297 | SPQRMatrixQTransposeReturnType<SPQRType> adjoint() const
|
---|
298 | {
|
---|
299 | return SPQRMatrixQTransposeReturnType<SPQRType>(m_spqr);
|
---|
300 | }
|
---|
301 | // To use for operations with the transpose of Q
|
---|
302 | SPQRMatrixQTransposeReturnType<SPQRType> transpose() const
|
---|
303 | {
|
---|
304 | return SPQRMatrixQTransposeReturnType<SPQRType>(m_spqr);
|
---|
305 | }
|
---|
306 | const SPQRType& m_spqr;
|
---|
307 | };
|
---|
308 |
|
---|
309 | template<typename SPQRType>
|
---|
310 | struct SPQRMatrixQTransposeReturnType{
|
---|
311 | SPQRMatrixQTransposeReturnType(const SPQRType& spqr) : m_spqr(spqr) {}
|
---|
312 | template<typename Derived>
|
---|
313 | SPQR_QProduct<SPQRType,Derived> operator*(const MatrixBase<Derived>& other)
|
---|
314 | {
|
---|
315 | return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(), true);
|
---|
316 | }
|
---|
317 | const SPQRType& m_spqr;
|
---|
318 | };
|
---|
319 |
|
---|
320 | namespace internal {
|
---|
321 |
|
---|
322 | template<typename _MatrixType, typename Rhs>
|
---|
323 | struct solve_retval<SPQR<_MatrixType>, Rhs>
|
---|
324 | : solve_retval_base<SPQR<_MatrixType>, Rhs>
|
---|
325 | {
|
---|
326 | typedef SPQR<_MatrixType> Dec;
|
---|
327 | EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
|
---|
328 |
|
---|
329 | template<typename Dest> void evalTo(Dest& dst) const
|
---|
330 | {
|
---|
331 | dec()._solve(rhs(),dst);
|
---|
332 | }
|
---|
333 | };
|
---|
334 |
|
---|
335 | } // end namespace internal
|
---|
336 |
|
---|
337 | }// End namespace Eigen
|
---|
338 | #endif
|
---|