10 #ifndef EIGEN_SUITESPARSEQRSUPPORT_H
11 #define EIGEN_SUITESPARSEQRSUPPORT_H
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;
20 template <
typename SPQRType>
struct traits<SPQRMatrixQReturnType<SPQRType> >
22 typedef typename SPQRType::MatrixType ReturnType;
24 template <
typename SPQRType>
struct traits<SPQRMatrixQTransposeReturnType<SPQRType> >
26 typedef typename SPQRType::MatrixType ReturnType;
28 template <
typename SPQRType,
typename Derived>
struct traits<SPQR_QProduct<SPQRType, Derived> >
30 typedef typename Derived::PlainObject ReturnType;
56 template<
typename _MatrixType>
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;
67 : m_isInitialized(false),
68 m_ordering(SPQR_ORDERING_DEFAULT),
69 m_allow_tol(SPQR_DEFAULT_TOL),
70 m_tolerance (NumTraits<Scalar>::epsilon())
72 cholmod_l_start(&m_cc);
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())
81 cholmod_l_start(&m_cc);
88 cholmod_l_finish(&m_cc);
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);
99 void compute(
const _MatrixType& matrix)
101 if(m_isInitialized) SPQR_free();
103 MatrixType mat(matrix);
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);
113 m_isInitialized =
false;
117 m_isInitialized =
true;
118 m_isRUpToDate =
false;
123 inline Index
rows()
const {
return m_H->nrow; }
128 inline Index
cols()
const {
return m_cR->ncol; }
134 template<
typename Rhs>
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());
143 template<
typename Rhs,
typename Dest>
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");
150 typename Dest::PlainObject y;
153 Index rk = this->
rank();
155 y.bottomRows(
cols()-rk).setZero();
166 eigen_assert(m_isInitialized &&
" The QR factorization should be computed first, call compute()");
168 m_R = viewAsEigen<Scalar,ColMajor, typename MatrixType::Index>(*m_cR);
169 m_isRUpToDate =
true;
176 return SPQRMatrixQReturnType<SPQR>(*this);
181 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
182 Index n = m_cR->ncol;
184 for(Index j = 0; j <n; j++) colsPerm.
indices()(j) = m_E[j];
194 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
195 return m_cc.SPQR_istat[4];
213 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
217 bool m_isInitialized;
219 bool m_factorizationIsOk;
220 mutable bool m_isRUpToDate;
224 RealScalar m_tolerance;
225 mutable cholmod_sparse *m_cR;
226 mutable MatrixType m_R;
228 mutable cholmod_sparse *m_H;
229 mutable Index *m_HPinv;
230 mutable cholmod_dense *m_HTau;
231 mutable Index m_rank;
232 mutable cholmod_common m_cc;
233 template<
typename ,
typename >
friend struct SPQR_QProduct;
236 template <
typename SPQRType,
typename Derived>
237 struct SPQR_QProduct : ReturnByValue<SPQR_QProduct<SPQRType,Derived> >
239 typedef typename SPQRType::Scalar Scalar;
240 typedef typename SPQRType::Index Index;
242 SPQR_QProduct(
const SPQRType& spqr,
const Derived& other,
bool transpose) : m_spqr(spqr),m_other(other),m_transpose(transpose) {}
244 inline Index rows()
const {
return m_transpose ? m_spqr.rows() : m_spqr.cols(); }
245 inline Index cols()
const {
return m_other.cols(); }
247 template<
typename ResType>
248 void evalTo(ResType& res)
const
252 int method = m_transpose ? SPQR_QTX : SPQR_QX;
253 cholmod_common *cc = m_spqr.cholmodCommon();
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);
259 const SPQRType& m_spqr;
260 const Derived& m_other;
264 template<
typename SPQRType>
265 struct SPQRMatrixQReturnType{
267 SPQRMatrixQReturnType(
const SPQRType& spqr) : m_spqr(spqr) {}
268 template<
typename Derived>
269 SPQR_QProduct<SPQRType, Derived>
operator*(
const MatrixBase<Derived>& other)
271 return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(),
false);
273 SPQRMatrixQTransposeReturnType<SPQRType> adjoint()
const
275 return SPQRMatrixQTransposeReturnType<SPQRType>(m_spqr);
278 SPQRMatrixQTransposeReturnType<SPQRType> transpose()
const
280 return SPQRMatrixQTransposeReturnType<SPQRType>(m_spqr);
282 const SPQRType& m_spqr;
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)
291 return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(),
true);
293 const SPQRType& m_spqr;
298 template<
typename _MatrixType,
typename Rhs>
299 struct solve_retval<SPQR<_MatrixType>, Rhs>
300 : solve_retval_base<SPQR<_MatrixType>, Rhs>
302 typedef SPQR<_MatrixType> Dec;
303 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
305 template<typename Dest>
void evalTo(Dest& dst)
const
307 dec()._solve(rhs(),dst);
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