11 #ifndef EIGEN_SUITESPARSEQRSUPPORT_H
12 #define EIGEN_SUITESPARSEQRSUPPORT_H
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;
21 template <
typename SPQRType>
struct traits<SPQRMatrixQReturnType<SPQRType> >
23 typedef typename SPQRType::MatrixType ReturnType;
25 template <
typename SPQRType>
struct traits<SPQRMatrixQTransposeReturnType<SPQRType> >
27 typedef typename SPQRType::MatrixType ReturnType;
29 template <
typename SPQRType,
typename Derived>
struct traits<SPQR_QProduct<SPQRType, Derived> >
31 typedef typename Derived::PlainObject ReturnType;
57 template<
typename _MatrixType>
58 class SPQR :
public SparseSolverBase<SPQR<_MatrixType> >
61 typedef SparseSolverBase<SPQR<_MatrixType> > Base;
62 using Base::m_isInitialized;
64 typedef typename _MatrixType::Scalar Scalar;
65 typedef typename _MatrixType::RealScalar RealScalar;
66 typedef UF_long StorageIndex ;
67 typedef SparseMatrix<Scalar, ColMajor, StorageIndex> MatrixType;
68 typedef Map<PermutationMatrix<Dynamic, Dynamic, StorageIndex> > PermutationType;
71 : m_ordering(SPQR_ORDERING_DEFAULT), m_allow_tol(SPQR_DEFAULT_TOL), m_tolerance (NumTraits<Scalar>::epsilon()), m_useDefaultThreshold(true)
73 cholmod_l_start(&m_cc);
76 explicit SPQR(
const _MatrixType& matrix)
77 : m_ordering(SPQR_ORDERING_DEFAULT), m_allow_tol(SPQR_DEFAULT_TOL), m_tolerance (NumTraits<Scalar>::epsilon()), m_useDefaultThreshold(true)
79 cholmod_l_start(&m_cc);
86 cholmod_l_finish(&m_cc);
90 cholmod_l_free_sparse(&m_H, &m_cc);
91 cholmod_l_free_sparse(&m_cR, &m_cc);
92 cholmod_l_free_dense(&m_HTau, &m_cc);
97 void compute(
const _MatrixType& matrix)
99 if(m_isInitialized) SPQR_free();
101 MatrixType mat(matrix);
107 RealScalar pivotThreshold = m_tolerance;
108 if(m_useDefaultThreshold)
110 RealScalar max2Norm = 0.0;
111 for (
int j = 0; j < mat.cols(); j++) max2Norm = numext::maxi(max2Norm, mat.col(j).norm());
112 if(max2Norm==RealScalar(0))
113 max2Norm = RealScalar(1);
114 pivotThreshold = 20 * (mat.rows() + mat.cols()) * max2Norm * NumTraits<RealScalar>::epsilon();
118 A = viewAsCholmod(mat);
119 Index col = matrix.cols();
120 m_rank = SuiteSparseQR<Scalar>(m_ordering, pivotThreshold, col, &A,
121 &m_cR, &m_E, &m_H, &m_HPinv, &m_HTau, &m_cc);
126 m_isInitialized =
false;
130 m_isInitialized =
true;
131 m_isRUpToDate =
false;
136 inline Index
rows()
const {
return m_cR->nrow; }
141 inline Index
cols()
const {
return m_cR->ncol; }
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, y2;
154 Index rk = this->
rank();
156 y.resize((std::max)(
cols(),Index(y.rows())),y.cols());
162 for(Index i = 0; i < rk; ++i) dest.
row(m_E[i]) = y.row(i);
163 for(Index i = rk; i <
cols(); ++i) dest.
row(m_E[i]).setZero();
175 eigen_assert(m_isInitialized &&
" The QR factorization should be computed first, call compute()");
177 m_R = viewAsEigen<Scalar,ColMajor, typename MatrixType::StorageIndex>(*m_cR);
178 m_isRUpToDate =
true;
185 return SPQRMatrixQReturnType<SPQR>(*this);
190 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
191 return PermutationType(m_E, m_cR->ncol);
199 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
200 return m_cc.SPQR_istat[4];
207 m_useDefaultThreshold =
false;
222 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
227 bool m_factorizationIsOk;
228 mutable bool m_isRUpToDate;
232 RealScalar m_tolerance;
233 mutable cholmod_sparse *m_cR;
234 mutable MatrixType m_R;
235 mutable StorageIndex *m_E;
236 mutable cholmod_sparse *m_H;
237 mutable StorageIndex *m_HPinv;
238 mutable cholmod_dense *m_HTau;
239 mutable Index m_rank;
240 mutable cholmod_common m_cc;
241 bool m_useDefaultThreshold;
242 template<
typename ,
typename >
friend struct SPQR_QProduct;
245 template <
typename SPQRType,
typename Derived>
246 struct SPQR_QProduct : ReturnByValue<SPQR_QProduct<SPQRType,Derived> >
248 typedef typename SPQRType::Scalar Scalar;
249 typedef typename SPQRType::StorageIndex StorageIndex;
251 SPQR_QProduct(
const SPQRType& spqr,
const Derived& other,
bool transpose) : m_spqr(spqr),m_other(other),m_transpose(transpose) {}
253 inline Index rows()
const {
return m_transpose ? m_spqr.rows() : m_spqr.cols(); }
254 inline Index cols()
const {
return m_other.cols(); }
256 template<
typename ResType>
257 void evalTo(ResType& res)
const
261 int method = m_transpose ? SPQR_QTX : SPQR_QX;
262 cholmod_common *cc = m_spqr.cholmodCommon();
263 y_cd = viewAsCholmod(m_other.const_cast_derived());
264 x_cd = SuiteSparseQR_qmult<Scalar>(method, m_spqr.m_H, m_spqr.m_HTau, m_spqr.m_HPinv, &y_cd, cc);
265 res = Matrix<Scalar,ResType::RowsAtCompileTime,ResType::ColsAtCompileTime>::Map(reinterpret_cast<Scalar*>(x_cd->x), x_cd->nrow, x_cd->ncol);
266 cholmod_l_free_dense(&x_cd, cc);
268 const SPQRType& m_spqr;
269 const Derived& m_other;
273 template<
typename SPQRType>
274 struct SPQRMatrixQReturnType{
276 SPQRMatrixQReturnType(
const SPQRType& spqr) : m_spqr(spqr) {}
277 template<
typename Derived>
278 SPQR_QProduct<SPQRType, Derived> operator*(
const MatrixBase<Derived>& other)
280 return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(),
false);
282 SPQRMatrixQTransposeReturnType<SPQRType> adjoint()
const
284 return SPQRMatrixQTransposeReturnType<SPQRType>(m_spqr);
287 SPQRMatrixQTransposeReturnType<SPQRType> transpose()
const
289 return SPQRMatrixQTransposeReturnType<SPQRType>(m_spqr);
291 const SPQRType& m_spqr;
294 template<
typename SPQRType>
295 struct SPQRMatrixQTransposeReturnType{
296 SPQRMatrixQTransposeReturnType(
const SPQRType& spqr) : m_spqr(spqr) {}
297 template<
typename Derived>
298 SPQR_QProduct<SPQRType,Derived> operator*(
const MatrixBase<Derived>& other)
300 return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(),
true);
302 const SPQRType& m_spqr;
cholmod_common * cholmodCommon() const
Definition: SuiteSparseQRSupport.h:212
Index rows() const
Definition: SuiteSparseQRSupport.h:136
const Solve< SPQR< _MatrixType >, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: SparseSolverBase.h:74
RowXpr row(Index i)
Definition: DenseBase.h:797
Index cols() const
Definition: SuiteSparseQRSupport.h:141
void setPivotThreshold(const RealScalar &tol)
Set the tolerance tol to treat columns with 2-norm < =tol as zero.
Definition: SuiteSparseQRSupport.h:205
PermutationType colsPermutation() const
Get the permutation that was applied to columns of A.
Definition: SuiteSparseQRSupport.h:188
Definition: Constants.h:426
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SuiteSparseQRSupport.h:220
Index rank() const
Definition: SuiteSparseQRSupport.h:197
Definition: Constants.h:424
Block< Derived > topLeftCorner(Index cRows, Index cCols)
Definition: SparseMatrixBase.h:164
Definition: Eigen_Colamd.h:54
const MatrixType matrixR() const
Definition: SuiteSparseQRSupport.h:173
void setSPQROrdering(int ord)
Set the fill-reducing ordering method to be used.
Definition: SuiteSparseQRSupport.h:203
SPQRMatrixQReturnType< SPQR > matrixQ() const
Get an expression of the matrix Q.
Definition: SuiteSparseQRSupport.h:183
Sparse QR factorization based on SuiteSparseQR library.
Definition: SuiteSparseQRSupport.h:16
ComputationInfo
Definition: Constants.h:422
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48