Eigen  3.2.91
SuiteSparseQRSupport.h
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 
57 template<typename _MatrixType>
58 class SPQR : public SparseSolverBase<SPQR<_MatrixType> >
59 {
60  protected:
61  typedef SparseSolverBase<SPQR<_MatrixType> > Base;
62  using Base::m_isInitialized;
63  public:
64  typedef typename _MatrixType::Scalar Scalar;
65  typedef typename _MatrixType::RealScalar RealScalar;
66  typedef UF_long StorageIndex ;
67  typedef SparseMatrix<Scalar, ColMajor, StorageIndex> MatrixType;
68  typedef Map<PermutationMatrix<Dynamic, Dynamic, StorageIndex> > PermutationType;
69  public:
70  SPQR()
71  : m_ordering(SPQR_ORDERING_DEFAULT), m_allow_tol(SPQR_DEFAULT_TOL), m_tolerance (NumTraits<Scalar>::epsilon()), m_useDefaultThreshold(true)
72  {
73  cholmod_l_start(&m_cc);
74  }
75 
76  explicit SPQR(const _MatrixType& matrix)
77  : m_ordering(SPQR_ORDERING_DEFAULT), m_allow_tol(SPQR_DEFAULT_TOL), m_tolerance (NumTraits<Scalar>::epsilon()), m_useDefaultThreshold(true)
78  {
79  cholmod_l_start(&m_cc);
80  compute(matrix);
81  }
82 
83  ~SPQR()
84  {
85  SPQR_free();
86  cholmod_l_finish(&m_cc);
87  }
88  void SPQR_free()
89  {
90  cholmod_l_free_sparse(&m_H, &m_cc);
91  cholmod_l_free_sparse(&m_cR, &m_cc);
92  cholmod_l_free_dense(&m_HTau, &m_cc);
93  std::free(m_E);
94  std::free(m_HPinv);
95  }
96 
97  void compute(const _MatrixType& matrix)
98  {
99  if(m_isInitialized) SPQR_free();
100 
101  MatrixType mat(matrix);
102 
103  /* Compute the default threshold as in MatLab, see:
104  * Tim Davis, "Algorithm 915, SuiteSparseQR: Multifrontal Multithreaded Rank-Revealing
105  * Sparse QR Factorization, ACM Trans. on Math. Soft. 38(1), 2011, Page 8:3
106  */
107  RealScalar pivotThreshold = m_tolerance;
108  if(m_useDefaultThreshold)
109  {
110  RealScalar max2Norm = 0.0;
111  for (int j = 0; j < mat.cols(); j++) max2Norm = numext::maxi(max2Norm, mat.col(j).norm());
112  if(max2Norm==RealScalar(0))
113  max2Norm = RealScalar(1);
114  pivotThreshold = 20 * (mat.rows() + mat.cols()) * max2Norm * NumTraits<RealScalar>::epsilon();
115  }
116 
117  cholmod_sparse A;
118  A = viewAsCholmod(mat);
119  Index col = matrix.cols();
120  m_rank = SuiteSparseQR<Scalar>(m_ordering, pivotThreshold, col, &A,
121  &m_cR, &m_E, &m_H, &m_HPinv, &m_HTau, &m_cc);
122 
123  if (!m_cR)
124  {
125  m_info = NumericalIssue;
126  m_isInitialized = false;
127  return;
128  }
129  m_info = Success;
130  m_isInitialized = true;
131  m_isRUpToDate = false;
132  }
136  inline Index rows() const {return m_cR->nrow; }
137 
141  inline Index cols() const { return m_cR->ncol; }
142 
143  template<typename Rhs, typename Dest>
144  void _solve_impl(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, y2;
151  y = matrixQ().transpose() * b;
152 
153  // Solves with the triangular matrix R
154  Index rk = this->rank();
155  y2 = y;
156  y.resize((std::max)(cols(),Index(y.rows())),y.cols());
157  y.topRows(rk) = this->matrixR().topLeftCorner(rk, rk).template triangularView<Upper>().solve(y2.topRows(rk));
158 
159  // Apply the column permutation
160  // colsPermutation() performs a copy of the permutation,
161  // so let's apply it manually:
162  for(Index i = 0; i < rk; ++i) dest.row(m_E[i]) = y.row(i);
163  for(Index i = rk; i < cols(); ++i) dest.row(m_E[i]).setZero();
164 
165 // y.bottomRows(y.rows()-rk).setZero();
166 // dest = colsPermutation() * y.topRows(cols());
167 
168  m_info = Success;
169  }
170 
173  const MatrixType matrixR() const
174  {
175  eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()");
176  if(!m_isRUpToDate) {
177  m_R = viewAsEigen<Scalar,ColMajor, typename MatrixType::StorageIndex>(*m_cR);
178  m_isRUpToDate = true;
179  }
180  return m_R;
181  }
183  SPQRMatrixQReturnType<SPQR> matrixQ() const
184  {
185  return SPQRMatrixQReturnType<SPQR>(*this);
186  }
188  PermutationType colsPermutation() const
189  {
190  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
191  return PermutationType(m_E, m_cR->ncol);
192  }
197  Index rank() const
198  {
199  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
200  return m_cc.SPQR_istat[4];
201  }
203  void setSPQROrdering(int ord) { m_ordering = ord;}
205  void setPivotThreshold(const RealScalar& tol)
206  {
207  m_useDefaultThreshold = false;
208  m_tolerance = tol;
209  }
210 
212  cholmod_common *cholmodCommon() const { return &m_cc; }
213 
214 
221  {
222  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
223  return m_info;
224  }
225  protected:
226  bool m_analysisIsOk;
227  bool m_factorizationIsOk;
228  mutable bool m_isRUpToDate;
229  mutable ComputationInfo m_info;
230  int m_ordering; // Ordering method to use, see SPQR's manual
231  int m_allow_tol; // Allow to use some tolerance during numerical factorization.
232  RealScalar m_tolerance; // treat columns with 2-norm below this tolerance as zero
233  mutable cholmod_sparse *m_cR; // The sparse R factor in cholmod format
234  mutable MatrixType m_R; // The sparse matrix R in Eigen format
235  mutable StorageIndex *m_E; // The permutation applied to columns
236  mutable cholmod_sparse *m_H; //The householder vectors
237  mutable StorageIndex *m_HPinv; // The row permutation of H
238  mutable cholmod_dense *m_HTau; // The Householder coefficients
239  mutable Index m_rank; // The rank of the matrix
240  mutable cholmod_common m_cc; // Workspace and parameters
241  bool m_useDefaultThreshold; // Use default threshold
242  template<typename ,typename > friend struct SPQR_QProduct;
243 };
244 
245 template <typename SPQRType, typename Derived>
246 struct SPQR_QProduct : ReturnByValue<SPQR_QProduct<SPQRType,Derived> >
247 {
248  typedef typename SPQRType::Scalar Scalar;
249  typedef typename SPQRType::StorageIndex StorageIndex;
250  //Define the constructor to get reference to argument types
251  SPQR_QProduct(const SPQRType& spqr, const Derived& other, bool transpose) : m_spqr(spqr),m_other(other),m_transpose(transpose) {}
252 
253  inline Index rows() const { return m_transpose ? m_spqr.rows() : m_spqr.cols(); }
254  inline Index cols() const { return m_other.cols(); }
255  // Assign to a vector
256  template<typename ResType>
257  void evalTo(ResType& res) const
258  {
259  cholmod_dense y_cd;
260  cholmod_dense *x_cd;
261  int method = m_transpose ? SPQR_QTX : SPQR_QX;
262  cholmod_common *cc = m_spqr.cholmodCommon();
263  y_cd = viewAsCholmod(m_other.const_cast_derived());
264  x_cd = SuiteSparseQR_qmult<Scalar>(method, m_spqr.m_H, m_spqr.m_HTau, m_spqr.m_HPinv, &y_cd, cc);
265  res = Matrix<Scalar,ResType::RowsAtCompileTime,ResType::ColsAtCompileTime>::Map(reinterpret_cast<Scalar*>(x_cd->x), x_cd->nrow, x_cd->ncol);
266  cholmod_l_free_dense(&x_cd, cc);
267  }
268  const SPQRType& m_spqr;
269  const Derived& m_other;
270  bool m_transpose;
271 
272 };
273 template<typename SPQRType>
274 struct SPQRMatrixQReturnType{
275 
276  SPQRMatrixQReturnType(const SPQRType& spqr) : m_spqr(spqr) {}
277  template<typename Derived>
278  SPQR_QProduct<SPQRType, Derived> operator*(const MatrixBase<Derived>& other)
279  {
280  return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(),false);
281  }
282  SPQRMatrixQTransposeReturnType<SPQRType> adjoint() const
283  {
284  return SPQRMatrixQTransposeReturnType<SPQRType>(m_spqr);
285  }
286  // To use for operations with the transpose of Q
287  SPQRMatrixQTransposeReturnType<SPQRType> transpose() const
288  {
289  return SPQRMatrixQTransposeReturnType<SPQRType>(m_spqr);
290  }
291  const SPQRType& m_spqr;
292 };
293 
294 template<typename SPQRType>
295 struct SPQRMatrixQTransposeReturnType{
296  SPQRMatrixQTransposeReturnType(const SPQRType& spqr) : m_spqr(spqr) {}
297  template<typename Derived>
298  SPQR_QProduct<SPQRType,Derived> operator*(const MatrixBase<Derived>& other)
299  {
300  return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(), true);
301  }
302  const SPQRType& m_spqr;
303 };
304 
305 }// End namespace Eigen
306 #endif
cholmod_common * cholmodCommon() const
Definition: SuiteSparseQRSupport.h:212
Index rows() const
Definition: SuiteSparseQRSupport.h:136
const Solve< SPQR< _MatrixType >, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: SparseSolverBase.h:74
Definition: LDLT.h:16
RowXpr row(Index i)
Definition: DenseBase.h:797
Index cols() const
Definition: SuiteSparseQRSupport.h:141
void setPivotThreshold(const RealScalar &tol)
Set the tolerance tol to treat columns with 2-norm < =tol as zero.
Definition: SuiteSparseQRSupport.h:205
PermutationType colsPermutation() const
Get the permutation that was applied to columns of A.
Definition: SuiteSparseQRSupport.h:188
Definition: Constants.h:426
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SuiteSparseQRSupport.h:220
Index rank() const
Definition: SuiteSparseQRSupport.h:197
Definition: Constants.h:424
Block< Derived > topLeftCorner(Index cRows, Index cCols)
Definition: SparseMatrixBase.h:164
Definition: Eigen_Colamd.h:54
const MatrixType matrixR() const
Definition: SuiteSparseQRSupport.h:173
void setSPQROrdering(int ord)
Set the fill-reducing ordering method to be used.
Definition: SuiteSparseQRSupport.h:203
SPQRMatrixQReturnType< SPQR > matrixQ() const
Get an expression of the matrix Q.
Definition: SuiteSparseQRSupport.h:183
Sparse QR factorization based on SuiteSparseQR library.
Definition: SuiteSparseQRSupport.h:16
ComputationInfo
Definition: Constants.h:422
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48