Eigen  3.4.90 (git rev 67eeba6e720c5745abc77ae6c92ce0a44aa7b7ae)
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 #include "./InternalHeaderCheck.h"
15 
16 namespace Eigen {
17 
18  template<typename MatrixType> class SPQR;
19  template<typename SPQRType> struct SPQRMatrixQReturnType;
20  template<typename SPQRType> struct SPQRMatrixQTransposeReturnType;
21  template <typename SPQRType, typename Derived> struct SPQR_QProduct;
22  namespace internal {
23  template <typename SPQRType> struct traits<SPQRMatrixQReturnType<SPQRType> >
24  {
25  typedef typename SPQRType::MatrixType ReturnType;
26  };
27  template <typename SPQRType> struct traits<SPQRMatrixQTransposeReturnType<SPQRType> >
28  {
29  typedef typename SPQRType::MatrixType ReturnType;
30  };
31  template <typename SPQRType, typename Derived> struct traits<SPQR_QProduct<SPQRType, Derived> >
32  {
33  typedef typename Derived::PlainObject ReturnType;
34  };
35  } // End namespace internal
36 
61 template<typename MatrixType_>
62 class SPQR : public SparseSolverBase<SPQR<MatrixType_> >
63 {
64  protected:
66  using Base::m_isInitialized;
67  public:
68  typedef typename MatrixType_::Scalar Scalar;
69  typedef typename MatrixType_::RealScalar RealScalar;
70  typedef SuiteSparse_long StorageIndex ;
73  enum {
74  ColsAtCompileTime = Dynamic,
75  MaxColsAtCompileTime = Dynamic
76  };
77  public:
78  SPQR()
79  : m_analysisIsOk(false),
80  m_factorizationIsOk(false),
81  m_isRUpToDate(false),
82  m_ordering(SPQR_ORDERING_DEFAULT),
83  m_allow_tol(SPQR_DEFAULT_TOL),
84  m_tolerance (NumTraits<Scalar>::epsilon()),
85  m_cR(0),
86  m_E(0),
87  m_H(0),
88  m_HPinv(0),
89  m_HTau(0),
90  m_useDefaultThreshold(true)
91  {
92  cholmod_l_start(&m_cc);
93  }
94 
95  explicit SPQR(const MatrixType_& matrix)
96  : m_analysisIsOk(false),
97  m_factorizationIsOk(false),
98  m_isRUpToDate(false),
99  m_ordering(SPQR_ORDERING_DEFAULT),
100  m_allow_tol(SPQR_DEFAULT_TOL),
101  m_tolerance (NumTraits<Scalar>::epsilon()),
102  m_cR(0),
103  m_E(0),
104  m_H(0),
105  m_HPinv(0),
106  m_HTau(0),
107  m_useDefaultThreshold(true)
108  {
109  cholmod_l_start(&m_cc);
110  compute(matrix);
111  }
112 
113  ~SPQR()
114  {
115  SPQR_free();
116  cholmod_l_finish(&m_cc);
117  }
118  void SPQR_free()
119  {
120  cholmod_l_free_sparse(&m_H, &m_cc);
121  cholmod_l_free_sparse(&m_cR, &m_cc);
122  cholmod_l_free_dense(&m_HTau, &m_cc);
123  std::free(m_E);
124  std::free(m_HPinv);
125  }
126 
127  void compute(const MatrixType_& matrix)
128  {
129  if(m_isInitialized) SPQR_free();
130 
131  MatrixType mat(matrix);
132 
133  /* Compute the default threshold as in MatLab, see:
134  * Tim Davis, "Algorithm 915, SuiteSparseQR: Multifrontal Multithreaded Rank-Revealing
135  * Sparse QR Factorization, ACM Trans. on Math. Soft. 38(1), 2011, Page 8:3
136  */
137  RealScalar pivotThreshold = m_tolerance;
138  if(m_useDefaultThreshold)
139  {
140  RealScalar max2Norm = 0.0;
141  for (int j = 0; j < mat.cols(); j++) max2Norm = numext::maxi(max2Norm, mat.col(j).norm());
142  if(numext::is_exactly_zero(max2Norm))
143  max2Norm = RealScalar(1);
144  pivotThreshold = 20 * (mat.rows() + mat.cols()) * max2Norm * NumTraits<RealScalar>::epsilon();
145  }
146  cholmod_sparse A;
147  A = viewAsCholmod(mat);
148  m_rows = matrix.rows();
149  Index col = matrix.cols();
150  m_rank = SuiteSparseQR<Scalar>(m_ordering, pivotThreshold, col, &A,
151  &m_cR, &m_E, &m_H, &m_HPinv, &m_HTau, &m_cc);
152 
153  if (!m_cR)
154  {
155  m_info = NumericalIssue;
156  m_isInitialized = false;
157  return;
158  }
159  m_info = Success;
160  m_isInitialized = true;
161  m_isRUpToDate = false;
162  }
166  inline Index rows() const {return m_rows; }
167 
171  inline Index cols() const { return m_cR->ncol; }
172 
173  template<typename Rhs, typename Dest>
174  void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
175  {
176  eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()");
177  eigen_assert(b.cols()==1 && "This method is for vectors only");
178 
179  //Compute Q^T * b
180  typename Dest::PlainObject y, y2;
181  y = matrixQ().transpose() * b;
182 
183  // Solves with the triangular matrix R
184  Index rk = this->rank();
185  y2 = y;
186  y.resize((std::max)(cols(),Index(y.rows())),y.cols());
187  y.topRows(rk) = this->matrixR().topLeftCorner(rk, rk).template triangularView<Upper>().solve(y2.topRows(rk));
188 
189  // Apply the column permutation
190  // colsPermutation() performs a copy of the permutation,
191  // so let's apply it manually:
192  for(Index i = 0; i < rk; ++i) dest.row(m_E[i]) = y.row(i);
193  for(Index i = rk; i < cols(); ++i) dest.row(m_E[i]).setZero();
194 
195 // y.bottomRows(y.rows()-rk).setZero();
196 // dest = colsPermutation() * y.topRows(cols());
197 
198  m_info = Success;
199  }
200 
203  const MatrixType matrixR() const
204  {
205  eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()");
206  if(!m_isRUpToDate) {
207  m_R = viewAsEigen<Scalar,ColMajor, typename MatrixType::StorageIndex>(*m_cR);
208  m_isRUpToDate = true;
209  }
210  return m_R;
211  }
213  SPQRMatrixQReturnType<SPQR> matrixQ() const
214  {
215  return SPQRMatrixQReturnType<SPQR>(*this);
216  }
219  {
220  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
221  return PermutationType(m_E, m_cR->ncol);
222  }
227  Index rank() const
228  {
229  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
230  return m_cc.SPQR_istat[4];
231  }
233  void setSPQROrdering(int ord) { m_ordering = ord;}
235  void setPivotThreshold(const RealScalar& tol)
236  {
237  m_useDefaultThreshold = false;
238  m_tolerance = tol;
239  }
240 
242  cholmod_common *cholmodCommon() const { return &m_cc; }
243 
244 
251  {
252  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
253  return m_info;
254  }
255  protected:
256  bool m_analysisIsOk;
257  bool m_factorizationIsOk;
258  mutable bool m_isRUpToDate;
259  mutable ComputationInfo m_info;
260  int m_ordering; // Ordering method to use, see SPQR's manual
261  int m_allow_tol; // Allow to use some tolerance during numerical factorization.
262  RealScalar m_tolerance; // treat columns with 2-norm below this tolerance as zero
263  mutable cholmod_sparse *m_cR = nullptr; // The sparse R factor in cholmod format
264  mutable MatrixType m_R; // The sparse matrix R in Eigen format
265  mutable StorageIndex *m_E = nullptr; // The permutation applied to columns
266  mutable cholmod_sparse *m_H = nullptr; //The householder vectors
267  mutable StorageIndex *m_HPinv = nullptr; // The row permutation of H
268  mutable cholmod_dense *m_HTau = nullptr; // The Householder coefficients
269  mutable Index m_rank; // The rank of the matrix
270  mutable cholmod_common m_cc; // Workspace and parameters
271  bool m_useDefaultThreshold; // Use default threshold
272  Index m_rows;
273  template<typename ,typename > friend struct SPQR_QProduct;
274 };
275 
276 template <typename SPQRType, typename Derived>
277 struct SPQR_QProduct : ReturnByValue<SPQR_QProduct<SPQRType,Derived> >
278 {
279  typedef typename SPQRType::Scalar Scalar;
280  typedef typename SPQRType::StorageIndex StorageIndex;
281  //Define the constructor to get reference to argument types
282  SPQR_QProduct(const SPQRType& spqr, const Derived& other, bool transpose) : m_spqr(spqr),m_other(other),m_transpose(transpose) {}
283 
284  inline Index rows() const { return m_transpose ? m_spqr.rows() : m_spqr.cols(); }
285  inline Index cols() const { return m_other.cols(); }
286  // Assign to a vector
287  template<typename ResType>
288  void evalTo(ResType& res) const
289  {
290  cholmod_dense y_cd;
291  cholmod_dense *x_cd;
292  int method = m_transpose ? SPQR_QTX : SPQR_QX;
293  cholmod_common *cc = m_spqr.cholmodCommon();
294  y_cd = viewAsCholmod(m_other.const_cast_derived());
295  x_cd = SuiteSparseQR_qmult<Scalar>(method, m_spqr.m_H, m_spqr.m_HTau, m_spqr.m_HPinv, &y_cd, cc);
296  res = Matrix<Scalar,ResType::RowsAtCompileTime,ResType::ColsAtCompileTime>::Map(reinterpret_cast<Scalar*>(x_cd->x), x_cd->nrow, x_cd->ncol);
297  cholmod_l_free_dense(&x_cd, cc);
298  }
299  const SPQRType& m_spqr;
300  const Derived& m_other;
301  bool m_transpose;
302 
303 };
304 template<typename SPQRType>
305 struct SPQRMatrixQReturnType{
306 
307  SPQRMatrixQReturnType(const SPQRType& spqr) : m_spqr(spqr) {}
308  template<typename Derived>
309  SPQR_QProduct<SPQRType, Derived> operator*(const MatrixBase<Derived>& other)
310  {
311  return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(),false);
312  }
313  SPQRMatrixQTransposeReturnType<SPQRType> adjoint() const
314  {
315  return SPQRMatrixQTransposeReturnType<SPQRType>(m_spqr);
316  }
317  // To use for operations with the transpose of Q
318  SPQRMatrixQTransposeReturnType<SPQRType> transpose() const
319  {
320  return SPQRMatrixQTransposeReturnType<SPQRType>(m_spqr);
321  }
322  const SPQRType& m_spqr;
323 };
324 
325 template<typename SPQRType>
326 struct SPQRMatrixQTransposeReturnType{
327  SPQRMatrixQTransposeReturnType(const SPQRType& spqr) : m_spqr(spqr) {}
328  template<typename Derived>
329  SPQR_QProduct<SPQRType,Derived> operator*(const MatrixBase<Derived>& other)
330  {
331  return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(), true);
332  }
333  const SPQRType& m_spqr;
334 };
335 
336 }// End namespace Eigen
337 #endif
Derived & setZero()
Definition: CwiseNullaryOp.h:548
EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Definition: EigenBase.h:65
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:98
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:52
static ConstMapType Map(const Scalar *data)
Definition: PlainObjectBase.h:642
Sparse QR factorization based on SuiteSparseQR library.
Definition: SuiteSparseQRSupport.h:63
void setPivotThreshold(const RealScalar &tol)
Set the tolerance tol to treat columns with 2-norm < =tol as zero.
Definition: SuiteSparseQRSupport.h:235
Index rank() const
Definition: SuiteSparseQRSupport.h:227
cholmod_common * cholmodCommon() const
Definition: SuiteSparseQRSupport.h:242
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SuiteSparseQRSupport.h:250
Index rows() const
Definition: SuiteSparseQRSupport.h:166
const MatrixType matrixR() const
Definition: SuiteSparseQRSupport.h:203
Index cols() const
Definition: SuiteSparseQRSupport.h:171
SPQRMatrixQReturnType< SPQR > matrixQ() const
Get an expression of the matrix Q.
Definition: SuiteSparseQRSupport.h:213
PermutationType colsPermutation() const
Get the permutation that was applied to columns of A.
Definition: SuiteSparseQRSupport.h:218
void setSPQROrdering(int ord)
Set the fill-reducing ordering method to be used.
Definition: SuiteSparseQRSupport.h:233
Index cols() const
Definition: SparseMatrix.h:142
Index rows() const
Definition: SparseMatrix.h:140
A base class for sparse solvers.
Definition: SparseSolverBase.h:70
ComputationInfo
Definition: Constants.h:442
@ NumericalIssue
Definition: Constants.h:446
@ Success
Definition: Constants.h:444
Namespace containing all symbols from the Eigen library.
Definition: Core:139
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
cholmod_sparse viewAsCholmod(Ref< SparseMatrix< Scalar_, Options_, StorageIndex_ > > mat)
Definition: CholmodSupport.h:60
const Product< MatrixDerived, PermutationDerived, AliasFreeProduct > operator*(const MatrixBase< MatrixDerived > &matrix, const PermutationBase< PermutationDerived > &permutation)
Definition: PermutationMatrix.h:517
const int Dynamic
Definition: Constants.h:24
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:231