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 ;
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);
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);
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;
151 y = matrixQ().transpose() * b;
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();
157 dest.
topRows(cols()) = colsPermutation() * y.topRows(cols());
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;
183 PermutationType colsPerm(n);
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();
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);
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>
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>
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>
303 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
305 template<
typename Dest>
void evalTo(Dest& dst)
const 307 dec()._solve(rhs(),dst);
Index cols() const
Definition: SuiteSparseQRSupport.h:128
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SuiteSparseQRSupport.h:211
const internal::solve_retval< SPQR, Rhs > solve(const MatrixBase< Rhs > &B) const
Definition: SuiteSparseQRSupport.h:135
Definition: Constants.h:378
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:88
Index rank() const
Definition: SuiteSparseQRSupport.h:192
Permutation matrix.
Definition: PermutationMatrix.h:283
const IndicesType & indices() const
Definition: PermutationMatrix.h:358
const MatrixType matrixR() const
Definition: SuiteSparseQRSupport.h:164
void setSPQROrdering(int ord)
Set the fill-reducing ordering method to be used.
Definition: SuiteSparseQRSupport.h:198
PermutationType colsPermutation() const
Get the permutation that was applied to columns of A.
Definition: SuiteSparseQRSupport.h:179
SPQRMatrixQReturnType< SPQR > matrixQ() const
Get an expression of the matrix Q.
Definition: SuiteSparseQRSupport.h:174
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
cholmod_common * cholmodCommon() const
Definition: SuiteSparseQRSupport.h:203
Definition: Eigen_Colamd.h:54
Index rows() const
Definition: SuiteSparseQRSupport.h:123
Sparse QR factorization based on SuiteSparseQR library.
Definition: SuiteSparseQRSupport.h:15
Definition: Constants.h:376
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:127
ComputationInfo
Definition: Constants.h:374
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48