11 #ifndef EIGEN_SUITESPARSEQRSUPPORT_H
12 #define EIGEN_SUITESPARSEQRSUPPORT_H
14 #include "./InternalHeaderCheck.h"
18 template<
typename MatrixType>
class SPQR;
19 template<
typename SPQRType>
struct SPQRMatrixQReturnType;
20 template<
typename SPQRType>
struct SPQRMatrixQTransposeReturnType;
21 template <
typename SPQRType,
typename Derived>
struct SPQR_QProduct;
23 template <
typename SPQRType>
struct traits<SPQRMatrixQReturnType<SPQRType> >
25 typedef typename SPQRType::MatrixType ReturnType;
27 template <
typename SPQRType>
struct traits<SPQRMatrixQTransposeReturnType<SPQRType> >
29 typedef typename SPQRType::MatrixType ReturnType;
31 template <
typename SPQRType,
typename Derived>
struct traits<SPQR_QProduct<SPQRType, Derived> >
33 typedef typename Derived::PlainObject ReturnType;
61 template<
typename MatrixType_>
66 using Base::m_isInitialized;
68 typedef typename MatrixType_::Scalar Scalar;
69 typedef typename MatrixType_::RealScalar RealScalar;
70 typedef SuiteSparse_long StorageIndex ;
79 : m_analysisIsOk(
false),
80 m_factorizationIsOk(
false),
82 m_ordering(SPQR_ORDERING_DEFAULT),
83 m_allow_tol(SPQR_DEFAULT_TOL),
90 m_useDefaultThreshold(
true)
92 cholmod_l_start(&m_cc);
95 explicit SPQR(
const MatrixType_& matrix)
96 : m_analysisIsOk(
false),
97 m_factorizationIsOk(
false),
99 m_ordering(SPQR_ORDERING_DEFAULT),
100 m_allow_tol(SPQR_DEFAULT_TOL),
107 m_useDefaultThreshold(
true)
109 cholmod_l_start(&m_cc);
116 cholmod_l_finish(&m_cc);
120 cholmod_l_free_sparse(&m_H, &m_cc);
121 cholmod_l_free_sparse(&m_cR, &m_cc);
122 cholmod_l_free_dense(&m_HTau, &m_cc);
127 void compute(
const MatrixType_& matrix)
129 if(m_isInitialized) SPQR_free();
137 RealScalar pivotThreshold = m_tolerance;
138 if(m_useDefaultThreshold)
140 RealScalar max2Norm = 0.0;
141 for (
int j = 0; j < mat.
cols(); j++) max2Norm = numext::maxi(max2Norm, mat.col(j).norm());
142 if(numext::is_exactly_zero(max2Norm))
143 max2Norm = RealScalar(1);
148 m_rows = matrix.rows();
149 Index col = matrix.cols();
150 m_rank = SuiteSparseQR<Scalar>(m_ordering, pivotThreshold, col, &A,
151 &m_cR, &m_E, &m_H, &m_HPinv, &m_HTau, &m_cc);
156 m_isInitialized =
false;
160 m_isInitialized =
true;
161 m_isRUpToDate =
false;
173 template<
typename Rhs,
typename Dest>
176 eigen_assert(m_isInitialized &&
" The QR factorization should be computed first, call compute()");
177 eigen_assert(b.
cols()==1 &&
"This method is for vectors only");
180 typename Dest::PlainObject y, y2;
186 y.resize((std::max)(
cols(),
Index(y.rows())),y.cols());
187 y.topRows(rk) = this->
matrixR().topLeftCorner(rk, rk).template triangularView<Upper>().solve(y2.topRows(rk));
192 for(
Index i = 0; i < rk; ++i) dest.row(m_E[i]) = y.row(i);
205 eigen_assert(m_isInitialized &&
" The QR factorization should be computed first, call compute()");
207 m_R = viewAsEigen<Scalar,ColMajor, typename MatrixType::StorageIndex>(*m_cR);
208 m_isRUpToDate =
true;
215 return SPQRMatrixQReturnType<SPQR>(*
this);
220 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
229 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
230 return m_cc.SPQR_istat[4];
237 m_useDefaultThreshold =
false;
252 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
257 bool m_factorizationIsOk;
258 mutable bool m_isRUpToDate;
262 RealScalar m_tolerance;
263 mutable cholmod_sparse *m_cR =
nullptr;
264 mutable MatrixType m_R;
265 mutable StorageIndex *m_E =
nullptr;
266 mutable cholmod_sparse *m_H =
nullptr;
267 mutable StorageIndex *m_HPinv =
nullptr;
268 mutable cholmod_dense *m_HTau =
nullptr;
269 mutable Index m_rank;
270 mutable cholmod_common m_cc;
271 bool m_useDefaultThreshold;
273 template<
typename ,
typename >
friend struct SPQR_QProduct;
276 template <
typename SPQRType,
typename Derived>
277 struct SPQR_QProduct : ReturnByValue<SPQR_QProduct<SPQRType,Derived> >
279 typedef typename SPQRType::Scalar Scalar;
280 typedef typename SPQRType::StorageIndex StorageIndex;
282 SPQR_QProduct(
const SPQRType& spqr,
const Derived& other,
bool transpose) : m_spqr(spqr),m_other(other),m_transpose(transpose) {}
284 inline Index rows()
const {
return m_transpose ? m_spqr.rows() : m_spqr.cols(); }
285 inline Index cols()
const {
return m_other.cols(); }
287 template<
typename ResType>
288 void evalTo(ResType& res)
const
292 int method = m_transpose ? SPQR_QTX : SPQR_QX;
293 cholmod_common *cc = m_spqr.cholmodCommon();
295 x_cd = SuiteSparseQR_qmult<Scalar>(method, m_spqr.m_H, m_spqr.m_HTau, m_spqr.m_HPinv, &y_cd, cc);
297 cholmod_l_free_dense(&x_cd, cc);
299 const SPQRType& m_spqr;
300 const Derived& m_other;
304 template<
typename SPQRType>
305 struct SPQRMatrixQReturnType{
307 SPQRMatrixQReturnType(
const SPQRType& spqr) : m_spqr(spqr) {}
308 template<
typename Derived>
309 SPQR_QProduct<SPQRType, Derived>
operator*(
const MatrixBase<Derived>& other)
311 return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(),
false);
313 SPQRMatrixQTransposeReturnType<SPQRType> adjoint()
const
315 return SPQRMatrixQTransposeReturnType<SPQRType>(m_spqr);
318 SPQRMatrixQTransposeReturnType<SPQRType> transpose()
const
320 return SPQRMatrixQTransposeReturnType<SPQRType>(m_spqr);
322 const SPQRType& m_spqr;
325 template<
typename SPQRType>
326 struct SPQRMatrixQTransposeReturnType{
327 SPQRMatrixQTransposeReturnType(
const SPQRType& spqr) : m_spqr(spqr) {}
328 template<
typename Derived>
329 SPQR_QProduct<SPQRType,Derived>
operator*(
const MatrixBase<Derived>& other)
331 return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(),
true);
333 const SPQRType& m_spqr;
Derived & setZero()
Definition: CwiseNullaryOp.h:548
EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Definition: EigenBase.h:65
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:98
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:52
static ConstMapType Map(const Scalar *data)
Definition: PlainObjectBase.h:642
Sparse QR factorization based on SuiteSparseQR library.
Definition: SuiteSparseQRSupport.h:63
void setPivotThreshold(const RealScalar &tol)
Set the tolerance tol to treat columns with 2-norm < =tol as zero.
Definition: SuiteSparseQRSupport.h:235
Index rank() const
Definition: SuiteSparseQRSupport.h:227
cholmod_common * cholmodCommon() const
Definition: SuiteSparseQRSupport.h:242
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SuiteSparseQRSupport.h:250
Index rows() const
Definition: SuiteSparseQRSupport.h:166
const MatrixType matrixR() const
Definition: SuiteSparseQRSupport.h:203
Index cols() const
Definition: SuiteSparseQRSupport.h:171
SPQRMatrixQReturnType< SPQR > matrixQ() const
Get an expression of the matrix Q.
Definition: SuiteSparseQRSupport.h:213
PermutationType colsPermutation() const
Get the permutation that was applied to columns of A.
Definition: SuiteSparseQRSupport.h:218
void setSPQROrdering(int ord)
Set the fill-reducing ordering method to be used.
Definition: SuiteSparseQRSupport.h:233
Index cols() const
Definition: SparseMatrix.h:142
Index rows() const
Definition: SparseMatrix.h:140
A base class for sparse solvers.
Definition: SparseSolverBase.h:70
ComputationInfo
Definition: Constants.h:442
@ NumericalIssue
Definition: Constants.h:446
@ Success
Definition: Constants.h:444
Namespace containing all symbols from the Eigen library.
Definition: Core:139
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
cholmod_sparse viewAsCholmod(Ref< SparseMatrix< Scalar_, Options_, StorageIndex_ > > mat)
Definition: CholmodSupport.h:60
const Product< MatrixDerived, PermutationDerived, AliasFreeProduct > operator*(const MatrixBase< MatrixDerived > &matrix, const PermutationBase< PermutationDerived > &permutation)
Definition: PermutationMatrix.h:517
const int Dynamic
Definition: Constants.h:24
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:231