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;
59 template<
typename _MatrixType>
64 using Base::m_isInitialized;
66 typedef typename _MatrixType::Scalar Scalar;
67 typedef typename _MatrixType::RealScalar RealScalar;
68 typedef SuiteSparse_long StorageIndex ;
77 : m_analysisIsOk(false),
78 m_factorizationIsOk(false),
80 m_ordering(SPQR_ORDERING_DEFAULT),
81 m_allow_tol(SPQR_DEFAULT_TOL),
82 m_tolerance (
NumTraits<Scalar>::epsilon()),
88 m_useDefaultThreshold(true)
90 cholmod_l_start(&m_cc);
93 explicit SPQR(
const _MatrixType& matrix)
94 : m_analysisIsOk(false),
95 m_factorizationIsOk(false),
97 m_ordering(SPQR_ORDERING_DEFAULT),
98 m_allow_tol(SPQR_DEFAULT_TOL),
99 m_tolerance (
NumTraits<Scalar>::epsilon()),
105 m_useDefaultThreshold(true)
107 cholmod_l_start(&m_cc);
114 cholmod_l_finish(&m_cc);
118 cholmod_l_free_sparse(&m_H, &m_cc);
119 cholmod_l_free_sparse(&m_cR, &m_cc);
120 cholmod_l_free_dense(&m_HTau, &m_cc);
125 void compute(
const _MatrixType& matrix)
127 if(m_isInitialized) SPQR_free();
129 MatrixType mat(matrix);
135 RealScalar pivotThreshold = m_tolerance;
136 if(m_useDefaultThreshold)
138 RealScalar max2Norm = 0.0;
139 for (
int j = 0; j < mat.cols(); j++) max2Norm = numext::maxi(max2Norm, mat.col(j).norm());
140 if(max2Norm==RealScalar(0))
141 max2Norm = RealScalar(1);
146 m_rows = matrix.rows();
147 Index col = matrix.cols();
148 m_rank = SuiteSparseQR<Scalar>(m_ordering, pivotThreshold, col, &A,
149 &m_cR, &m_E, &m_H, &m_HPinv, &m_HTau, &m_cc);
154 m_isInitialized =
false;
158 m_isInitialized =
true;
159 m_isRUpToDate =
false;
171 template<
typename Rhs,
typename Dest>
174 eigen_assert(m_isInitialized &&
" The QR factorization should be computed first, call compute()");
175 eigen_assert(b.
cols()==1 &&
"This method is for vectors only");
178 typename Dest::PlainObject y, y2;
179 y = matrixQ().transpose() * b;
182 Index rk = this->rank();
184 y.resize((std::max)(cols(),
Index(y.rows())),y.cols());
185 y.topRows(rk) = this->matrixR().topLeftCorner(rk, rk).template triangularView<Upper>().solve(y2.topRows(rk));
190 for(
Index i = 0; i < rk; ++i) dest.
row(m_E[i]) = y.row(i);
191 for(
Index i = rk; i < cols(); ++i) dest.
row(m_E[i]).setZero();
203 eigen_assert(m_isInitialized &&
" The QR factorization should be computed first, call compute()");
205 m_R = viewAsEigen<Scalar,ColMajor, typename MatrixType::StorageIndex>(*m_cR);
206 m_isRUpToDate =
true;
213 return SPQRMatrixQReturnType<SPQR>(*this);
218 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
219 return PermutationType(m_E, m_cR->ncol);
227 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
228 return m_cc.SPQR_istat[4];
235 m_useDefaultThreshold =
false;
250 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
255 bool m_factorizationIsOk;
256 mutable bool m_isRUpToDate;
260 RealScalar m_tolerance;
261 mutable cholmod_sparse *m_cR;
262 mutable MatrixType m_R;
263 mutable StorageIndex *m_E;
264 mutable cholmod_sparse *m_H;
265 mutable StorageIndex *m_HPinv;
266 mutable cholmod_dense *m_HTau;
267 mutable Index m_rank;
268 mutable cholmod_common m_cc;
269 bool m_useDefaultThreshold;
271 template<
typename ,
typename >
friend struct SPQR_QProduct;
274 template <
typename SPQRType,
typename Derived>
275 struct SPQR_QProduct : ReturnByValue<SPQR_QProduct<SPQRType,Derived> >
277 typedef typename SPQRType::Scalar Scalar;
278 typedef typename SPQRType::StorageIndex StorageIndex;
280 SPQR_QProduct(
const SPQRType& spqr,
const Derived& other,
bool transpose) : m_spqr(spqr),m_other(other),m_transpose(transpose) {}
282 inline Index rows()
const {
return m_transpose ? m_spqr.rows() : m_spqr.cols(); }
283 inline Index cols()
const {
return m_other.cols(); }
285 template<
typename ResType>
286 void evalTo(ResType& res)
const 290 int method = m_transpose ? SPQR_QTX : SPQR_QX;
291 cholmod_common *cc = m_spqr.cholmodCommon();
293 x_cd = SuiteSparseQR_qmult<Scalar>(method, m_spqr.m_H, m_spqr.m_HTau, m_spqr.m_HPinv, &y_cd, cc);
295 cholmod_l_free_dense(&x_cd, cc);
297 const SPQRType& m_spqr;
298 const Derived& m_other;
302 template<
typename SPQRType>
303 struct SPQRMatrixQReturnType{
305 SPQRMatrixQReturnType(
const SPQRType& spqr) : m_spqr(spqr) {}
306 template<
typename Derived>
309 return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.
derived(),
false);
311 SPQRMatrixQTransposeReturnType<SPQRType> adjoint()
const 313 return SPQRMatrixQTransposeReturnType<SPQRType>(m_spqr);
316 SPQRMatrixQTransposeReturnType<SPQRType> transpose()
const 318 return SPQRMatrixQTransposeReturnType<SPQRType>(m_spqr);
320 const SPQRType& m_spqr;
323 template<
typename SPQRType>
324 struct SPQRMatrixQTransposeReturnType{
325 SPQRMatrixQTransposeReturnType(
const SPQRType& spqr) : m_spqr(spqr) {}
326 template<
typename Derived>
329 return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.
derived(),
true);
331 const SPQRType& m_spqr;
const MatrixType matrixR() const
Definition: SuiteSparseQRSupport.h:201
Index rank() const
Definition: SuiteSparseQRSupport.h:225
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
A base class for sparse solvers.
Definition: SparseSolverBase.h:67
Namespace containing all symbols from the Eigen library.
Definition: Core:141
Index cols() const
Definition: SuiteSparseQRSupport.h:169
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:232
RowXpr row(Index i)
Definition: DenseBase.h:1119
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SuiteSparseQRSupport.h:248
Derived & derived()
Definition: EigenBase.h:46
void setPivotThreshold(const RealScalar &tol)
Set the tolerance tol to treat columns with 2-norm < =tol as zero.
Definition: SuiteSparseQRSupport.h:233
const Product< MatrixDerived, PermutationDerived, AliasFreeProduct > operator*(const MatrixBase< MatrixDerived > &matrix, const PermutationBase< PermutationDerived > &permutation)
Definition: PermutationMatrix.h:515
Definition: Constants.h:444
SPQRMatrixQReturnType< SPQR > matrixQ() const
Get an expression of the matrix Q.
Definition: SuiteSparseQRSupport.h:211
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Definition: EigenBase.h:63
cholmod_sparse viewAsCholmod(Ref< SparseMatrix< _Scalar, _Options, _StorageIndex > > mat)
Definition: CholmodSupport.h:58
Definition: Constants.h:442
Index rows() const
Definition: SuiteSparseQRSupport.h:164
Definition: Eigen_Colamd.h:50
PermutationType colsPermutation() const
Get the permutation that was applied to columns of A.
Definition: SuiteSparseQRSupport.h:216
void setSPQROrdering(int ord)
Set the fill-reducing ordering method to be used.
Definition: SuiteSparseQRSupport.h:231
Sparse QR factorization based on SuiteSparseQR library.
Definition: SuiteSparseQRSupport.h:16
const int Dynamic
Definition: Constants.h:22
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:178
ComputationInfo
Definition: Constants.h:440
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
cholmod_common * cholmodCommon() const
Definition: SuiteSparseQRSupport.h:240