Eigen  3.3.3
SparseQR.h
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2012-2013 Desire Nuentsa <desire.nuentsa_wakam@inria.fr>
00005 // Copyright (C) 2012-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
00006 //
00007 // This Source Code Form is subject to the terms of the Mozilla
00008 // Public License v. 2.0. If a copy of the MPL was not distributed
00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00010 
00011 #ifndef EIGEN_SPARSE_QR_H
00012 #define EIGEN_SPARSE_QR_H
00013 
00014 namespace Eigen {
00015 
00016 template<typename MatrixType, typename OrderingType> class SparseQR;
00017 template<typename SparseQRType> struct SparseQRMatrixQReturnType;
00018 template<typename SparseQRType> struct SparseQRMatrixQTransposeReturnType;
00019 template<typename SparseQRType, typename Derived> struct SparseQR_QProduct;
00020 namespace internal {
00021   template <typename SparseQRType> struct traits<SparseQRMatrixQReturnType<SparseQRType> >
00022   {
00023     typedef typename SparseQRType::MatrixType ReturnType;
00024     typedef typename ReturnType::StorageIndex StorageIndex;
00025     typedef typename ReturnType::StorageKind StorageKind;
00026     enum {
00027       RowsAtCompileTime = Dynamic,
00028       ColsAtCompileTime = Dynamic
00029     };
00030   };
00031   template <typename SparseQRType> struct traits<SparseQRMatrixQTransposeReturnType<SparseQRType> >
00032   {
00033     typedef typename SparseQRType::MatrixType ReturnType;
00034   };
00035   template <typename SparseQRType, typename Derived> struct traits<SparseQR_QProduct<SparseQRType, Derived> >
00036   {
00037     typedef typename Derived::PlainObject ReturnType;
00038   };
00039 } // End namespace internal
00040 
00070 template<typename _MatrixType, typename _OrderingType>
00071 class SparseQR : public SparseSolverBase<SparseQR<_MatrixType,_OrderingType> >
00072 {
00073   protected:
00074     typedef SparseSolverBase<SparseQR<_MatrixType,_OrderingType> > Base;
00075     using Base::m_isInitialized;
00076   public:
00077     using Base::_solve_impl;
00078     typedef _MatrixType MatrixType;
00079     typedef _OrderingType OrderingType;
00080     typedef typename MatrixType::Scalar Scalar;
00081     typedef typename MatrixType::RealScalar RealScalar;
00082     typedef typename MatrixType::StorageIndex StorageIndex;
00083     typedef SparseMatrix<Scalar,ColMajor,StorageIndex> QRMatrixType;
00084     typedef Matrix<StorageIndex, Dynamic, 1> IndexVector;
00085     typedef Matrix<Scalar, Dynamic, 1> ScalarVector;
00086     typedef PermutationMatrix<Dynamic, Dynamic, StorageIndex> PermutationType;
00087 
00088     enum {
00089       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00090       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00091     };
00092     
00093   public:
00094     SparseQR () :  m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false)
00095     { }
00096     
00103     explicit SparseQR(const MatrixType& mat) : m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false)
00104     {
00105       compute(mat);
00106     }
00107     
00114     void compute(const MatrixType& mat)
00115     {
00116       analyzePattern(mat);
00117       factorize(mat);
00118     }
00119     void analyzePattern(const MatrixType& mat);
00120     void factorize(const MatrixType& mat);
00121     
00124     inline Index rows() const { return m_pmat.rows(); }
00125     
00128     inline Index cols() const { return m_pmat.cols();}
00129     
00143     const QRMatrixType& matrixR() const { return m_R; }
00144     
00149     Index rank() const
00150     {
00151       eigen_assert(m_isInitialized && "The factorization should be called first, use compute()");
00152       return m_nonzeropivots; 
00153     }
00154     
00173     SparseQRMatrixQReturnType<SparseQR> matrixQ() const 
00174     { return SparseQRMatrixQReturnType<SparseQR>(*this); }
00175     
00179     const PermutationType& colsPermutation() const
00180     { 
00181       eigen_assert(m_isInitialized && "Decomposition is not initialized.");
00182       return m_outputPerm_c;
00183     }
00184     
00188     std::string lastErrorMessage() const { return m_lastError; }
00189     
00191     template<typename Rhs, typename Dest>
00192     bool _solve_impl(const MatrixBase<Rhs> &B, MatrixBase<Dest> &dest) const
00193     {
00194       eigen_assert(m_isInitialized && "The factorization should be called first, use compute()");
00195       eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix");
00196 
00197       Index rank = this->rank();
00198       
00199       // Compute Q^T * b;
00200       typename Dest::PlainObject y, b;
00201       y = this->matrixQ().transpose() * B; 
00202       b = y;
00203       
00204       // Solve with the triangular matrix R
00205       y.resize((std::max<Index>)(cols(),y.rows()),y.cols());
00206       y.topRows(rank) = this->matrixR().topLeftCorner(rank, rank).template triangularView<Upper>().solve(b.topRows(rank));
00207       y.bottomRows(y.rows()-rank).setZero();
00208       
00209       // Apply the column permutation
00210       if (m_perm_c.size())  dest = colsPermutation() * y.topRows(cols());
00211       else                  dest = y.topRows(cols());
00212       
00213       m_info = Success;
00214       return true;
00215     }
00216 
00222     void setPivotThreshold(const RealScalar& threshold)
00223     {
00224       m_useDefaultThreshold = false;
00225       m_threshold = threshold;
00226     }
00227     
00232     template<typename Rhs>
00233     inline const Solve<SparseQR, Rhs> solve(const MatrixBase<Rhs>& B) const 
00234     {
00235       eigen_assert(m_isInitialized && "The factorization should be called first, use compute()");
00236       eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix");
00237       return Solve<SparseQR, Rhs>(*this, B.derived());
00238     }
00239     template<typename Rhs>
00240     inline const Solve<SparseQR, Rhs> solve(const SparseMatrixBase<Rhs>& B) const
00241     {
00242           eigen_assert(m_isInitialized && "The factorization should be called first, use compute()");
00243           eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix");
00244           return Solve<SparseQR, Rhs>(*this, B.derived());
00245     }
00246     
00255     ComputationInfo info() const
00256     {
00257       eigen_assert(m_isInitialized && "Decomposition is not initialized.");
00258       return m_info;
00259     }
00260 
00261 
00263     inline void _sort_matrix_Q()
00264     {
00265       if(this->m_isQSorted) return;
00266       // The matrix Q is sorted during the transposition
00267       SparseMatrix<Scalar, RowMajor, Index> mQrm(this->m_Q);
00268       this->m_Q = mQrm;
00269       this->m_isQSorted = true;
00270     }
00271 
00272     
00273   protected:
00274     bool m_analysisIsok;
00275     bool m_factorizationIsok;
00276     mutable ComputationInfo m_info;
00277     std::string m_lastError;
00278     QRMatrixType m_pmat;            // Temporary matrix
00279     QRMatrixType m_R;               // The triangular factor matrix
00280     QRMatrixType m_Q;               // The orthogonal reflectors
00281     ScalarVector m_hcoeffs;         // The Householder coefficients
00282     PermutationType m_perm_c;       // Fill-reducing  Column  permutation
00283     PermutationType m_pivotperm;    // The permutation for rank revealing
00284     PermutationType m_outputPerm_c; // The final column permutation
00285     RealScalar m_threshold;         // Threshold to determine null Householder reflections
00286     bool m_useDefaultThreshold;     // Use default threshold
00287     Index m_nonzeropivots;          // Number of non zero pivots found
00288     IndexVector m_etree;            // Column elimination tree
00289     IndexVector m_firstRowElt;      // First element in each row
00290     bool m_isQSorted;               // whether Q is sorted or not
00291     bool m_isEtreeOk;               // whether the elimination tree match the initial input matrix
00292     
00293     template <typename, typename > friend struct SparseQR_QProduct;
00294     
00295 };
00296 
00306 template <typename MatrixType, typename OrderingType>
00307 void SparseQR<MatrixType,OrderingType>::analyzePattern(const MatrixType& mat)
00308 {
00309   eigen_assert(mat.isCompressed() && "SparseQR requires a sparse matrix in compressed mode. Call .makeCompressed() before passing it to SparseQR");
00310   // Copy to a column major matrix if the input is rowmajor
00311   typename internal::conditional<MatrixType::IsRowMajor,QRMatrixType,const MatrixType&>::type matCpy(mat);
00312   // Compute the column fill reducing ordering
00313   OrderingType ord; 
00314   ord(matCpy, m_perm_c); 
00315   Index n = mat.cols();
00316   Index m = mat.rows();
00317   Index diagSize = (std::min)(m,n);
00318   
00319   if (!m_perm_c.size())
00320   {
00321     m_perm_c.resize(n);
00322     m_perm_c.indices().setLinSpaced(n, 0,StorageIndex(n-1));
00323   }
00324   
00325   // Compute the column elimination tree of the permuted matrix
00326   m_outputPerm_c = m_perm_c.inverse();
00327   internal::coletree(matCpy, m_etree, m_firstRowElt, m_outputPerm_c.indices().data());
00328   m_isEtreeOk = true;
00329   
00330   m_R.resize(m, n);
00331   m_Q.resize(m, diagSize);
00332   
00333   // Allocate space for nonzero elements : rough estimation
00334   m_R.reserve(2*mat.nonZeros()); //FIXME Get a more accurate estimation through symbolic factorization with the etree
00335   m_Q.reserve(2*mat.nonZeros());
00336   m_hcoeffs.resize(diagSize);
00337   m_analysisIsok = true;
00338 }
00339 
00347 template <typename MatrixType, typename OrderingType>
00348 void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat)
00349 {
00350   using std::abs;
00351   
00352   eigen_assert(m_analysisIsok && "analyzePattern() should be called before this step");
00353   StorageIndex m = StorageIndex(mat.rows());
00354   StorageIndex n = StorageIndex(mat.cols());
00355   StorageIndex diagSize = (std::min)(m,n);
00356   IndexVector mark((std::max)(m,n)); mark.setConstant(-1);  // Record the visited nodes
00357   IndexVector Ridx(n), Qidx(m);                             // Store temporarily the row indexes for the current column of R and Q
00358   Index nzcolR, nzcolQ;                                     // Number of nonzero for the current column of R and Q
00359   ScalarVector tval(m);                                     // The dense vector used to compute the current column
00360   RealScalar pivotThreshold = m_threshold;
00361   
00362   m_R.setZero();
00363   m_Q.setZero();
00364   m_pmat = mat;
00365   if(!m_isEtreeOk)
00366   {
00367     m_outputPerm_c = m_perm_c.inverse();
00368     internal::coletree(m_pmat, m_etree, m_firstRowElt, m_outputPerm_c.indices().data());
00369     m_isEtreeOk = true;
00370   }
00371 
00372   m_pmat.uncompress(); // To have the innerNonZeroPtr allocated
00373   
00374   // Apply the fill-in reducing permutation lazily:
00375   {
00376     // If the input is row major, copy the original column indices,
00377     // otherwise directly use the input matrix
00378     // 
00379     IndexVector originalOuterIndicesCpy;
00380     const StorageIndex *originalOuterIndices = mat.outerIndexPtr();
00381     if(MatrixType::IsRowMajor)
00382     {
00383       originalOuterIndicesCpy = IndexVector::Map(m_pmat.outerIndexPtr(),n+1);
00384       originalOuterIndices = originalOuterIndicesCpy.data();
00385     }
00386     
00387     for (int i = 0; i < n; i++)
00388     {
00389       Index p = m_perm_c.size() ? m_perm_c.indices()(i) : i;
00390       m_pmat.outerIndexPtr()[p] = originalOuterIndices[i]; 
00391       m_pmat.innerNonZeroPtr()[p] = originalOuterIndices[i+1] - originalOuterIndices[i]; 
00392     }
00393   }
00394   
00395   /* Compute the default threshold as in MatLab, see:
00396    * Tim Davis, "Algorithm 915, SuiteSparseQR: Multifrontal Multithreaded Rank-Revealing
00397    * Sparse QR Factorization, ACM Trans. on Math. Soft. 38(1), 2011, Page 8:3 
00398    */
00399   if(m_useDefaultThreshold) 
00400   {
00401     RealScalar max2Norm = 0.0;
00402     for (int j = 0; j < n; j++) max2Norm = numext::maxi(max2Norm, m_pmat.col(j).norm());
00403     if(max2Norm==RealScalar(0))
00404       max2Norm = RealScalar(1);
00405     pivotThreshold = 20 * (m + n) * max2Norm * NumTraits<RealScalar>::epsilon();
00406   }
00407   
00408   // Initialize the numerical permutation
00409   m_pivotperm.setIdentity(n);
00410   
00411   StorageIndex nonzeroCol = 0; // Record the number of valid pivots
00412   m_Q.startVec(0);
00413 
00414   // Left looking rank-revealing QR factorization: compute a column of R and Q at a time
00415   for (StorageIndex col = 0; col < n; ++col)
00416   {
00417     mark.setConstant(-1);
00418     m_R.startVec(col);
00419     mark(nonzeroCol) = col;
00420     Qidx(0) = nonzeroCol;
00421     nzcolR = 0; nzcolQ = 1;
00422     bool found_diag = nonzeroCol>=m;
00423     tval.setZero(); 
00424     
00425     // Symbolic factorization: find the nonzero locations of the column k of the factors R and Q, i.e.,
00426     // all the nodes (with indexes lower than rank) reachable through the column elimination tree (etree) rooted at node k.
00427     // Note: if the diagonal entry does not exist, then its contribution must be explicitly added,
00428     // thus the trick with found_diag that permits to do one more iteration on the diagonal element if this one has not been found.
00429     for (typename QRMatrixType::InnerIterator itp(m_pmat, col); itp || !found_diag; ++itp)
00430     {
00431       StorageIndex curIdx = nonzeroCol;
00432       if(itp) curIdx = StorageIndex(itp.row());
00433       if(curIdx == nonzeroCol) found_diag = true;
00434       
00435       // Get the nonzeros indexes of the current column of R
00436       StorageIndex st = m_firstRowElt(curIdx); // The traversal of the etree starts here
00437       if (st < 0 )
00438       {
00439         m_lastError = "Empty row found during numerical factorization";
00440         m_info = InvalidInput;
00441         return;
00442       }
00443 
00444       // Traverse the etree 
00445       Index bi = nzcolR;
00446       for (; mark(st) != col; st = m_etree(st))
00447       {
00448         Ridx(nzcolR) = st;  // Add this row to the list,
00449         mark(st) = col;     // and mark this row as visited
00450         nzcolR++;
00451       }
00452 
00453       // Reverse the list to get the topological ordering
00454       Index nt = nzcolR-bi;
00455       for(Index i = 0; i < nt/2; i++) std::swap(Ridx(bi+i), Ridx(nzcolR-i-1));
00456        
00457       // Copy the current (curIdx,pcol) value of the input matrix
00458       if(itp) tval(curIdx) = itp.value();
00459       else    tval(curIdx) = Scalar(0);
00460       
00461       // Compute the pattern of Q(:,k)
00462       if(curIdx > nonzeroCol && mark(curIdx) != col ) 
00463       {
00464         Qidx(nzcolQ) = curIdx;  // Add this row to the pattern of Q,
00465         mark(curIdx) = col;     // and mark it as visited
00466         nzcolQ++;
00467       }
00468     }
00469 
00470     // Browse all the indexes of R(:,col) in reverse order
00471     for (Index i = nzcolR-1; i >= 0; i--)
00472     {
00473       Index curIdx = Ridx(i);
00474       
00475       // Apply the curIdx-th householder vector to the current column (temporarily stored into tval)
00476       Scalar tdot(0);
00477       
00478       // First compute q' * tval
00479       tdot = m_Q.col(curIdx).dot(tval);
00480 
00481       tdot *= m_hcoeffs(curIdx);
00482       
00483       // Then update tval = tval - q * tau
00484       // FIXME: tval -= tdot * m_Q.col(curIdx) should amount to the same (need to check/add support for efficient "dense ?= sparse")
00485       for (typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq)
00486         tval(itq.row()) -= itq.value() * tdot;
00487 
00488       // Detect fill-in for the current column of Q
00489       if(m_etree(Ridx(i)) == nonzeroCol)
00490       {
00491         for (typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq)
00492         {
00493           StorageIndex iQ = StorageIndex(itq.row());
00494           if (mark(iQ) != col)
00495           {
00496             Qidx(nzcolQ++) = iQ;  // Add this row to the pattern of Q,
00497             mark(iQ) = col;       // and mark it as visited
00498           }
00499         }
00500       }
00501     } // End update current column
00502     
00503     Scalar tau = RealScalar(0);
00504     RealScalar beta = 0;
00505     
00506     if(nonzeroCol < diagSize)
00507     {
00508       // Compute the Householder reflection that eliminate the current column
00509       // FIXME this step should call the Householder module.
00510       Scalar c0 = nzcolQ ? tval(Qidx(0)) : Scalar(0);
00511       
00512       // First, the squared norm of Q((col+1):m, col)
00513       RealScalar sqrNorm = 0.;
00514       for (Index itq = 1; itq < nzcolQ; ++itq) sqrNorm += numext::abs2(tval(Qidx(itq)));
00515       if(sqrNorm == RealScalar(0) && numext::imag(c0) == RealScalar(0))
00516       {
00517         beta = numext::real(c0);
00518         tval(Qidx(0)) = 1;
00519       }
00520       else
00521       {
00522         using std::sqrt;
00523         beta = sqrt(numext::abs2(c0) + sqrNorm);
00524         if(numext::real(c0) >= RealScalar(0))
00525           beta = -beta;
00526         tval(Qidx(0)) = 1;
00527         for (Index itq = 1; itq < nzcolQ; ++itq)
00528           tval(Qidx(itq)) /= (c0 - beta);
00529         tau = numext::conj((beta-c0) / beta);
00530           
00531       }
00532     }
00533 
00534     // Insert values in R
00535     for (Index  i = nzcolR-1; i >= 0; i--)
00536     {
00537       Index curIdx = Ridx(i);
00538       if(curIdx < nonzeroCol) 
00539       {
00540         m_R.insertBackByOuterInnerUnordered(col, curIdx) = tval(curIdx);
00541         tval(curIdx) = Scalar(0.);
00542       }
00543     }
00544 
00545     if(nonzeroCol < diagSize && abs(beta) >= pivotThreshold)
00546     {
00547       m_R.insertBackByOuterInner(col, nonzeroCol) = beta;
00548       // The householder coefficient
00549       m_hcoeffs(nonzeroCol) = tau;
00550       // Record the householder reflections
00551       for (Index itq = 0; itq < nzcolQ; ++itq)
00552       {
00553         Index iQ = Qidx(itq);
00554         m_Q.insertBackByOuterInnerUnordered(nonzeroCol,iQ) = tval(iQ);
00555         tval(iQ) = Scalar(0.);
00556       }
00557       nonzeroCol++;
00558       if(nonzeroCol<diagSize)
00559         m_Q.startVec(nonzeroCol);
00560     }
00561     else
00562     {
00563       // Zero pivot found: move implicitly this column to the end
00564       for (Index j = nonzeroCol; j < n-1; j++) 
00565         std::swap(m_pivotperm.indices()(j), m_pivotperm.indices()[j+1]);
00566       
00567       // Recompute the column elimination tree
00568       internal::coletree(m_pmat, m_etree, m_firstRowElt, m_pivotperm.indices().data());
00569       m_isEtreeOk = false;
00570     }
00571   }
00572   
00573   m_hcoeffs.tail(diagSize-nonzeroCol).setZero();
00574   
00575   // Finalize the column pointers of the sparse matrices R and Q
00576   m_Q.finalize();
00577   m_Q.makeCompressed();
00578   m_R.finalize();
00579   m_R.makeCompressed();
00580   m_isQSorted = false;
00581 
00582   m_nonzeropivots = nonzeroCol;
00583   
00584   if(nonzeroCol<n)
00585   {
00586     // Permute the triangular factor to put the 'dead' columns to the end
00587     QRMatrixType tempR(m_R);
00588     m_R = tempR * m_pivotperm;
00589     
00590     // Update the column permutation
00591     m_outputPerm_c = m_outputPerm_c * m_pivotperm;
00592   }
00593   
00594   m_isInitialized = true; 
00595   m_factorizationIsok = true;
00596   m_info = Success;
00597 }
00598 
00599 template <typename SparseQRType, typename Derived>
00600 struct SparseQR_QProduct : ReturnByValue<SparseQR_QProduct<SparseQRType, Derived> >
00601 {
00602   typedef typename SparseQRType::QRMatrixType MatrixType;
00603   typedef typename SparseQRType::Scalar Scalar;
00604   // Get the references 
00605   SparseQR_QProduct(const SparseQRType& qr, const Derived& other, bool transpose) : 
00606   m_qr(qr),m_other(other),m_transpose(transpose) {}
00607   inline Index rows() const { return m_transpose ? m_qr.rows() : m_qr.cols(); }
00608   inline Index cols() const { return m_other.cols(); }
00609   
00610   // Assign to a vector
00611   template<typename DesType>
00612   void evalTo(DesType& res) const
00613   {
00614     Index m = m_qr.rows();
00615     Index n = m_qr.cols();
00616     Index diagSize = (std::min)(m,n);
00617     res = m_other;
00618     if (m_transpose)
00619     {
00620       eigen_assert(m_qr.m_Q.rows() == m_other.rows() && "Non conforming object sizes");
00621       //Compute res = Q' * other column by column
00622       for(Index j = 0; j < res.cols(); j++){
00623         for (Index k = 0; k < diagSize; k++)
00624         {
00625           Scalar tau = Scalar(0);
00626           tau = m_qr.m_Q.col(k).dot(res.col(j));
00627           if(tau==Scalar(0)) continue;
00628           tau = tau * m_qr.m_hcoeffs(k);
00629           res.col(j) -= tau * m_qr.m_Q.col(k);
00630         }
00631       }
00632     }
00633     else
00634     {
00635       eigen_assert(m_qr.m_Q.rows() == m_other.rows() && "Non conforming object sizes");
00636       // Compute res = Q * other column by column
00637       for(Index j = 0; j < res.cols(); j++)
00638       {
00639         for (Index k = diagSize-1; k >=0; k--)
00640         {
00641           Scalar tau = Scalar(0);
00642           tau = m_qr.m_Q.col(k).dot(res.col(j));
00643           if(tau==Scalar(0)) continue;
00644           tau = tau * m_qr.m_hcoeffs(k);
00645           res.col(j) -= tau * m_qr.m_Q.col(k);
00646         }
00647       }
00648     }
00649   }
00650   
00651   const SparseQRType& m_qr;
00652   const Derived& m_other;
00653   bool m_transpose;
00654 };
00655 
00656 template<typename SparseQRType>
00657 struct SparseQRMatrixQReturnType : public EigenBase<SparseQRMatrixQReturnType<SparseQRType> >
00658 {  
00659   typedef typename SparseQRType::Scalar Scalar;
00660   typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
00661   enum {
00662     RowsAtCompileTime = Dynamic,
00663     ColsAtCompileTime = Dynamic
00664   };
00665   explicit SparseQRMatrixQReturnType(const SparseQRType& qr) : m_qr(qr) {}
00666   template<typename Derived>
00667   SparseQR_QProduct<SparseQRType, Derived> operator*(const MatrixBase<Derived>& other)
00668   {
00669     return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(),false);
00670   }
00671   SparseQRMatrixQTransposeReturnType<SparseQRType> adjoint() const
00672   {
00673     return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
00674   }
00675   inline Index rows() const { return m_qr.rows(); }
00676   inline Index cols() const { return (std::min)(m_qr.rows(),m_qr.cols()); }
00677   // To use for operations with the transpose of Q
00678   SparseQRMatrixQTransposeReturnType<SparseQRType> transpose() const
00679   {
00680     return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
00681   }
00682   const SparseQRType& m_qr;
00683 };
00684 
00685 template<typename SparseQRType>
00686 struct SparseQRMatrixQTransposeReturnType
00687 {
00688   explicit SparseQRMatrixQTransposeReturnType(const SparseQRType& qr) : m_qr(qr) {}
00689   template<typename Derived>
00690   SparseQR_QProduct<SparseQRType,Derived> operator*(const MatrixBase<Derived>& other)
00691   {
00692     return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(), true);
00693   }
00694   const SparseQRType& m_qr;
00695 };
00696 
00697 namespace internal {
00698   
00699 template<typename SparseQRType>
00700 struct evaluator_traits<SparseQRMatrixQReturnType<SparseQRType> >
00701 {
00702   typedef typename SparseQRType::MatrixType MatrixType;
00703   typedef typename storage_kind_to_evaluator_kind<typename MatrixType::StorageKind>::Kind Kind;
00704   typedef SparseShape Shape;
00705 };
00706 
00707 template< typename DstXprType, typename SparseQRType>
00708 struct Assignment<DstXprType, SparseQRMatrixQReturnType<SparseQRType>, internal::assign_op<typename DstXprType::Scalar,typename DstXprType::Scalar>, Sparse2Sparse>
00709 {
00710   typedef SparseQRMatrixQReturnType<SparseQRType> SrcXprType;
00711   typedef typename DstXprType::Scalar Scalar;
00712   typedef typename DstXprType::StorageIndex StorageIndex;
00713   static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,Scalar> &/*func*/)
00714   {
00715     typename DstXprType::PlainObject idMat(src.m_qr.rows(), src.m_qr.rows());
00716     idMat.setIdentity();
00717     // Sort the sparse householder reflectors if needed
00718     const_cast<SparseQRType *>(&src.m_qr)->_sort_matrix_Q();
00719     dst = SparseQR_QProduct<SparseQRType, DstXprType>(src.m_qr, idMat, false);
00720   }
00721 };
00722 
00723 template< typename DstXprType, typename SparseQRType>
00724 struct Assignment<DstXprType, SparseQRMatrixQReturnType<SparseQRType>, internal::assign_op<typename DstXprType::Scalar,typename DstXprType::Scalar>, Sparse2Dense>
00725 {
00726   typedef SparseQRMatrixQReturnType<SparseQRType> SrcXprType;
00727   typedef typename DstXprType::Scalar Scalar;
00728   typedef typename DstXprType::StorageIndex StorageIndex;
00729   static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,Scalar> &/*func*/)
00730   {
00731     dst = src.m_qr.matrixQ() * DstXprType::Identity(src.m_qr.rows(), src.m_qr.rows());
00732   }
00733 };
00734 
00735 } // end namespace internal
00736 
00737 } // end namespace Eigen
00738 
00739 #endif
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends