Eigen  3.2.2
 All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
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 //
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 
56 template<typename _MatrixType>
57 class SPQR
58 {
59  public:
60  typedef typename _MatrixType::Scalar Scalar;
61  typedef typename _MatrixType::RealScalar RealScalar;
62  typedef UF_long Index ;
63  typedef SparseMatrix<Scalar, ColMajor, Index> MatrixType;
64  typedef PermutationMatrix<Dynamic, Dynamic> PermutationType;
65  public:
66  SPQR()
67  : m_isInitialized(false),
68  m_ordering(SPQR_ORDERING_DEFAULT),
69  m_allow_tol(SPQR_DEFAULT_TOL),
70  m_tolerance (NumTraits<Scalar>::epsilon())
71  {
72  cholmod_l_start(&m_cc);
73  }
74 
75  SPQR(const _MatrixType& matrix)
76  : m_isInitialized(false),
77  m_ordering(SPQR_ORDERING_DEFAULT),
78  m_allow_tol(SPQR_DEFAULT_TOL),
79  m_tolerance (NumTraits<Scalar>::epsilon())
80  {
81  cholmod_l_start(&m_cc);
82  compute(matrix);
83  }
84 
85  ~SPQR()
86  {
87  SPQR_free();
88  cholmod_l_finish(&m_cc);
89  }
90  void SPQR_free()
91  {
92  cholmod_l_free_sparse(&m_H, &m_cc);
93  cholmod_l_free_sparse(&m_cR, &m_cc);
94  cholmod_l_free_dense(&m_HTau, &m_cc);
95  std::free(m_E);
96  std::free(m_HPinv);
97  }
98 
99  void compute(const _MatrixType& matrix)
100  {
101  if(m_isInitialized) SPQR_free();
102 
103  MatrixType mat(matrix);
104  cholmod_sparse A;
105  A = viewAsCholmod(mat);
106  Index col = matrix.cols();
107  m_rank = SuiteSparseQR<Scalar>(m_ordering, m_tolerance, col, &A,
108  &m_cR, &m_E, &m_H, &m_HPinv, &m_HTau, &m_cc);
109 
110  if (!m_cR)
111  {
112  m_info = NumericalIssue;
113  m_isInitialized = false;
114  return;
115  }
116  m_info = Success;
117  m_isInitialized = true;
118  m_isRUpToDate = false;
119  }
123  inline Index rows() const {return m_H->nrow; }
124 
128  inline Index cols() const { return m_cR->ncol; }
129 
134  template<typename Rhs>
135  inline const internal::solve_retval<SPQR, Rhs> solve(const MatrixBase<Rhs>& B) const
136  {
137  eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()");
138  eigen_assert(this->rows()==B.rows()
139  && "SPQR::solve(): invalid number of rows of the right hand side matrix B");
140  return internal::solve_retval<SPQR, Rhs>(*this, B.derived());
141  }
142 
143  template<typename Rhs, typename Dest>
144  void _solve(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;
151  y = matrixQ().transpose() * b;
152  // Solves with the triangular matrix R
153  Index rk = this->rank();
154  y.topRows(rk) = this->matrixR().topLeftCorner(rk, rk).template triangularView<Upper>().solve(y.topRows(rk));
155  y.bottomRows(cols()-rk).setZero();
156  // Apply the column permutation
157  dest.topRows(cols()) = colsPermutation() * y.topRows(cols());
158 
159  m_info = Success;
160  }
161 
164  const MatrixType matrixR() const
165  {
166  eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()");
167  if(!m_isRUpToDate) {
168  m_R = viewAsEigen<Scalar,ColMajor, typename MatrixType::Index>(*m_cR);
169  m_isRUpToDate = true;
170  }
171  return m_R;
172  }
174  SPQRMatrixQReturnType<SPQR> matrixQ() const
175  {
176  return SPQRMatrixQReturnType<SPQR>(*this);
177  }
180  {
181  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
182  Index n = m_cR->ncol;
183  PermutationType colsPerm(n);
184  for(Index j = 0; j <n; j++) colsPerm.indices()(j) = m_E[j];
185  return colsPerm;
186 
187  }
192  Index rank() const
193  {
194  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
195  return m_cc.SPQR_istat[4];
196  }
198  void setSPQROrdering(int ord) { m_ordering = ord;}
200  void setPivotThreshold(const RealScalar& tol) { m_tolerance = tol; }
201 
203  cholmod_common *cholmodCommon() const { return &m_cc; }
204 
205 
212  {
213  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
214  return m_info;
215  }
216  protected:
217  bool m_isInitialized;
218  bool m_analysisIsOk;
219  bool m_factorizationIsOk;
220  mutable bool m_isRUpToDate;
221  mutable ComputationInfo m_info;
222  int m_ordering; // Ordering method to use, see SPQR's manual
223  int m_allow_tol; // Allow to use some tolerance during numerical factorization.
224  RealScalar m_tolerance; // treat columns with 2-norm below this tolerance as zero
225  mutable cholmod_sparse *m_cR; // The sparse R factor in cholmod format
226  mutable MatrixType m_R; // The sparse matrix R in Eigen format
227  mutable Index *m_E; // The permutation applied to columns
228  mutable cholmod_sparse *m_H; //The householder vectors
229  mutable Index *m_HPinv; // The row permutation of H
230  mutable cholmod_dense *m_HTau; // The Householder coefficients
231  mutable Index m_rank; // The rank of the matrix
232  mutable cholmod_common m_cc; // Workspace and parameters
233  template<typename ,typename > friend struct SPQR_QProduct;
234 };
235 
236 template <typename SPQRType, typename Derived>
237 struct SPQR_QProduct : ReturnByValue<SPQR_QProduct<SPQRType,Derived> >
238 {
239  typedef typename SPQRType::Scalar Scalar;
240  typedef typename SPQRType::Index Index;
241  //Define the constructor to get reference to argument types
242  SPQR_QProduct(const SPQRType& spqr, const Derived& other, bool transpose) : m_spqr(spqr),m_other(other),m_transpose(transpose) {}
243 
244  inline Index rows() const { return m_transpose ? m_spqr.rows() : m_spqr.cols(); }
245  inline Index cols() const { return m_other.cols(); }
246  // Assign to a vector
247  template<typename ResType>
248  void evalTo(ResType& res) const
249  {
250  cholmod_dense y_cd;
251  cholmod_dense *x_cd;
252  int method = m_transpose ? SPQR_QTX : SPQR_QX;
253  cholmod_common *cc = m_spqr.cholmodCommon();
254  y_cd = viewAsCholmod(m_other.const_cast_derived());
255  x_cd = SuiteSparseQR_qmult<Scalar>(method, m_spqr.m_H, m_spqr.m_HTau, m_spqr.m_HPinv, &y_cd, cc);
256  res = Matrix<Scalar,ResType::RowsAtCompileTime,ResType::ColsAtCompileTime>::Map(reinterpret_cast<Scalar*>(x_cd->x), x_cd->nrow, x_cd->ncol);
257  cholmod_l_free_dense(&x_cd, cc);
258  }
259  const SPQRType& m_spqr;
260  const Derived& m_other;
261  bool m_transpose;
262 
263 };
264 template<typename SPQRType>
265 struct SPQRMatrixQReturnType{
266 
267  SPQRMatrixQReturnType(const SPQRType& spqr) : m_spqr(spqr) {}
268  template<typename Derived>
269  SPQR_QProduct<SPQRType, Derived> operator*(const MatrixBase<Derived>& other)
270  {
271  return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(),false);
272  }
273  SPQRMatrixQTransposeReturnType<SPQRType> adjoint() const
274  {
275  return SPQRMatrixQTransposeReturnType<SPQRType>(m_spqr);
276  }
277  // To use for operations with the transpose of Q
278  SPQRMatrixQTransposeReturnType<SPQRType> transpose() const
279  {
280  return SPQRMatrixQTransposeReturnType<SPQRType>(m_spqr);
281  }
282  const SPQRType& m_spqr;
283 };
284 
285 template<typename SPQRType>
286 struct SPQRMatrixQTransposeReturnType{
287  SPQRMatrixQTransposeReturnType(const SPQRType& spqr) : m_spqr(spqr) {}
288  template<typename Derived>
289  SPQR_QProduct<SPQRType,Derived> operator*(const MatrixBase<Derived>& other)
290  {
291  return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(), true);
292  }
293  const SPQRType& m_spqr;
294 };
295 
296 namespace internal {
297 
298 template<typename _MatrixType, typename Rhs>
299 struct solve_retval<SPQR<_MatrixType>, Rhs>
300  : solve_retval_base<SPQR<_MatrixType>, Rhs>
301 {
302  typedef SPQR<_MatrixType> Dec;
303  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
304 
305  template<typename Dest> void evalTo(Dest& dst) const
306  {
307  dec()._solve(rhs(),dst);
308  }
309 };
310 
311 } // end namespace internal
312 
313 }// End namespace Eigen
314 #endif
const IndicesType & indices() const
Definition: PermutationMatrix.h:358
const internal::solve_retval< SPQR, Rhs > solve(const MatrixBase< Rhs > &B) const
Definition: SuiteSparseQRSupport.h:135
Definition: Constants.h:378
Index rows() const
Definition: SuiteSparseQRSupport.h:123
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SuiteSparseQRSupport.h:211
const internal::permut_matrix_product_retval< PermutationDerived, Derived, OnTheRight > operator*(const MatrixBase< Derived > &matrix, const PermutationBase< PermutationDerived > &permutation)
Definition: PermutationMatrix.h:510
SPQRMatrixQReturnType< SPQR > matrixQ() const
Get an expression of the matrix Q.
Definition: SuiteSparseQRSupport.h:174
cholmod_common * cholmodCommon() const
Definition: SuiteSparseQRSupport.h:203
Permutation matrix.
Definition: PermutationMatrix.h:283
cholmod_sparse viewAsCholmod(SparseMatrix< _Scalar, _Options, _Index > &mat)
Definition: CholmodSupport.h:52
void setSPQROrdering(int ord)
Set the fill-reducing ordering method to be used.
Definition: SuiteSparseQRSupport.h:198
void setPivotThreshold(const RealScalar &tol)
Set the tolerance tol to treat columns with 2-norm < =tol as zero.
Definition: SuiteSparseQRSupport.h:200
RowsBlockXpr topRows(Index n)
Definition: DenseBase.h:381
Index cols() const
Definition: SuiteSparseQRSupport.h:128
const MatrixType matrixR() const
Definition: SuiteSparseQRSupport.h:164
Sparse QR factorization based on SuiteSparseQR library.
Definition: SuiteSparseQRSupport.h:15
Index rank() const
Definition: SuiteSparseQRSupport.h:192
Definition: Constants.h:376
Block< Derived > topLeftCorner(Index cRows, Index cCols)
Definition: SparseMatrixBase.h:157
ComputationInfo
Definition: Constants.h:374
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
PermutationType colsPermutation() const
Get the permutation that was applied to columns of A.
Definition: SuiteSparseQRSupport.h:179