SuiteSparseQRSupport.h
Go to the documentation of this file.
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 
59 template<typename _MatrixType>
60 class SPQR : public SparseSolverBase<SPQR<_MatrixType> >
61 {
62  protected:
64  using Base::m_isInitialized;
65  public:
66  typedef typename _MatrixType::Scalar Scalar;
71  enum {
72  ColsAtCompileTime = Dynamic,
73  MaxColsAtCompileTime = Dynamic
74  };
75  public:
76  SPQR()
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  }
81 
82  explicit SPQR(const _MatrixType& matrix)
83  : m_ordering(SPQR_ORDERING_DEFAULT), m_allow_tol(SPQR_DEFAULT_TOL), m_tolerance (NumTraits<Scalar>::epsilon()), m_useDefaultThreshold(true)
84  {
85  cholmod_l_start(&m_cc);
86  compute(matrix);
87  }
88 
90  {
91  SPQR_free();
92  cholmod_l_finish(&m_cc);
93  }
94  void SPQR_free()
95  {
96  cholmod_l_free_sparse(&m_H, &m_cc);
97  cholmod_l_free_sparse(&m_cR, &m_cc);
98  cholmod_l_free_dense(&m_HTau, &m_cc);
99  std::free(m_E);
100  std::free(m_HPinv);
101  }
102 
103  void compute(const _MatrixType& matrix)
104  {
105  if(m_isInitialized) SPQR_free();
106 
107  MatrixType mat(matrix);
108 
109  /* Compute the default threshold as in MatLab, see:
110  * Tim Davis, "Algorithm 915, SuiteSparseQR: Multifrontal Multithreaded Rank-Revealing
111  * Sparse QR Factorization, ACM Trans. on Math. Soft. 38(1), 2011, Page 8:3
112  */
113  RealScalar pivotThreshold = m_tolerance;
114  if(m_useDefaultThreshold)
115  {
116  RealScalar max2Norm = 0.0;
117  for (int j = 0; j < mat.cols(); j++) max2Norm = numext::maxi(max2Norm, mat.col(j).norm());
118  if(max2Norm==RealScalar(0))
119  max2Norm = RealScalar(1);
120  pivotThreshold = 20 * (mat.rows() + mat.cols()) * max2Norm * NumTraits<RealScalar>::epsilon();
121  }
122  cholmod_sparse A;
123  A = viewAsCholmod(mat);
124  m_rows = matrix.rows();
125  Index col = matrix.cols();
126  m_rank = SuiteSparseQR<Scalar>(m_ordering, pivotThreshold, col, &A,
127  &m_cR, &m_E, &m_H, &m_HPinv, &m_HTau, &m_cc);
128 
129  if (!m_cR)
130  {
131  m_info = NumericalIssue;
132  m_isInitialized = false;
133  return;
134  }
135  m_info = Success;
136  m_isInitialized = true;
137  m_isRUpToDate = false;
138  }
142  inline Index rows() const {return m_rows; }
143 
147  inline Index cols() const { return m_cR->ncol; }
148 
149  template<typename Rhs, typename Dest>
150  void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
151  {
152  eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()");
153  eigen_assert(b.cols()==1 && "This method is for vectors only");
154 
155  //Compute Q^T * b
156  typename Dest::PlainObject y, y2;
157  y = matrixQ().transpose() * b;
158 
159  // Solves with the triangular matrix R
160  Index rk = this->rank();
161  y2 = y;
162  y.resize((std::max)(cols(),Index(y.rows())),y.cols());
163  y.topRows(rk) = this->matrixR().topLeftCorner(rk, rk).template triangularView<Upper>().solve(y2.topRows(rk));
164 
165  // Apply the column permutation
166  // colsPermutation() performs a copy of the permutation,
167  // so let's apply it manually:
168  for(Index i = 0; i < rk; ++i) dest.row(m_E[i]) = y.row(i);
169  for(Index i = rk; i < cols(); ++i) dest.row(m_E[i]).setZero();
170 
171 // y.bottomRows(y.rows()-rk).setZero();
172 // dest = colsPermutation() * y.topRows(cols());
173 
174  m_info = Success;
175  }
176 
179  const MatrixType matrixR() const
180  {
181  eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()");
182  if(!m_isRUpToDate) {
183  m_R = viewAsEigen<Scalar,ColMajor, typename MatrixType::StorageIndex>(*m_cR);
184  m_isRUpToDate = true;
185  }
186  return m_R;
187  }
190  {
191  return SPQRMatrixQReturnType<SPQR>(*this);
192  }
194  PermutationType colsPermutation() const
195  {
196  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
197  return PermutationType(m_E, m_cR->ncol);
198  }
203  Index rank() const
204  {
205  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
206  return m_cc.SPQR_istat[4];
207  }
209  void setSPQROrdering(int ord) { m_ordering = ord;}
211  void setPivotThreshold(const RealScalar& tol)
212  {
213  m_useDefaultThreshold = false;
214  m_tolerance = tol;
215  }
216 
218  cholmod_common *cholmodCommon() const { return &m_cc; }
219 
220 
227  {
228  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
229  return m_info;
230  }
231  protected:
234  mutable bool m_isRUpToDate;
236  int m_ordering; // Ordering method to use, see SPQR's manual
237  int m_allow_tol; // Allow to use some tolerance during numerical factorization.
238  RealScalar m_tolerance; // treat columns with 2-norm below this tolerance as zero
239  mutable cholmod_sparse *m_cR; // The sparse R factor in cholmod format
240  mutable MatrixType m_R; // The sparse matrix R in Eigen format
241  mutable StorageIndex *m_E; // The permutation applied to columns
242  mutable cholmod_sparse *m_H; //The householder vectors
243  mutable StorageIndex *m_HPinv; // The row permutation of H
244  mutable cholmod_dense *m_HTau; // The Householder coefficients
245  mutable Index m_rank; // The rank of the matrix
246  mutable cholmod_common m_cc; // Workspace and parameters
247  bool m_useDefaultThreshold; // Use default threshold
249  template<typename ,typename > friend struct SPQR_QProduct;
250 };
251 
252 template <typename SPQRType, typename Derived>
253 struct SPQR_QProduct : ReturnByValue<SPQR_QProduct<SPQRType,Derived> >
254 {
255  typedef typename SPQRType::Scalar Scalar;
256  typedef typename SPQRType::StorageIndex StorageIndex;
257  //Define the constructor to get reference to argument types
258  SPQR_QProduct(const SPQRType& spqr, const Derived& other, bool transpose) : m_spqr(spqr),m_other(other),m_transpose(transpose) {}
259 
260  inline Index rows() const { return m_transpose ? m_spqr.rows() : m_spqr.cols(); }
261  inline Index cols() const { return m_other.cols(); }
262  // Assign to a vector
263  template<typename ResType>
264  void evalTo(ResType& res) const
265  {
266  cholmod_dense y_cd;
267  cholmod_dense *x_cd;
268  int method = m_transpose ? SPQR_QTX : SPQR_QX;
269  cholmod_common *cc = m_spqr.cholmodCommon();
270  y_cd = viewAsCholmod(m_other.const_cast_derived());
271  x_cd = SuiteSparseQR_qmult<Scalar>(method, m_spqr.m_H, m_spqr.m_HTau, m_spqr.m_HPinv, &y_cd, cc);
272  res = Matrix<Scalar,ResType::RowsAtCompileTime,ResType::ColsAtCompileTime>::Map(reinterpret_cast<Scalar*>(x_cd->x), x_cd->nrow, x_cd->ncol);
273  cholmod_l_free_dense(&x_cd, cc);
274  }
275  const SPQRType& m_spqr;
276  const Derived& m_other;
277  bool m_transpose;
278 
279 };
280 template<typename SPQRType>
281 struct SPQRMatrixQReturnType{
282 
283  SPQRMatrixQReturnType(const SPQRType& spqr) : m_spqr(spqr) {}
284  template<typename Derived>
286  {
287  return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(),false);
288  }
290  {
292  }
293  // To use for operations with the transpose of Q
295  {
297  }
298  const SPQRType& m_spqr;
299 };
300 
301 template<typename SPQRType>
303  SPQRMatrixQTransposeReturnType(const SPQRType& spqr) : m_spqr(spqr) {}
304  template<typename Derived>
306  {
307  return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(), true);
308  }
309  const SPQRType& m_spqr;
310 };
311 
312 }// End namespace Eigen
313 #endif
SPQR_QProduct< SPQRType, Derived > operator*(const MatrixBase< Derived > &other)
cholmod_common * cholmodCommon() const
SCALAR Scalar
Definition: bench_gemm.cpp:33
SuiteSparse_long StorageIndex
#define max(a, b)
Definition: datatypes.h:20
SPQRMatrixQTransposeReturnType< SPQRType > transpose() const
Scalar * b
Definition: benchVecAdd.cpp:17
Map< PermutationMatrix< Dynamic, Dynamic, StorageIndex > > PermutationType
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
SPQR(const _MatrixType &matrix)
A base class for sparse solvers.
SPQRMatrixQReturnType(const SPQRType &spqr)
Index rows() const
StorageIndex * m_HPinv
ComputationInfo m_info
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
MatrixXf MatrixType
SPQR_QProduct< SPQRType, Derived > operator*(const MatrixBase< Derived > &other)
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:150
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T maxi(const T &x, const T &y)
_MatrixType::Scalar Scalar
SparseMatrix< Scalar, ColMajor, StorageIndex > MatrixType
static double epsilon
Definition: testRot3.cpp:39
Index rows() const
Definition: SparseMatrix.h:136
RealScalar m_tolerance
Index cols() const
cholmod_common m_cc
cout<< "Here is the matrix m:"<< endl<< m<< endl;Matrix< ptrdiff_t, 3, 1 > res
void setPivotThreshold(const RealScalar &tol)
Set the tolerance tol to treat columns with 2-norm < =tol as zero.
cholmod_dense * m_HTau
void _solve_impl(const MatrixBase< Rhs > &b, MatrixBase< Dest > &dest) const
#define SuiteSparse_long
PermutationType colsPermutation() const
Get the permutation that was applied to columns of A.
Index cols() const
Definition: SparseMatrix.h:138
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
ComputationInfo info() const
Reports whether previous computation was successful.
SPQR_QProduct(const SPQRType &spqr, const Derived &other, bool transpose)
#define eigen_assert(x)
Definition: Macros.h:579
Index rank() const
cholmod_sparse * m_cR
void compute(const _MatrixType &matrix)
EIGEN_DEVICE_FUNC ColXpr col(Index i)
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:34
cholmod_sparse viewAsCholmod(Ref< SparseMatrix< _Scalar, _Options, _StorageIndex > > mat)
SparseSolverBase< SPQR< _MatrixType > > Base
const MatrixType matrixR() const
cholmod_sparse * m_H
void setSPQROrdering(int ord)
Set the fill-reducing ordering method to be used.
SPQRMatrixQReturnType< SPQR > matrixQ() const
Get an expression of the matrix Q.
SPQRMatrixQTransposeReturnType(const SPQRType &spqr)
Sparse QR factorization based on SuiteSparseQR library.
SPQRType::StorageIndex StorageIndex
m col(1)
void evalTo(ResType &res) const
SPQRMatrixQTransposeReturnType< SPQRType > adjoint() const
const G double tol
Definition: Group.h:83
const int Dynamic
Definition: Constants.h:21
StorageIndex * m_E
EIGEN_DEVICE_FUNC RowXpr row(Index i)
Definition: DenseBase.h:860
EIGEN_DONT_INLINE void compute(Solver &solver, const MatrixType &A)
Map< Matrix< T, Dynamic, Dynamic, ColMajor >, 0, OuterStride<> > matrix(T *data, int rows, int cols, int stride)
ComputationInfo
Definition: Constants.h:430
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
std::ptrdiff_t j
_MatrixType::RealScalar RealScalar


gtsam
Author(s):
autogenerated on Sat May 8 2021 02:44:51