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:
65  public:
66  typedef typename _MatrixType::Scalar Scalar;
71  enum {
74  };
75  public:
76  SPQR()
77  : m_analysisIsOk(false),
78  m_factorizationIsOk(false),
79  m_isRUpToDate(false),
80  m_ordering(SPQR_ORDERING_DEFAULT),
81  m_allow_tol(SPQR_DEFAULT_TOL),
83  m_cR(0),
84  m_E(0),
85  m_H(0),
86  m_HPinv(0),
87  m_HTau(0),
89  {
90  cholmod_l_start(&m_cc);
91  }
92 
93  explicit SPQR(const _MatrixType& matrix)
94  : m_analysisIsOk(false),
95  m_factorizationIsOk(false),
96  m_isRUpToDate(false),
97  m_ordering(SPQR_ORDERING_DEFAULT),
98  m_allow_tol(SPQR_DEFAULT_TOL),
100  m_cR(0),
101  m_E(0),
102  m_H(0),
103  m_HPinv(0),
104  m_HTau(0),
106  {
107  cholmod_l_start(&m_cc);
108  compute(matrix);
109  }
110 
112  {
113  SPQR_free();
114  cholmod_l_finish(&m_cc);
115  }
116  void SPQR_free()
117  {
118  cholmod_l_free_sparse(&m_H, &m_cc);
119  cholmod_l_free_sparse(&m_cR, &m_cc);
120  cholmod_l_free_dense(&m_HTau, &m_cc);
121  std::free(m_E);
122  std::free(m_HPinv);
123  }
124 
125  void compute(const _MatrixType& matrix)
126  {
128 
130 
131  /* Compute the default threshold as in MatLab, see:
132  * Tim Davis, "Algorithm 915, SuiteSparseQR: Multifrontal Multithreaded Rank-Revealing
133  * Sparse QR Factorization, ACM Trans. on Math. Soft. 38(1), 2011, Page 8:3
134  */
135  RealScalar pivotThreshold = m_tolerance;
137  {
138  RealScalar max2Norm = 0.0;
139  for (int j = 0; j < mat.cols(); j++) max2Norm = numext::maxi(max2Norm, mat.col(j).norm());
140  if(max2Norm==RealScalar(0))
141  max2Norm = RealScalar(1);
142  pivotThreshold = 20 * (mat.rows() + mat.cols()) * max2Norm * NumTraits<RealScalar>::epsilon();
143  }
144  cholmod_sparse A;
145  A = viewAsCholmod(mat);
146  m_rows = matrix.rows();
147  Index col = matrix.cols();
148  m_rank = SuiteSparseQR<Scalar>(m_ordering, pivotThreshold, col, &A,
149  &m_cR, &m_E, &m_H, &m_HPinv, &m_HTau, &m_cc);
150 
151  if (!m_cR)
152  {
154  m_isInitialized = false;
155  return;
156  }
157  m_info = Success;
158  m_isInitialized = true;
159  m_isRUpToDate = false;
160  }
164  inline Index rows() const {return m_rows; }
165 
169  inline Index cols() const { return m_cR->ncol; }
170 
171  template<typename Rhs, typename Dest>
172  void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
173  {
174  eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()");
175  eigen_assert(b.cols()==1 && "This method is for vectors only");
176 
177  //Compute Q^T * b
178  typename Dest::PlainObject y, y2;
179  y = matrixQ().transpose() * b;
180 
181  // Solves with the triangular matrix R
182  Index rk = this->rank();
183  y2 = y;
184  y.resize((std::max)(cols(),Index(y.rows())),y.cols());
185  y.topRows(rk) = this->matrixR().topLeftCorner(rk, rk).template triangularView<Upper>().solve(y2.topRows(rk));
186 
187  // Apply the column permutation
188  // colsPermutation() performs a copy of the permutation,
189  // so let's apply it manually:
190  for(Index i = 0; i < rk; ++i) dest.row(m_E[i]) = y.row(i);
191  for(Index i = rk; i < cols(); ++i) dest.row(m_E[i]).setZero();
192 
193 // y.bottomRows(y.rows()-rk).setZero();
194 // dest = colsPermutation() * y.topRows(cols());
195 
196  m_info = Success;
197  }
198 
201  const MatrixType matrixR() const
202  {
203  eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()");
204  if(!m_isRUpToDate) {
205  m_R = viewAsEigen<Scalar,ColMajor, typename MatrixType::StorageIndex>(*m_cR);
206  m_isRUpToDate = true;
207  }
208  return m_R;
209  }
212  {
213  return SPQRMatrixQReturnType<SPQR>(*this);
214  }
217  {
218  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
219  return PermutationType(m_E, m_cR->ncol);
220  }
225  Index rank() const
226  {
227  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
228  return m_cc.SPQR_istat[4];
229  }
231  void setSPQROrdering(int ord) { m_ordering = ord;}
234  {
235  m_useDefaultThreshold = false;
236  m_tolerance = tol;
237  }
238 
240  cholmod_common *cholmodCommon() const { return &m_cc; }
241 
242 
249  {
250  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
251  return m_info;
252  }
253  protected:
256  mutable bool m_isRUpToDate;
258  int m_ordering; // Ordering method to use, see SPQR's manual
259  int m_allow_tol; // Allow to use some tolerance during numerical factorization.
260  RealScalar m_tolerance; // treat columns with 2-norm below this tolerance as zero
261  mutable cholmod_sparse *m_cR; // The sparse R factor in cholmod format
262  mutable MatrixType m_R; // The sparse matrix R in Eigen format
263  mutable StorageIndex *m_E; // The permutation applied to columns
264  mutable cholmod_sparse *m_H; //The householder vectors
265  mutable StorageIndex *m_HPinv; // The row permutation of H
266  mutable cholmod_dense *m_HTau; // The Householder coefficients
267  mutable Index m_rank; // The rank of the matrix
268  mutable cholmod_common m_cc; // Workspace and parameters
269  bool m_useDefaultThreshold; // Use default threshold
271  template<typename ,typename > friend struct SPQR_QProduct;
272 };
273 
274 template <typename SPQRType, typename Derived>
275 struct SPQR_QProduct : ReturnByValue<SPQR_QProduct<SPQRType,Derived> >
276 {
277  typedef typename SPQRType::Scalar Scalar;
278  typedef typename SPQRType::StorageIndex StorageIndex;
279  //Define the constructor to get reference to argument types
280  SPQR_QProduct(const SPQRType& spqr, const Derived& other, bool transpose) : m_spqr(spqr),m_other(other),m_transpose(transpose) {}
281 
282  inline Index rows() const { return m_transpose ? m_spqr.rows() : m_spqr.cols(); }
283  inline Index cols() const { return m_other.cols(); }
284  // Assign to a vector
285  template<typename ResType>
286  void evalTo(ResType& res) const
287  {
288  cholmod_dense y_cd;
289  cholmod_dense *x_cd;
290  int method = m_transpose ? SPQR_QTX : SPQR_QX;
291  cholmod_common *cc = m_spqr.cholmodCommon();
292  y_cd = viewAsCholmod(m_other.const_cast_derived());
293  x_cd = SuiteSparseQR_qmult<Scalar>(method, m_spqr.m_H, m_spqr.m_HTau, m_spqr.m_HPinv, &y_cd, cc);
294  res = Matrix<Scalar,ResType::RowsAtCompileTime,ResType::ColsAtCompileTime>::Map(reinterpret_cast<Scalar*>(x_cd->x), x_cd->nrow, x_cd->ncol);
295  cholmod_l_free_dense(&x_cd, cc);
296  }
297  const SPQRType& m_spqr;
298  const Derived& m_other;
299  bool m_transpose;
300 
301 };
302 template<typename SPQRType>
303 struct SPQRMatrixQReturnType{
304 
305  SPQRMatrixQReturnType(const SPQRType& spqr) : m_spqr(spqr) {}
306  template<typename Derived>
308  {
309  return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(),false);
310  }
312  {
314  }
315  // To use for operations with the transpose of Q
317  {
319  }
320  const SPQRType& m_spqr;
321 };
322 
323 template<typename SPQRType>
325  SPQRMatrixQTransposeReturnType(const SPQRType& spqr) : m_spqr(spqr) {}
326  template<typename Derived>
328  {
329  return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(), true);
330  }
331  const SPQRType& m_spqr;
332 };
333 
334 }// End namespace Eigen
335 #endif
Eigen::SPQR::~SPQR
~SPQR()
Definition: SuiteSparseQRSupport.h:111
Eigen::SPQR::cholmodCommon
cholmod_common * cholmodCommon() const
Definition: SuiteSparseQRSupport.h:240
Eigen::SPQR::compute
void compute(const _MatrixType &matrix)
Definition: SuiteSparseQRSupport.h:125
Eigen::NumericalIssue
@ NumericalIssue
Definition: Constants.h:444
Eigen::SPQR::PermutationType
Map< PermutationMatrix< Dynamic, Dynamic, StorageIndex > > PermutationType
Definition: SuiteSparseQRSupport.h:70
Eigen::SPQR_QProduct::rows
Index rows() const
Definition: SuiteSparseQRSupport.h:282
Eigen
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
Eigen::SPQR_QProduct::m_transpose
bool m_transpose
Definition: SuiteSparseQRSupport.h:299
Eigen::SparseMatrix< Scalar, ColMajor, StorageIndex >
Eigen::DenseBase::setZero
EIGEN_DEVICE_FUNC Derived & setZero()
Definition: CwiseNullaryOp.h:546
Eigen::ReturnByValue
Definition: ReturnByValue.h:50
Eigen::SPQR_QProduct::SPQR_QProduct
SPQR_QProduct(const SPQRType &spqr, const Derived &other, bool transpose)
Definition: SuiteSparseQRSupport.h:280
col
m col(1)
Eigen::SPQR::rank
Index rank() const
Definition: SuiteSparseQRSupport.h:225
MatrixType
MatrixXf MatrixType
Definition: benchmark-blocking-sizes.cpp:52
Eigen::internal::traits< SPQRMatrixQReturnType< SPQRType > >::ReturnType
SPQRType::MatrixType ReturnType
Definition: SuiteSparseQRSupport.h:23
b
Scalar * b
Definition: benchVecAdd.cpp:17
Eigen::SPQR::matrixQ
SPQRMatrixQReturnType< SPQR > matrixQ() const
Get an expression of the matrix Q.
Definition: SuiteSparseQRSupport.h:211
Eigen::SPQR::m_factorizationIsOk
bool m_factorizationIsOk
Definition: SuiteSparseQRSupport.h:255
eigen_assert
#define eigen_assert(x)
Definition: Macros.h:1037
Eigen::SPQR::m_info
ComputationInfo m_info
Definition: SuiteSparseQRSupport.h:257
Eigen::SPQRMatrixQTransposeReturnType::m_spqr
const SPQRType & m_spqr
Definition: SuiteSparseQRSupport.h:331
Eigen::SPQRMatrixQReturnType::operator*
SPQR_QProduct< SPQRType, Derived > operator*(const MatrixBase< Derived > &other)
Definition: SuiteSparseQRSupport.h:307
Eigen::SPQR::SPQR
SPQR()
Definition: SuiteSparseQRSupport.h:76
Eigen::Success
@ Success
Definition: Constants.h:442
Eigen::SPQR::setSPQROrdering
void setSPQROrdering(int ord)
Set the fill-reducing ordering method to be used.
Definition: SuiteSparseQRSupport.h:231
Eigen::SPQR::SPQR_free
void SPQR_free()
Definition: SuiteSparseQRSupport.h:116
Eigen::SPQR::m_R
MatrixType m_R
Definition: SuiteSparseQRSupport.h:262
Eigen::SPQR_QProduct::m_other
const Derived & m_other
Definition: SuiteSparseQRSupport.h:298
mat
MatrixXf mat
Definition: Tutorial_AdvancedInitialization_CommaTemporary.cpp:1
res
cout<< "Here is the matrix m:"<< endl<< m<< endl;Matrix< ptrdiff_t, 3, 1 > res
Definition: PartialRedux_count.cpp:3
Eigen::internal::traits< SPQR_QProduct< SPQRType, Derived > >::ReturnType
Derived::PlainObject ReturnType
Definition: SuiteSparseQRSupport.h:31
Eigen::SPQRMatrixQTransposeReturnType::operator*
SPQR_QProduct< SPQRType, Derived > operator*(const MatrixBase< Derived > &other)
Definition: SuiteSparseQRSupport.h:327
Eigen::SPQRMatrixQReturnType::transpose
SPQRMatrixQTransposeReturnType< SPQRType > transpose() const
Definition: SuiteSparseQRSupport.h:316
Eigen::SPQR::setPivotThreshold
void setPivotThreshold(const RealScalar &tol)
Set the tolerance tol to treat columns with 2-norm < =tol as zero.
Definition: SuiteSparseQRSupport.h:233
Eigen::viewAsCholmod
cholmod_sparse viewAsCholmod(Ref< SparseMatrix< _Scalar, _Options, _StorageIndex > > mat)
Definition: CholmodSupport.h:58
Eigen::SPQR::MatrixType
SparseMatrix< Scalar, ColMajor, StorageIndex > MatrixType
Definition: SuiteSparseQRSupport.h:69
Eigen::SPQRMatrixQTransposeReturnType
Definition: SuiteSparseQRSupport.h:18
epsilon
static double epsilon
Definition: testRot3.cpp:37
Eigen::SPQR::colsPermutation
PermutationType colsPermutation() const
Get the permutation that was applied to columns of A.
Definition: SuiteSparseQRSupport.h:216
A
Definition: test_numpy_dtypes.cpp:298
Eigen::SPQR::m_rows
Index m_rows
Definition: SuiteSparseQRSupport.h:270
Eigen::SPQR::m_rank
Index m_rank
Definition: SuiteSparseQRSupport.h:267
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2
Eigen::Dynamic
const int Dynamic
Definition: Constants.h:22
Eigen::SPQR::Scalar
_MatrixType::Scalar Scalar
Definition: SuiteSparseQRSupport.h:66
Eigen::SPQR::RealScalar
_MatrixType::RealScalar RealScalar
Definition: SuiteSparseQRSupport.h:67
Eigen::SPQR::Base
SparseSolverBase< SPQR< _MatrixType > > Base
Definition: SuiteSparseQRSupport.h:63
Eigen::SPQR_QProduct::m_spqr
const SPQRType & m_spqr
Definition: SuiteSparseQRSupport.h:297
Eigen::SPQR::m_useDefaultThreshold
bool m_useDefaultThreshold
Definition: SuiteSparseQRSupport.h:269
Eigen::SPQR::info
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SuiteSparseQRSupport.h:248
Eigen::Map
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
Eigen::SPQR::m_isInitialized
bool m_isInitialized
Definition: SparseSolverBase.h:119
SuiteSparse_long
#define SuiteSparse_long
Definition: SuiteSparse_config.h:62
Eigen::SPQR::rows
Index rows() const
Definition: SuiteSparseQRSupport.h:164
Eigen::PlainObjectBase< Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > >::Map
static ConstMapType Map(const Scalar *data)
Definition: PlainObjectBase.h:644
y
Scalar * y
Definition: level1_cplx_impl.h:124
matrix
Map< Matrix< T, Dynamic, Dynamic, ColMajor >, 0, OuterStride<> > matrix(T *data, int rows, int cols, int stride)
Definition: gtsam/3rdparty/Eigen/blas/common.h:110
Eigen::SPQR::m_cR
cholmod_sparse * m_cR
Definition: SuiteSparseQRSupport.h:261
Eigen::SPQR::cols
Index cols() const
Definition: SuiteSparseQRSupport.h:169
Eigen::SPQR::matrixR
const MatrixType matrixR() const
Definition: SuiteSparseQRSupport.h:201
Eigen::SPQR_QProduct::StorageIndex
SPQRType::StorageIndex StorageIndex
Definition: SuiteSparseQRSupport.h:278
RealScalar
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:47
Eigen::SPQR::ColsAtCompileTime
@ ColsAtCompileTime
Definition: SuiteSparseQRSupport.h:72
Eigen::SPQR::m_HPinv
StorageIndex * m_HPinv
Definition: SuiteSparseQRSupport.h:265
Eigen::SPQR
Sparse QR factorization based on SuiteSparseQR library.
Definition: SuiteSparseQRSupport.h:16
Eigen::SparseSolverBase
A base class for sparse solvers.
Definition: SparseSolverBase.h:67
Eigen::SPQRMatrixQReturnType::adjoint
SPQRMatrixQTransposeReturnType< SPQRType > adjoint() const
Definition: SuiteSparseQRSupport.h:311
Eigen::internal::traits< SPQRMatrixQTransposeReturnType< SPQRType > >::ReturnType
SPQRType::MatrixType ReturnType
Definition: SuiteSparseQRSupport.h:27
Eigen::SPQR::m_tolerance
RealScalar m_tolerance
Definition: SuiteSparseQRSupport.h:260
Eigen::internal::traits
Definition: ForwardDeclarations.h:17
Eigen::SPQR::m_cc
cholmod_common m_cc
Definition: SuiteSparseQRSupport.h:268
Eigen::SPQR::m_H
cholmod_sparse * m_H
Definition: SuiteSparseQRSupport.h:264
Eigen::SPQR::m_E
StorageIndex * m_E
Definition: SuiteSparseQRSupport.h:263
Eigen::SPQRMatrixQTransposeReturnType::SPQRMatrixQTransposeReturnType
SPQRMatrixQTransposeReturnType(const SPQRType &spqr)
Definition: SuiteSparseQRSupport.h:325
gtsam::tol
const G double tol
Definition: Group.h:79
Eigen::SPQR::m_analysisIsOk
bool m_analysisIsOk
Definition: SuiteSparseQRSupport.h:254
Eigen::SPQR::SPQR
SPQR(const _MatrixType &matrix)
Definition: SuiteSparseQRSupport.h:93
Eigen::SPQR::m_isRUpToDate
bool m_isRUpToDate
Definition: SuiteSparseQRSupport.h:256
internal
Definition: BandTriangularSolver.h:13
Eigen::numext::maxi
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T maxi(const T &x, const T &y)
Definition: Eigen/src/Core/MathFunctions.h:1093
Eigen::MatrixBase
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
Eigen::SPQR_QProduct::cols
Index cols() const
Definition: SuiteSparseQRSupport.h:283
Eigen::SPQRMatrixQReturnType
Definition: SuiteSparseQRSupport.h:17
Eigen::ComputationInfo
ComputationInfo
Definition: Constants.h:440
max
#define max(a, b)
Definition: datatypes.h:20
Eigen::SPQR::StorageIndex
SuiteSparse_long StorageIndex
Definition: SuiteSparseQRSupport.h:68
Eigen::SPQR::MaxColsAtCompileTime
@ MaxColsAtCompileTime
Definition: SuiteSparseQRSupport.h:73
Eigen::SPQRMatrixQReturnType::m_spqr
const SPQRType & m_spqr
Definition: SuiteSparseQRSupport.h:320
Eigen::SPQR_QProduct::evalTo
void evalTo(ResType &res) const
Definition: SuiteSparseQRSupport.h:286
Eigen::SPQR::m_HTau
cholmod_dense * m_HTau
Definition: SuiteSparseQRSupport.h:266
Eigen::SPQR::m_allow_tol
int m_allow_tol
Definition: SuiteSparseQRSupport.h:259
Eigen::SPQR_QProduct
Definition: SuiteSparseQRSupport.h:19
Eigen::SPQRMatrixQReturnType::SPQRMatrixQReturnType
SPQRMatrixQReturnType(const SPQRType &spqr)
Definition: SuiteSparseQRSupport.h:305
Eigen::SPQR_QProduct::Scalar
SPQRType::Scalar Scalar
Definition: SuiteSparseQRSupport.h:277
Eigen::NumTraits
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:232
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
pybind_wrapper_test_script.other
other
Definition: pybind_wrapper_test_script.py:42
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:46
Eigen::SPQR::_solve_impl
void _solve_impl(const MatrixBase< Rhs > &b, MatrixBase< Dest > &dest) const
Definition: SuiteSparseQRSupport.h:172
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
Eigen::SPQR::m_ordering
int m_ordering
Definition: SuiteSparseQRSupport.h:258
Eigen::SparseSolverBase::m_isInitialized
bool m_isInitialized
Definition: SparseSolverBase.h:119


gtsam
Author(s):
autogenerated on Sat Nov 16 2024 04:05:02