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