1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2012 Desire Nuentsa <desire.nuentsa_wakam@inria.fr>
5 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 
10 #ifndef EIGEN_SUITESPARSEQRSUPPORT_H
11 #define EIGEN_SUITESPARSEQRSUPPORT_H
12 
13 namespace Eigen {
14 
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;
19   namespace internal {
20     template <typename SPQRType> struct traits<SPQRMatrixQReturnType<SPQRType> >
21     {
22       typedef typename SPQRType::MatrixType ReturnType;
23     };
24     template <typename SPQRType> struct traits<SPQRMatrixQTransposeReturnType<SPQRType> >
25     {
26       typedef typename SPQRType::MatrixType ReturnType;
27     };
28     template <typename SPQRType, typename Derived> struct traits<SPQR_QProduct<SPQRType, Derived> >
29     {
30       typedef typename Derived::PlainObject ReturnType;
31     };
32   } // End namespace internal
33 
34 /**
35  * \ingroup SPQRSupport_Module
36  * \class SPQR
37  * \brief Sparse QR factorization based on SuiteSparseQR library
38  *
39  * This class is used to perform a multithreaded and multifrontal rank-revealing QR decomposition
40  * of sparse matrices. The result is then used to solve linear leasts_square systems.
41  * Clearly, a QR factorization is returned such that A*P = Q*R where :
42  *
43  * P is the column permutation. Use colsPermutation() to get it.
44  *
45  * Q is the orthogonal matrix represented as Householder reflectors.
46  * Use matrixQ() to get an expression and matrixQ().transpose() to get the transpose.
47  * You can then apply it to a vector.
48  *
49  * R is the sparse triangular factor. Use matrixQR() to get it as SparseMatrix.
50  * NOTE : The Index type of R is always UF_long. You can get it with SPQR::Index
51  *
52  * \tparam _MatrixType The type of the sparse matrix A, must be a column-major SparseMatrix<>
53  * NOTE
54  *
55  */
56 template<typename _MatrixType>
57 class SPQR
58 {
59   public:
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;
65   public:
66     SPQR()
67       : m_isInitialized(false),
68       m_ordering(SPQR_ORDERING_DEFAULT),
69       m_allow_tol(SPQR_DEFAULT_TOL),
70       m_tolerance (NumTraits<Scalar>::epsilon())
71     {
72       cholmod_l_start(&m_cc);
73     }
74 
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())
80     {
81       cholmod_l_start(&m_cc);
82       compute(matrix);
83     }
84 
85     ~SPQR()
86     {
87       SPQR_free();
88       cholmod_l_finish(&m_cc);
89     }
90     void SPQR_free()
91     {
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);
95       std::free(m_E);
96       std::free(m_HPinv);
97     }
98 
99     void compute(const _MatrixType& matrix)
100     {
101       if(m_isInitialized) SPQR_free();
102 
103       MatrixType mat(matrix);
104       cholmod_sparse A;
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);
109 
110       if (!m_cR)
111       {
112         m_info = NumericalIssue;
113         m_isInitialized = false;
114         return;
115       }
116       m_info = Success;
117       m_isInitialized = true;
118       m_isRUpToDate = false;
119     }
120     /**
121      * Get the number of rows of the input matrix and the Q matrix
122      */
123     inline Index rows() const {return m_H->nrow; }
124 
125     /**
126      * Get the number of columns of the input matrix.
127      */
128     inline Index cols() const { return m_cR->ncol; }
129 
130       /** \returns the solution X of \f$ A X = B \f$ using the current decomposition of A.
131       *
132       * \sa compute()
133       */
134     template<typename Rhs>
135     inline const internal::solve_retval<SPQR, Rhs> solve(const MatrixBase<Rhs>& B) const
136     {
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());
141     }
142 
143     template<typename Rhs, typename Dest>
144     void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
145     {
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");
148 
149       //Compute Q^T * b
150       typename Dest::PlainObject y;
151       y = matrixQ().transpose() * b;
152         // Solves with the triangular matrix R
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();
156       // Apply the column permutation
157       dest.topRows(cols()) = colsPermutation() * y.topRows(cols());
158 
159       m_info = Success;
160     }
161 
162     /** \returns the sparse triangular factor R. It is a sparse matrix
163      */
164     const MatrixType matrixR() const
165     {
166       eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()");
167       if(!m_isRUpToDate) {
168         m_R = viewAsEigen<Scalar,ColMajor, typename MatrixType::Index>(*m_cR);
169         m_isRUpToDate = true;
170       }
171       return m_R;
172     }
173     /// Get an expression of the matrix Q
174     SPQRMatrixQReturnType<SPQR> matrixQ() const
175     {
176       return SPQRMatrixQReturnType<SPQR>(*this);
177     }
178     /// Get the permutation that was applied to columns of A
179     PermutationType colsPermutation() const
180     {
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];
185       return colsPerm;
186 
187     }
188     /**
189      * Gets the rank of the matrix.
190      * It should be equal to matrixQR().cols if the matrix is full-rank
191      */
192     Index rank() const
193     {
194       eigen_assert(m_isInitialized && "Decomposition is not initialized.");
195       return m_cc.SPQR_istat[4];
196     }
197     /// Set the fill-reducing ordering method to be used
198     void setSPQROrdering(int ord) { m_ordering = ord;}
199     /// Set the tolerance tol to treat columns with 2-norm < =tol as zero
200     void setPivotThreshold(const RealScalar& tol) { m_tolerance = tol; }
201 
202     /** \returns a pointer to the SPQR workspace */
203     cholmod_common *cholmodCommon() const { return &m_cc; }
204 
205 
206     /** \brief Reports whether previous computation was successful.
207       *
208       * \returns \c Success if computation was succesful,
209       *          \c NumericalIssue if the sparse QR can not be computed
210       */
211     ComputationInfo info() const
212     {
213       eigen_assert(m_isInitialized && "Decomposition is not initialized.");
214       return m_info;
215     }
216   protected:
217     bool m_isInitialized;
218     bool m_analysisIsOk;
219     bool m_factorizationIsOk;
220     mutable bool m_isRUpToDate;
221     mutable ComputationInfo m_info;
222     int m_ordering; // Ordering method to use, see SPQR's manual
223     int m_allow_tol; // Allow to use some tolerance during numerical factorization.
224     RealScalar m_tolerance; // treat columns with 2-norm below this tolerance as zero
225     mutable cholmod_sparse *m_cR; // The sparse R factor in cholmod format
226     mutable MatrixType m_R; // The sparse matrix R in Eigen format
227     mutable Index *m_E; // The permutation applied to columns
228     mutable cholmod_sparse *m_H;  //The householder vectors
229     mutable Index *m_HPinv; // The row permutation of H
230     mutable cholmod_dense *m_HTau; // The Householder coefficients
231     mutable Index m_rank; // The rank of the matrix
232     mutable cholmod_common m_cc; // Workspace and parameters
233     template<typename ,typename > friend struct SPQR_QProduct;
234 };
235 
236 template <typename SPQRType, typename Derived>
237 struct SPQR_QProduct : ReturnByValue<SPQR_QProduct<SPQRType,Derived> >
238 {
239   typedef typename SPQRType::Scalar Scalar;
240   typedef typename SPQRType::Index Index;
241   //Define the constructor to get reference to argument types
242   SPQR_QProduct(const SPQRType& spqr, const Derived& other, bool transpose) : m_spqr(spqr),m_other(other),m_transpose(transpose) {}
243 
244   inline Index rows() const { return m_transpose ? m_spqr.rows() : m_spqr.cols(); }
245   inline Index cols() const { return m_other.cols(); }
246   // Assign to a vector
247   template<typename ResType>
248   void evalTo(ResType& res) const
249   {
250     cholmod_dense y_cd;
251     cholmod_dense *x_cd;
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);
256     res = Matrix<Scalar,ResType::RowsAtCompileTime,ResType::ColsAtCompileTime>::Map(reinterpret_cast<Scalar*>(x_cd->x), x_cd->nrow, x_cd->ncol);
257     cholmod_l_free_dense(&x_cd, cc);
258   }
259   const SPQRType& m_spqr;
260   const Derived& m_other;
261   bool m_transpose;
262 
263 };
264 template<typename SPQRType>
265 struct SPQRMatrixQReturnType{
266 
267   SPQRMatrixQReturnType(const SPQRType& spqr) : m_spqr(spqr) {}
268   template<typename Derived>
269   SPQR_QProduct<SPQRType, Derived> operator*(const MatrixBase<Derived>& other)
270   {
271     return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(),false);
272   }
273   SPQRMatrixQTransposeReturnType<SPQRType> adjoint() const
274   {
275     return SPQRMatrixQTransposeReturnType<SPQRType>(m_spqr);
276   }
277   // To use for operations with the transpose of Q
278   SPQRMatrixQTransposeReturnType<SPQRType> transpose() const
279   {
280     return SPQRMatrixQTransposeReturnType<SPQRType>(m_spqr);
281   }
282   const SPQRType& m_spqr;
283 };
284 
285 template<typename SPQRType>
286 struct SPQRMatrixQTransposeReturnType{
287   SPQRMatrixQTransposeReturnType(const SPQRType& spqr) : m_spqr(spqr) {}
288   template<typename Derived>
289   SPQR_QProduct<SPQRType,Derived> operator*(const MatrixBase<Derived>& other)
290   {
291     return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(), true);
292   }
293   const SPQRType& m_spqr;
294 };
295 
296 namespace internal {
297 
298 template<typename _MatrixType, typename Rhs>
299 struct solve_retval<SPQR<_MatrixType>, Rhs>
300   : solve_retval_base<SPQR<_MatrixType>, Rhs>
301 {
302   typedef SPQR<_MatrixType> Dec;
303   EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
304 
305   template<typename Dest> void evalTo(Dest& dst) const
306   {
307     dec()._solve(rhs(),dst);
308   }
309 };
310 
311 } // end namespace internal
312 
313 }// End namespace Eigen
314 #endif
315