1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2012-2013 Desire Nuentsa <desire.nuentsa_wakam@inria.fr>
5// Copyright (C) 2012-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_SPARSE_QR_H
12#define EIGEN_SPARSE_QR_H
13
14namespace Eigen {
15
16template<typename MatrixType, typename OrderingType> class SparseQR;
17template<typename SparseQRType> struct SparseQRMatrixQReturnType;
18template<typename SparseQRType> struct SparseQRMatrixQTransposeReturnType;
19template<typename SparseQRType, typename Derived> struct SparseQR_QProduct;
20namespace internal {
21  template <typename SparseQRType> struct traits<SparseQRMatrixQReturnType<SparseQRType> >
22  {
23    typedef typename SparseQRType::MatrixType ReturnType;
24    typedef typename ReturnType::Index Index;
25    typedef typename ReturnType::StorageKind StorageKind;
26  };
27  template <typename SparseQRType> struct traits<SparseQRMatrixQTransposeReturnType<SparseQRType> >
28  {
29    typedef typename SparseQRType::MatrixType ReturnType;
30  };
31  template <typename SparseQRType, typename Derived> struct traits<SparseQR_QProduct<SparseQRType, Derived> >
32  {
33    typedef typename Derived::PlainObject ReturnType;
34  };
35} // End namespace internal
36
37/**
38  * \ingroup SparseQR_Module
39  * \class SparseQR
40  * \brief Sparse left-looking rank-revealing QR factorization
41  *
42  * This class implements a left-looking rank-revealing QR decomposition
43  * of sparse matrices. When a column has a norm less than a given tolerance
44  * it is implicitly permuted to the end. The QR factorization thus obtained is
45  * given by A*P = Q*R where R is upper triangular or trapezoidal.
46  *
47  * P is the column permutation which is the product of the fill-reducing and the
48  * rank-revealing permutations. Use colsPermutation() to get it.
49  *
50  * Q is the orthogonal matrix represented as products of Householder reflectors.
51  * Use matrixQ() to get an expression and matrixQ().transpose() to get the transpose.
52  * You can then apply it to a vector.
53  *
54  * R is the sparse triangular or trapezoidal matrix. The later occurs when A is rank-deficient.
55  * matrixR().topLeftCorner(rank(), rank()) always returns a triangular factor of full rank.
56  *
57  * \tparam _MatrixType The type of the sparse matrix A, must be a column-major SparseMatrix<>
58  * \tparam _OrderingType The fill-reducing ordering method. See the \link OrderingMethods_Module
59  *  OrderingMethods \endlink module for the list of built-in and external ordering methods.
60  *
61  * \warning The input sparse matrix A must be in compressed mode (see SparseMatrix::makeCompressed()).
62  *
63  */
64template<typename _MatrixType, typename _OrderingType>
65class SparseQR
66{
67  public:
68    typedef _MatrixType MatrixType;
69    typedef _OrderingType OrderingType;
70    typedef typename MatrixType::Scalar Scalar;
71    typedef typename MatrixType::RealScalar RealScalar;
72    typedef typename MatrixType::Index Index;
73    typedef SparseMatrix<Scalar,ColMajor,Index> QRMatrixType;
74    typedef Matrix<Index, Dynamic, 1> IndexVector;
75    typedef Matrix<Scalar, Dynamic, 1> ScalarVector;
76    typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType;
77  public:
78    SparseQR () : m_isInitialized(false), m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false)
79    { }
80
81    /** Construct a QR factorization of the matrix \a mat.
82      *
83      * \warning The matrix \a mat must be in compressed mode (see SparseMatrix::makeCompressed()).
84      *
85      * \sa compute()
86      */
87    SparseQR(const MatrixType& mat) : m_isInitialized(false), m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false)
88    {
89      compute(mat);
90    }
91
92    /** Computes the QR factorization of the sparse matrix \a mat.
93      *
94      * \warning The matrix \a mat must be in compressed mode (see SparseMatrix::makeCompressed()).
95      *
96      * \sa analyzePattern(), factorize()
97      */
98    void compute(const MatrixType& mat)
99    {
100      analyzePattern(mat);
101      factorize(mat);
102    }
103    void analyzePattern(const MatrixType& mat);
104    void factorize(const MatrixType& mat);
105
106    /** \returns the number of rows of the represented matrix.
107      */
108    inline Index rows() const { return m_pmat.rows(); }
109
110    /** \returns the number of columns of the represented matrix.
111      */
112    inline Index cols() const { return m_pmat.cols();}
113
114    /** \returns a const reference to the \b sparse upper triangular matrix R of the QR factorization.
115      */
116    const QRMatrixType& matrixR() const { return m_R; }
117
118    /** \returns the number of non linearly dependent columns as determined by the pivoting threshold.
119      *
120      * \sa setPivotThreshold()
121      */
122    Index rank() const
123    {
124      eigen_assert(m_isInitialized && "The factorization should be called first, use compute()");
125      return m_nonzeropivots;
126    }
127
128    /** \returns an expression of the matrix Q as products of sparse Householder reflectors.
129    * The common usage of this function is to apply it to a dense matrix or vector
130    * \code
131    * VectorXd B1, B2;
132    * // Initialize B1
133    * B2 = matrixQ() * B1;
134    * \endcode
135    *
136    * To get a plain SparseMatrix representation of Q:
137    * \code
138    * SparseMatrix<double> Q;
139    * Q = SparseQR<SparseMatrix<double> >(A).matrixQ();
140    * \endcode
141    * Internally, this call simply performs a sparse product between the matrix Q
142    * and a sparse identity matrix. However, due to the fact that the sparse
143    * reflectors are stored unsorted, two transpositions are needed to sort
144    * them before performing the product.
145    */
146    SparseQRMatrixQReturnType<SparseQR> matrixQ() const
147    { return SparseQRMatrixQReturnType<SparseQR>(*this); }
148
149    /** \returns a const reference to the column permutation P that was applied to A such that A*P = Q*R
150      * It is the combination of the fill-in reducing permutation and numerical column pivoting.
151      */
152    const PermutationType& colsPermutation() const
153    {
154      eigen_assert(m_isInitialized && "Decomposition is not initialized.");
155      return m_outputPerm_c;
156    }
157
158    /** \returns A string describing the type of error.
159      * This method is provided to ease debugging, not to handle errors.
160      */
161    std::string lastErrorMessage() const { return m_lastError; }
162
163    /** \internal */
164    template<typename Rhs, typename Dest>
165    bool _solve(const MatrixBase<Rhs> &B, MatrixBase<Dest> &dest) const
166    {
167      eigen_assert(m_isInitialized && "The factorization should be called first, use compute()");
168      eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix");
169
170      Index rank = this->rank();
171
172      // Compute Q^T * b;
173      typename Dest::PlainObject y, b;
174      y = this->matrixQ().transpose() * B;
175      b = y;
176
177      // Solve with the triangular matrix R
178      y.resize((std::max)(cols(),Index(y.rows())),y.cols());
179      y.topRows(rank) = this->matrixR().topLeftCorner(rank, rank).template triangularView<Upper>().solve(b.topRows(rank));
180      y.bottomRows(y.rows()-rank).setZero();
181
182      // Apply the column permutation
183      if (m_perm_c.size())  dest = colsPermutation() * y.topRows(cols());
184      else                  dest = y.topRows(cols());
185
186      m_info = Success;
187      return true;
188    }
189
190
191    /** Sets the threshold that is used to determine linearly dependent columns during the factorization.
192      *
193      * In practice, if during the factorization the norm of the column that has to be eliminated is below
194      * this threshold, then the entire column is treated as zero, and it is moved at the end.
195      */
196    void setPivotThreshold(const RealScalar& threshold)
197    {
198      m_useDefaultThreshold = false;
199      m_threshold = threshold;
200    }
201
202    /** \returns the solution X of \f$ A X = B \f$ using the current decomposition of A.
203      *
204      * \sa compute()
205      */
206    template<typename Rhs>
207    inline const internal::solve_retval<SparseQR, Rhs> solve(const MatrixBase<Rhs>& B) const
208    {
209      eigen_assert(m_isInitialized && "The factorization should be called first, use compute()");
210      eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix");
211      return internal::solve_retval<SparseQR, Rhs>(*this, B.derived());
212    }
213    template<typename Rhs>
214    inline const internal::sparse_solve_retval<SparseQR, Rhs> solve(const SparseMatrixBase<Rhs>& B) const
215    {
216          eigen_assert(m_isInitialized && "The factorization should be called first, use compute()");
217          eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix");
218          return internal::sparse_solve_retval<SparseQR, Rhs>(*this, B.derived());
219    }
220
221    /** \brief Reports whether previous computation was successful.
222      *
223      * \returns \c Success if computation was successful,
224      *          \c NumericalIssue if the QR factorization reports a numerical problem
225      *          \c InvalidInput if the input matrix is invalid
226      *
227      * \sa iparm()
228      */
229    ComputationInfo info() const
230    {
231      eigen_assert(m_isInitialized && "Decomposition is not initialized.");
232      return m_info;
233    }
234
235  protected:
236    inline void sort_matrix_Q()
237    {
238      if(this->m_isQSorted) return;
239      // The matrix Q is sorted during the transposition
240      SparseMatrix<Scalar, RowMajor, Index> mQrm(this->m_Q);
241      this->m_Q = mQrm;
242      this->m_isQSorted = true;
243    }
244
245
246  protected:
247    bool m_isInitialized;
248    bool m_analysisIsok;
249    bool m_factorizationIsok;
250    mutable ComputationInfo m_info;
251    std::string m_lastError;
252    QRMatrixType m_pmat;            // Temporary matrix
253    QRMatrixType m_R;               // The triangular factor matrix
254    QRMatrixType m_Q;               // The orthogonal reflectors
255    ScalarVector m_hcoeffs;         // The Householder coefficients
256    PermutationType m_perm_c;       // Fill-reducing  Column  permutation
257    PermutationType m_pivotperm;    // The permutation for rank revealing
258    PermutationType m_outputPerm_c; // The final column permutation
259    RealScalar m_threshold;         // Threshold to determine null Householder reflections
260    bool m_useDefaultThreshold;     // Use default threshold
261    Index m_nonzeropivots;          // Number of non zero pivots found
262    IndexVector m_etree;            // Column elimination tree
263    IndexVector m_firstRowElt;      // First element in each row
264    bool m_isQSorted;               // whether Q is sorted or not
265
266    template <typename, typename > friend struct SparseQR_QProduct;
267    template <typename > friend struct SparseQRMatrixQReturnType;
268
269};
270
271/** \brief Preprocessing step of a QR factorization
272  *
273  * \warning The matrix \a mat must be in compressed mode (see SparseMatrix::makeCompressed()).
274  *
275  * In this step, the fill-reducing permutation is computed and applied to the columns of A
276  * and the column elimination tree is computed as well. Only the sparsity pattern of \a mat is exploited.
277  *
278  * \note In this step it is assumed that there is no empty row in the matrix \a mat.
279  */
280template <typename MatrixType, typename OrderingType>
281void SparseQR<MatrixType,OrderingType>::analyzePattern(const MatrixType& mat)
282{
283  eigen_assert(mat.isCompressed() && "SparseQR requires a sparse matrix in compressed mode. Call .makeCompressed() before passing it to SparseQR");
284  // Compute the column fill reducing ordering
285  OrderingType ord;
286  ord(mat, m_perm_c);
287  Index n = mat.cols();
288  Index m = mat.rows();
289  Index diagSize = (std::min)(m,n);
290
291  if (!m_perm_c.size())
292  {
293    m_perm_c.resize(n);
294    m_perm_c.indices().setLinSpaced(n, 0,n-1);
295  }
296
297  // Compute the column elimination tree of the permuted matrix
298  m_outputPerm_c = m_perm_c.inverse();
299  internal::coletree(mat, m_etree, m_firstRowElt, m_outputPerm_c.indices().data());
300
301  m_R.resize(m, n);
302  m_Q.resize(m, diagSize);
303
304  // Allocate space for nonzero elements : rough estimation
305  m_R.reserve(2*mat.nonZeros()); //FIXME Get a more accurate estimation through symbolic factorization with the etree
306  m_Q.reserve(2*mat.nonZeros());
307  m_hcoeffs.resize(diagSize);
308  m_analysisIsok = true;
309}
310
311/** \brief Performs the numerical QR factorization of the input matrix
312  *
313  * The function SparseQR::analyzePattern(const MatrixType&) must have been called beforehand with
314  * a matrix having the same sparsity pattern than \a mat.
315  *
316  * \param mat The sparse column-major matrix
317  */
318template <typename MatrixType, typename OrderingType>
319void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat)
320{
321  using std::abs;
322  using std::max;
323
324  eigen_assert(m_analysisIsok && "analyzePattern() should be called before this step");
325  Index m = mat.rows();
326  Index n = mat.cols();
327  Index diagSize = (std::min)(m,n);
328  IndexVector mark((std::max)(m,n)); mark.setConstant(-1);  // Record the visited nodes
329  IndexVector Ridx(n), Qidx(m);                             // Store temporarily the row indexes for the current column of R and Q
330  Index nzcolR, nzcolQ;                                     // Number of nonzero for the current column of R and Q
331  ScalarVector tval(m);                                     // The dense vector used to compute the current column
332  RealScalar pivotThreshold = m_threshold;
333
334  m_pmat = mat;
335  m_pmat.uncompress(); // To have the innerNonZeroPtr allocated
336  // Apply the fill-in reducing permutation lazily:
337  for (int i = 0; i < n; i++)
338  {
339    Index p = m_perm_c.size() ? m_perm_c.indices()(i) : i;
340    m_pmat.outerIndexPtr()[p] = mat.outerIndexPtr()[i];
341    m_pmat.innerNonZeroPtr()[p] = mat.outerIndexPtr()[i+1] - mat.outerIndexPtr()[i];
342  }
343
344  /* Compute the default threshold as in MatLab, see:
345   * Tim Davis, "Algorithm 915, SuiteSparseQR: Multifrontal Multithreaded Rank-Revealing
346   * Sparse QR Factorization, ACM Trans. on Math. Soft. 38(1), 2011, Page 8:3
347   */
348  if(m_useDefaultThreshold)
349  {
350    RealScalar max2Norm = 0.0;
351    for (int j = 0; j < n; j++) max2Norm = (max)(max2Norm, m_pmat.col(j).norm());
352    pivotThreshold = 20 * (m + n) * max2Norm * NumTraits<RealScalar>::epsilon();
353  }
354
355  // Initialize the numerical permutation
356  m_pivotperm.setIdentity(n);
357
358  Index nonzeroCol = 0; // Record the number of valid pivots
359  m_Q.startVec(0);
360
361  // Left looking rank-revealing QR factorization: compute a column of R and Q at a time
362  for (Index col = 0; col < n; ++col)
363  {
364    mark.setConstant(-1);
365    m_R.startVec(col);
366    mark(nonzeroCol) = col;
367    Qidx(0) = nonzeroCol;
368    nzcolR = 0; nzcolQ = 1;
369    bool found_diag = nonzeroCol>=m;
370    tval.setZero();
371
372    // Symbolic factorization: find the nonzero locations of the column k of the factors R and Q, i.e.,
373    // all the nodes (with indexes lower than rank) reachable through the column elimination tree (etree) rooted at node k.
374    // Note: if the diagonal entry does not exist, then its contribution must be explicitly added,
375    // thus the trick with found_diag that permits to do one more iteration on the diagonal element if this one has not been found.
376    for (typename MatrixType::InnerIterator itp(m_pmat, col); itp || !found_diag; ++itp)
377    {
378      Index curIdx = nonzeroCol;
379      if(itp) curIdx = itp.row();
380      if(curIdx == nonzeroCol) found_diag = true;
381
382      // Get the nonzeros indexes of the current column of R
383      Index st = m_firstRowElt(curIdx); // The traversal of the etree starts here
384      if (st < 0 )
385      {
386        m_lastError = "Empty row found during numerical factorization";
387        m_info = InvalidInput;
388        return;
389      }
390
391      // Traverse the etree
392      Index bi = nzcolR;
393      for (; mark(st) != col; st = m_etree(st))
394      {
395        Ridx(nzcolR) = st;  // Add this row to the list,
396        mark(st) = col;     // and mark this row as visited
397        nzcolR++;
398      }
399
400      // Reverse the list to get the topological ordering
401      Index nt = nzcolR-bi;
402      for(Index i = 0; i < nt/2; i++) std::swap(Ridx(bi+i), Ridx(nzcolR-i-1));
403
404      // Copy the current (curIdx,pcol) value of the input matrix
405      if(itp) tval(curIdx) = itp.value();
406      else    tval(curIdx) = Scalar(0);
407
408      // Compute the pattern of Q(:,k)
409      if(curIdx > nonzeroCol && mark(curIdx) != col )
410      {
411        Qidx(nzcolQ) = curIdx;  // Add this row to the pattern of Q,
412        mark(curIdx) = col;     // and mark it as visited
413        nzcolQ++;
414      }
415    }
416
417    // Browse all the indexes of R(:,col) in reverse order
418    for (Index i = nzcolR-1; i >= 0; i--)
419    {
420      Index curIdx = Ridx(i);
421
422      // Apply the curIdx-th householder vector to the current column (temporarily stored into tval)
423      Scalar tdot(0);
424
425      // First compute q' * tval
426      tdot = m_Q.col(curIdx).dot(tval);
427
428      tdot *= m_hcoeffs(curIdx);
429
430      // Then update tval = tval - q * tau
431      // FIXME: tval -= tdot * m_Q.col(curIdx) should amount to the same (need to check/add support for efficient "dense ?= sparse")
432      for (typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq)
433        tval(itq.row()) -= itq.value() * tdot;
434
435      // Detect fill-in for the current column of Q
436      if(m_etree(Ridx(i)) == nonzeroCol)
437      {
438        for (typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq)
439        {
440          Index iQ = itq.row();
441          if (mark(iQ) != col)
442          {
443            Qidx(nzcolQ++) = iQ;  // Add this row to the pattern of Q,
444            mark(iQ) = col;       // and mark it as visited
445          }
446        }
447      }
448    } // End update current column
449
450    Scalar tau;
451    RealScalar beta = 0;
452
453    if(nonzeroCol < diagSize)
454    {
455      // Compute the Householder reflection that eliminate the current column
456      // FIXME this step should call the Householder module.
457      Scalar c0 = nzcolQ ? tval(Qidx(0)) : Scalar(0);
458
459      // First, the squared norm of Q((col+1):m, col)
460      RealScalar sqrNorm = 0.;
461      for (Index itq = 1; itq < nzcolQ; ++itq) sqrNorm += numext::abs2(tval(Qidx(itq)));
462      if(sqrNorm == RealScalar(0) && numext::imag(c0) == RealScalar(0))
463      {
464        tau = RealScalar(0);
465        beta = numext::real(c0);
466        tval(Qidx(0)) = 1;
467      }
468      else
469      {
470        using std::sqrt;
471        beta = sqrt(numext::abs2(c0) + sqrNorm);
472        if(numext::real(c0) >= RealScalar(0))
473          beta = -beta;
474        tval(Qidx(0)) = 1;
475        for (Index itq = 1; itq < nzcolQ; ++itq)
476          tval(Qidx(itq)) /= (c0 - beta);
477        tau = numext::conj((beta-c0) / beta);
478
479      }
480    }
481
482    // Insert values in R
483    for (Index  i = nzcolR-1; i >= 0; i--)
484    {
485      Index curIdx = Ridx(i);
486      if(curIdx < nonzeroCol)
487      {
488        m_R.insertBackByOuterInnerUnordered(col, curIdx) = tval(curIdx);
489        tval(curIdx) = Scalar(0.);
490      }
491    }
492
493    if(nonzeroCol < diagSize && abs(beta) >= pivotThreshold)
494    {
495      m_R.insertBackByOuterInner(col, nonzeroCol) = beta;
496      // The householder coefficient
497      m_hcoeffs(nonzeroCol) = tau;
498      // Record the householder reflections
499      for (Index itq = 0; itq < nzcolQ; ++itq)
500      {
501        Index iQ = Qidx(itq);
502        m_Q.insertBackByOuterInnerUnordered(nonzeroCol,iQ) = tval(iQ);
503        tval(iQ) = Scalar(0.);
504      }
505      nonzeroCol++;
506      if(nonzeroCol<diagSize)
507        m_Q.startVec(nonzeroCol);
508    }
509    else
510    {
511      // Zero pivot found: move implicitly this column to the end
512      for (Index j = nonzeroCol; j < n-1; j++)
513        std::swap(m_pivotperm.indices()(j), m_pivotperm.indices()[j+1]);
514
515      // Recompute the column elimination tree
516      internal::coletree(m_pmat, m_etree, m_firstRowElt, m_pivotperm.indices().data());
517    }
518  }
519
520  m_hcoeffs.tail(diagSize-nonzeroCol).setZero();
521
522  // Finalize the column pointers of the sparse matrices R and Q
523  m_Q.finalize();
524  m_Q.makeCompressed();
525  m_R.finalize();
526  m_R.makeCompressed();
527  m_isQSorted = false;
528
529  m_nonzeropivots = nonzeroCol;
530
531  if(nonzeroCol<n)
532  {
533    // Permute the triangular factor to put the 'dead' columns to the end
534    MatrixType tempR(m_R);
535    m_R = tempR * m_pivotperm;
536
537    // Update the column permutation
538    m_outputPerm_c = m_outputPerm_c * m_pivotperm;
539  }
540
541  m_isInitialized = true;
542  m_factorizationIsok = true;
543  m_info = Success;
544}
545
546namespace internal {
547
548template<typename _MatrixType, typename OrderingType, typename Rhs>
549struct solve_retval<SparseQR<_MatrixType,OrderingType>, Rhs>
550  : solve_retval_base<SparseQR<_MatrixType,OrderingType>, Rhs>
551{
552  typedef SparseQR<_MatrixType,OrderingType> Dec;
553  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
554
555  template<typename Dest> void evalTo(Dest& dst) const
556  {
557    dec()._solve(rhs(),dst);
558  }
559};
560template<typename _MatrixType, typename OrderingType, typename Rhs>
561struct sparse_solve_retval<SparseQR<_MatrixType, OrderingType>, Rhs>
562 : sparse_solve_retval_base<SparseQR<_MatrixType, OrderingType>, Rhs>
563{
564  typedef SparseQR<_MatrixType, OrderingType> Dec;
565  EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec, Rhs)
566
567  template<typename Dest> void evalTo(Dest& dst) const
568  {
569    this->defaultEvalTo(dst);
570  }
571};
572} // end namespace internal
573
574template <typename SparseQRType, typename Derived>
575struct SparseQR_QProduct : ReturnByValue<SparseQR_QProduct<SparseQRType, Derived> >
576{
577  typedef typename SparseQRType::QRMatrixType MatrixType;
578  typedef typename SparseQRType::Scalar Scalar;
579  typedef typename SparseQRType::Index Index;
580  // Get the references
581  SparseQR_QProduct(const SparseQRType& qr, const Derived& other, bool transpose) :
582  m_qr(qr),m_other(other),m_transpose(transpose) {}
583  inline Index rows() const { return m_transpose ? m_qr.rows() : m_qr.cols(); }
584  inline Index cols() const { return m_other.cols(); }
585
586  // Assign to a vector
587  template<typename DesType>
588  void evalTo(DesType& res) const
589  {
590    Index m = m_qr.rows();
591    Index n = m_qr.cols();
592    Index diagSize = (std::min)(m,n);
593    res = m_other;
594    if (m_transpose)
595    {
596      eigen_assert(m_qr.m_Q.rows() == m_other.rows() && "Non conforming object sizes");
597      //Compute res = Q' * other column by column
598      for(Index j = 0; j < res.cols(); j++){
599        for (Index k = 0; k < diagSize; k++)
600        {
601          Scalar tau = Scalar(0);
602          tau = m_qr.m_Q.col(k).dot(res.col(j));
603          if(tau==Scalar(0)) continue;
604          tau = tau * m_qr.m_hcoeffs(k);
605          res.col(j) -= tau * m_qr.m_Q.col(k);
606        }
607      }
608    }
609    else
610    {
611      eigen_assert(m_qr.m_Q.rows() == m_other.rows() && "Non conforming object sizes");
612      // Compute res = Q * other column by column
613      for(Index j = 0; j < res.cols(); j++)
614      {
615        for (Index k = diagSize-1; k >=0; k--)
616        {
617          Scalar tau = Scalar(0);
618          tau = m_qr.m_Q.col(k).dot(res.col(j));
619          if(tau==Scalar(0)) continue;
620          tau = tau * m_qr.m_hcoeffs(k);
621          res.col(j) -= tau * m_qr.m_Q.col(k);
622        }
623      }
624    }
625  }
626
627  const SparseQRType& m_qr;
628  const Derived& m_other;
629  bool m_transpose;
630};
631
632template<typename SparseQRType>
633struct SparseQRMatrixQReturnType : public EigenBase<SparseQRMatrixQReturnType<SparseQRType> >
634{
635  typedef typename SparseQRType::Index Index;
636  typedef typename SparseQRType::Scalar Scalar;
637  typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
638  SparseQRMatrixQReturnType(const SparseQRType& qr) : m_qr(qr) {}
639  template<typename Derived>
640  SparseQR_QProduct<SparseQRType, Derived> operator*(const MatrixBase<Derived>& other)
641  {
642    return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(),false);
643  }
644  SparseQRMatrixQTransposeReturnType<SparseQRType> adjoint() const
645  {
646    return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
647  }
648  inline Index rows() const { return m_qr.rows(); }
649  inline Index cols() const { return (std::min)(m_qr.rows(),m_qr.cols()); }
650  // To use for operations with the transpose of Q
651  SparseQRMatrixQTransposeReturnType<SparseQRType> transpose() const
652  {
653    return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
654  }
655  template<typename Dest> void evalTo(MatrixBase<Dest>& dest) const
656  {
657    dest.derived() = m_qr.matrixQ() * Dest::Identity(m_qr.rows(), m_qr.rows());
658  }
659  template<typename Dest> void evalTo(SparseMatrixBase<Dest>& dest) const
660  {
661    Dest idMat(m_qr.rows(), m_qr.rows());
662    idMat.setIdentity();
663    // Sort the sparse householder reflectors if needed
664    const_cast<SparseQRType *>(&m_qr)->sort_matrix_Q();
665    dest.derived() = SparseQR_QProduct<SparseQRType, Dest>(m_qr, idMat, false);
666  }
667
668  const SparseQRType& m_qr;
669};
670
671template<typename SparseQRType>
672struct SparseQRMatrixQTransposeReturnType
673{
674  SparseQRMatrixQTransposeReturnType(const SparseQRType& qr) : m_qr(qr) {}
675  template<typename Derived>
676  SparseQR_QProduct<SparseQRType,Derived> operator*(const MatrixBase<Derived>& other)
677  {
678    return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(), true);
679  }
680  const SparseQRType& m_qr;
681};
682
683} // end namespace Eigen
684
685#endif
686