Home | History | Annotate | Download | only in SparseQR
      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 (at) inria.fr>
      5 // Copyright (C) 2012-2014 Gael Guennebaud <gael.guennebaud (at) 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 
     14 namespace Eigen {
     15 
     16 template<typename MatrixType, typename OrderingType> class SparseQR;
     17 template<typename SparseQRType> struct SparseQRMatrixQReturnType;
     18 template<typename SparseQRType> struct SparseQRMatrixQTransposeReturnType;
     19 template<typename SparseQRType, typename Derived> struct SparseQR_QProduct;
     20 namespace 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   */
     64 template<typename _MatrixType, typename _OrderingType>
     65 class 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   */
    280 template <typename MatrixType, typename OrderingType>
    281 void 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   */
    318 template <typename MatrixType, typename OrderingType>
    319 void 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 
    546 namespace internal {
    547 
    548 template<typename _MatrixType, typename OrderingType, typename Rhs>
    549 struct 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 };
    560 template<typename _MatrixType, typename OrderingType, typename Rhs>
    561 struct 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 
    574 template <typename SparseQRType, typename Derived>
    575 struct 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 
    632 template<typename SparseQRType>
    633 struct 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 
    671 template<typename SparseQRType>
    672 struct 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