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 //
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 ;
65  public:
66  SPQR()
67  : m_ordering(SPQR_ORDERING_DEFAULT),
68  m_allow_tol(SPQR_DEFAULT_TOL),
69  m_tolerance (NumTraits<Scalar>::epsilon())
70  {
71  cholmod_l_start(&m_cc);
72  }
73 
74  SPQR(const _MatrixType& matrix)
75  : m_ordering(SPQR_ORDERING_DEFAULT),
76  m_allow_tol(SPQR_DEFAULT_TOL),
77  m_tolerance (NumTraits<Scalar>::epsilon())
78  {
79  cholmod_l_start(&m_cc);
80  compute(matrix);
81  }
82 
84  {
85  // Calls SuiteSparseQR_free()
86  cholmod_l_free_sparse(&m_H, &m_cc);
87  cholmod_l_free_sparse(&m_cR, &m_cc);
88  cholmod_l_free_dense(&m_HTau, &m_cc);
89  std::free(m_E);
90  std::free(m_HPinv);
91  cholmod_l_finish(&m_cc);
92  }
93  void compute(const _MatrixType& matrix)
94  {
95  MatrixType mat(matrix);
96  cholmod_sparse A;
97  A = viewAsCholmod(mat);
98  Index col = matrix.cols();
99  m_rank = SuiteSparseQR<Scalar>(m_ordering, m_tolerance, col, &A,
100  &m_cR, &m_E, &m_H, &m_HPinv, &m_HTau, &m_cc);
101 
102  if (!m_cR)
103  {
104  m_info = NumericalIssue;
105  m_isInitialized = false;
106  return;
107  }
108  m_info = Success;
109  m_isInitialized = true;
110  m_isRUpToDate = false;
111  }
115  inline Index rows() const {return m_H->nrow; }
116 
120  inline Index cols() const { return m_cR->ncol; }
121 
126  template<typename Rhs>
128  {
129  eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()");
130  eigen_assert(this->rows()==B.rows()
131  && "SPQR::solve(): invalid number of rows of the right hand side matrix B");
132  return internal::solve_retval<SPQR, Rhs>(*this, B.derived());
133  }
134 
135  template<typename Rhs, typename Dest>
136  void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
137  {
138  eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()");
139  eigen_assert(b.cols()==1 && "This method is for vectors only");
140 
141  //Compute Q^T * b
142  Dest y;
143  y = matrixQ().transpose() * b;
144  // Solves with the triangular matrix R
145  Index rk = this->rank();
146  y.topRows(rk) = this->matrixR().topLeftCorner(rk, rk).template triangularView<Upper>().solve(y.topRows(rk));
147  y.bottomRows(cols()-rk).setZero();
148  // Apply the column permutation
149  dest.topRows(cols()) = colsPermutation() * y.topRows(cols());
150 
151  m_info = Success;
152  }
153 
156  const MatrixType matrixR() const
157  {
158  eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()");
159  if(!m_isRUpToDate) {
160  m_R = viewAsEigen<Scalar,ColMajor, typename MatrixType::Index>(*m_cR);
161  m_isRUpToDate = true;
162  }
163  return m_R;
164  }
167  {
168  return SPQRMatrixQReturnType<SPQR>(*this);
169  }
171  PermutationType colsPermutation() const
172  {
173  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
174  Index n = m_cR->ncol;
175  PermutationType colsPerm(n);
176  for(Index j = 0; j <n; j++) colsPerm.indices()(j) = m_E[j];
177  return colsPerm;
178 
179  }
184  Index rank() const
185  {
186  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
187  return m_cc.SPQR_istat[4];
188  }
190  void setSPQROrdering(int ord) { m_ordering = ord;}
192  void setPivotThreshold(const RealScalar& tol) { m_tolerance = tol; }
193 
195  cholmod_common *cholmodCommon() const { return &m_cc; }
196 
197 
204  {
205  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
206  return m_info;
207  }
208  protected:
212  mutable bool m_isRUpToDate;
214  int m_ordering; // Ordering method to use, see SPQR's manual
215  int m_allow_tol; // Allow to use some tolerance during numerical factorization.
216  RealScalar m_tolerance; // treat columns with 2-norm below this tolerance as zero
217  mutable cholmod_sparse *m_cR; // The sparse R factor in cholmod format
218  mutable MatrixType m_R; // The sparse matrix R in Eigen format
219  mutable Index *m_E; // The permutation applied to columns
220  mutable cholmod_sparse *m_H; //The householder vectors
221  mutable Index *m_HPinv; // The row permutation of H
222  mutable cholmod_dense *m_HTau; // The Householder coefficients
223  mutable Index m_rank; // The rank of the matrix
224  mutable cholmod_common m_cc; // Workspace and parameters
225  template<typename ,typename > friend struct SPQR_QProduct;
226 };
227 
228 template <typename SPQRType, typename Derived>
229 struct SPQR_QProduct : ReturnByValue<SPQR_QProduct<SPQRType,Derived> >
230 {
231  typedef typename SPQRType::Scalar Scalar;
232  typedef typename SPQRType::Index Index;
233  //Define the constructor to get reference to argument types
234  SPQR_QProduct(const SPQRType& spqr, const Derived& other, bool transpose) : m_spqr(spqr),m_other(other),m_transpose(transpose) {}
235 
236  inline Index rows() const { return m_transpose ? m_spqr.rows() : m_spqr.cols(); }
237  inline Index cols() const { return m_other.cols(); }
238  // Assign to a vector
239  template<typename ResType>
240  void evalTo(ResType& res) const
241  {
242  cholmod_dense y_cd;
243  cholmod_dense *x_cd;
244  int method = m_transpose ? SPQR_QTX : SPQR_QX;
245  cholmod_common *cc = m_spqr.cholmodCommon();
246  y_cd = viewAsCholmod(m_other.const_cast_derived());
247  x_cd = SuiteSparseQR_qmult<Scalar>(method, m_spqr.m_H, m_spqr.m_HTau, m_spqr.m_HPinv, &y_cd, cc);
248  res = Matrix<Scalar,ResType::RowsAtCompileTime,ResType::ColsAtCompileTime>::Map(reinterpret_cast<Scalar*>(x_cd->x), x_cd->nrow, x_cd->ncol);
249  cholmod_l_free_dense(&x_cd, cc);
250  }
251  const SPQRType& m_spqr;
252  const Derived& m_other;
253  bool m_transpose;
254 
255 };
256 template<typename SPQRType>
257 struct SPQRMatrixQReturnType{
258 
259  SPQRMatrixQReturnType(const SPQRType& spqr) : m_spqr(spqr) {}
260  template<typename Derived>
262  {
263  return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(),false);
264  }
266  {
268  }
269  // To use for operations with the transpose of Q
271  {
273  }
274  const SPQRType& m_spqr;
275 };
276 
277 template<typename SPQRType>
279  SPQRMatrixQTransposeReturnType(const SPQRType& spqr) : m_spqr(spqr) {}
280  template<typename Derived>
282  {
283  return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(), true);
284  }
285  const SPQRType& m_spqr;
286 };
287 
288 namespace internal {
289 
290 template<typename _MatrixType, typename Rhs>
291 struct solve_retval<SPQR<_MatrixType>, Rhs>
292  : solve_retval_base<SPQR<_MatrixType>, Rhs>
293 {
295  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
296 
297  template<typename Dest> void evalTo(Dest& dst) const
298  {
299  dec()._solve(rhs(),dst);
300  }
301 };
302 
303 } // end namespace internal
304 
305 }// End namespace Eigen
306 #endif
SPQR_QProduct< SPQRType, Derived > operator*(const MatrixBase< Derived > &other)
cholmod_common * cholmodCommon() const
SparseMatrix< Scalar, ColMajor, Index > MatrixType
PermutationMatrix< Dynamic, Dynamic > PermutationType
SPQRMatrixQTransposeReturnType< SPQRType > transpose() const
void _solve(const MatrixBase< Rhs > &b, MatrixBase< Dest > &dest) const
RowsBlockXpr topRows(Index n)
Definition: DenseBase.h:381
SPQR(const _MatrixType &matrix)
SPQRMatrixQReturnType(const SPQRType &spqr)
Index rows() const
ComputationInfo m_info
iterative scaling algorithm to equilibrate rows and column norms in matrices
Definition: matrix.hpp:471
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:88
_MatrixType::Scalar Scalar
const internal::solve_retval< SPQR, Rhs > solve(const MatrixBase< Rhs > &B) const
RealScalar m_tolerance
Index cols() const
cholmod_common m_cc
void setPivotThreshold(const RealScalar &tol)
Set the tolerance tol to treat columns with 2-norm < =tol as zero.
cholmod_dense * m_HTau
cholmod_sparse viewAsCholmod(SparseMatrix< _Scalar, _Options, _Index > &mat)
PermutationType colsPermutation() const
Get the permutation that was applied to columns of A.
ComputationInfo info() const
Reports whether previous computation was successful.
SPQR_QProduct(const SPQRType &spqr, const Derived &other, bool transpose)
Index rank() const
cholmod_sparse * m_cR
void compute(const _MatrixType &matrix)
void rhs(const real_t *x, real_t *f)
const MatrixType matrixR() const
cholmod_sparse * m_H
void setSPQROrdering(int ord)
Set the fill-reducing ordering method to be used.
const IndicesType & indices() const
SPQRMatrixQReturnType< SPQR > matrixQ() const
Get an expression of the matrix Q.
SPQRMatrixQTransposeReturnType(const SPQRType &spqr)
Sparse QR factorization based on SuiteSparseQR library.
void evalTo(ResType &res) const
#define EIGEN_MAKE_SOLVE_HELPERS(DecompositionType, Rhs)
Definition: Solve.h:61
SPQRMatrixQTransposeReturnType< SPQRType > adjoint() const
ColXpr col(Index i)
Definition: BlockMethods.h:708
#define eigen_assert(x)
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
_MatrixType::RealScalar RealScalar


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Mon Jun 10 2019 12:35:11