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_ordering(SPQR_ORDERING_DEFAULT),
68 m_allow_tol(SPQR_DEFAULT_TOL),
69 m_tolerance (NumTraits<Scalar>::epsilon())
71 cholmod_l_start(&m_cc);
74 SPQR(
const _MatrixType& matrix)
75 : m_ordering(SPQR_ORDERING_DEFAULT),
76 m_allow_tol(SPQR_DEFAULT_TOL),
77 m_tolerance (NumTraits<Scalar>::epsilon())
79 cholmod_l_start(&m_cc);
86 cholmod_l_free_sparse(&m_H, &m_cc);
87 cholmod_l_free_sparse(&m_cR, &m_cc);
88 cholmod_l_free_dense(&m_HTau, &m_cc);
91 cholmod_l_finish(&m_cc);
93 void compute(
const _MatrixType& matrix)
95 MatrixType mat(matrix);
98 Index col = matrix.cols();
99 m_rank = SuiteSparseQR<Scalar>(m_ordering, m_tolerance, col, &A,
100 &m_cR, &m_E, &m_H, &m_HPinv, &m_HTau, &m_cc);
105 m_isInitialized =
false;
109 m_isInitialized =
true;
110 m_isRUpToDate =
false;
115 inline Index
rows()
const {
return m_H->nrow; }
120 inline Index
cols()
const {
return m_cR->ncol; }
126 template<
typename Rhs>
129 eigen_assert(m_isInitialized &&
" The QR factorization should be computed first, call compute()");
130 eigen_assert(this->
rows()==B.rows()
131 &&
"SPQR::solve(): invalid number of rows of the right hand side matrix B");
132 return internal::solve_retval<SPQR, Rhs>(*
this, B.derived());
135 template<
typename Rhs,
typename Dest>
138 eigen_assert(m_isInitialized &&
" The QR factorization should be computed first, call compute()");
139 eigen_assert(b.cols()==1 &&
"This method is for vectors only");
145 Index rk = this->
rank();
147 y.bottomRows(
cols()-rk).setZero();
158 eigen_assert(m_isInitialized &&
" The QR factorization should be computed first, call compute()");
160 m_R = viewAsEigen<Scalar,ColMajor, typename MatrixType::Index>(*m_cR);
161 m_isRUpToDate =
true;
168 return SPQRMatrixQReturnType<SPQR>(*this);
173 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
174 Index n = m_cR->ncol;
176 for(Index j = 0; j <n; j++) colsPerm.
indices()(j) = m_E[j];
186 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
187 return m_cc.SPQR_istat[4];
205 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
209 bool m_isInitialized;
211 bool m_factorizationIsOk;
212 mutable bool m_isRUpToDate;
216 RealScalar m_tolerance;
217 mutable cholmod_sparse *m_cR;
218 mutable MatrixType m_R;
220 mutable cholmod_sparse *m_H;
221 mutable Index *m_HPinv;
222 mutable cholmod_dense *m_HTau;
223 mutable Index m_rank;
224 mutable cholmod_common m_cc;
225 template<
typename ,
typename >
friend struct SPQR_QProduct;
228 template <
typename SPQRType,
typename Derived>
229 struct SPQR_QProduct : ReturnByValue<SPQR_QProduct<SPQRType,Derived> >
231 typedef typename SPQRType::Scalar Scalar;
232 typedef typename SPQRType::Index Index;
234 SPQR_QProduct(
const SPQRType& spqr,
const Derived& other,
bool transpose) : m_spqr(spqr),m_other(other),m_transpose(transpose) {}
236 inline Index rows()
const {
return m_transpose ? m_spqr.rows() : m_spqr.cols(); }
237 inline Index cols()
const {
return m_other.cols(); }
239 template<
typename ResType>
240 void evalTo(ResType& res)
const
244 int method = m_transpose ? SPQR_QTX : SPQR_QX;
245 cholmod_common *cc = m_spqr.cholmodCommon();
247 x_cd = SuiteSparseQR_qmult<Scalar>(method, m_spqr.m_H, m_spqr.m_HTau, m_spqr.m_HPinv, &y_cd, cc);
248 res = Matrix<Scalar,ResType::RowsAtCompileTime,ResType::ColsAtCompileTime>::Map(reinterpret_cast<Scalar*>(x_cd->x), x_cd->nrow, x_cd->ncol);
249 cholmod_l_free_dense(&x_cd, cc);
251 const SPQRType& m_spqr;
252 const Derived& m_other;
256 template<
typename SPQRType>
257 struct SPQRMatrixQReturnType{
259 SPQRMatrixQReturnType(
const SPQRType& spqr) : m_spqr(spqr) {}
260 template<
typename Derived>
261 SPQR_QProduct<SPQRType, Derived>
operator*(
const MatrixBase<Derived>& other)
263 return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(),
false);
265 SPQRMatrixQTransposeReturnType<SPQRType> adjoint()
const
267 return SPQRMatrixQTransposeReturnType<SPQRType>(m_spqr);
270 SPQRMatrixQTransposeReturnType<SPQRType> transpose()
const
272 return SPQRMatrixQTransposeReturnType<SPQRType>(m_spqr);
274 const SPQRType& m_spqr;
277 template<
typename SPQRType>
278 struct SPQRMatrixQTransposeReturnType{
279 SPQRMatrixQTransposeReturnType(
const SPQRType& spqr) : m_spqr(spqr) {}
280 template<
typename Derived>
281 SPQR_QProduct<SPQRType,Derived>
operator*(
const MatrixBase<Derived>& other)
283 return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(),
true);
285 const SPQRType& m_spqr;
290 template<
typename _MatrixType,
typename Rhs>
291 struct solve_retval<SPQR<_MatrixType>, Rhs>
292 : solve_retval_base<SPQR<_MatrixType>, Rhs>
294 typedef SPQR<_MatrixType> Dec;
295 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
297 template<typename Dest>
void evalTo(Dest& dst)
const
299 dec()._solve(rhs(),dst);