EIC code displayed by LXR master/​include/​Rivet/​Math/​eigen3/​src/​SparseQR/​SparseQR.h	Source navigation Diff markup Identifier search General search
	Version: master test Revert   Change

File indexing completed on 2025-04-19 09:06:44

 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra.
 //
 // Copyright (C) 2012-2013 Desire Nuentsa <desire.nuentsa_wakam@inria.fr>
 // Copyright (C) 2012-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
 //
 // This Source Code Form is subject to the terms of the Mozilla
 // Public License v. 2.0. If a copy of the MPL was not distributed
 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
 #ifndef EIGEN_SPARSE_QR_H
 #define EIGEN_SPARSE_QR_H
 
 namespace RivetEigen {
 
 template<typename MatrixType, typename OrderingType> class SparseQR;
 template<typename SparseQRType> struct SparseQRMatrixQReturnType;
 template<typename SparseQRType> struct SparseQRMatrixQTransposeReturnType;
 template<typename SparseQRType, typename Derived> struct SparseQR_QProduct;
 namespace internal {
   template <typename SparseQRType> struct traits<SparseQRMatrixQReturnType<SparseQRType> >
   {
     typedef typename SparseQRType::MatrixType ReturnType;
     typedef typename ReturnType::StorageIndex StorageIndex;
     typedef typename ReturnType::StorageKind StorageKind;
     enum {
       RowsAtCompileTime = Dynamic,
       ColsAtCompileTime = Dynamic
     };
   };
   template <typename SparseQRType> struct traits<SparseQRMatrixQTransposeReturnType<SparseQRType> >
   {
     typedef typename SparseQRType::MatrixType ReturnType;
   };
   template <typename SparseQRType, typename Derived> struct traits<SparseQR_QProduct<SparseQRType, Derived> >
   {
     typedef typename Derived::PlainObject ReturnType;
   };
 } // End namespace internal
 
 /**
   * \ingroup SparseQR_Module
   * \class SparseQR
   * \brief Sparse left-looking QR factorization with numerical column pivoting
   * 
   * This class implements a left-looking QR decomposition of sparse matrices
   * with numerical column pivoting.
   * When a column has a norm less than a given tolerance
   * it is implicitly permuted to the end. The QR factorization thus obtained is 
   * given by A*P = Q*R where R is upper triangular or trapezoidal. 
   * 
   * P is the column permutation which is the product of the fill-reducing and the
   * numerical permutations. Use colsPermutation() to get it.
   * 
   * Q is the orthogonal matrix represented as products of Householder reflectors. 
   * Use matrixQ() to get an expression and matrixQ().adjoint() to get the adjoint.
   * You can then apply it to a vector.
   * 
   * R is the sparse triangular or trapezoidal matrix. The later occurs when A is rank-deficient.
   * matrixR().topLeftCorner(rank(), rank()) always returns a triangular factor of full rank.
   * 
   * \tparam _MatrixType The type of the sparse matrix A, must be a column-major SparseMatrix<>
   * \tparam _OrderingType The fill-reducing ordering method. See the \link OrderingMethods_Module 
   *  OrderingMethods \endlink module for the list of built-in and external ordering methods.
   * 
   * \implsparsesolverconcept
   *
   * The numerical pivoting strategy and default threshold are the same as in SuiteSparse QR, and
   * detailed in the following paper:
   * <i>
   * Tim Davis, "Algorithm 915, SuiteSparseQR: Multifrontal Multithreaded Rank-Revealing
   * Sparse QR Factorization, ACM Trans. on Math. Soft. 38(1), 2011.
   * </i>
   * Even though it is qualified as "rank-revealing", this strategy might fail for some 
   * rank deficient problems. When this class is used to solve linear or least-square problems
   * it is thus strongly recommended to check the accuracy of the computed solution. If it
   * failed, it usually helps to increase the threshold with setPivotThreshold.
   * 
   * \warning The input sparse matrix A must be in compressed mode (see SparseMatrix::makeCompressed()).
   * \warning For complex matrices matrixQ().transpose() will actually return the adjoint matrix.
   * 
   */
 template<typename _MatrixType, typename _OrderingType>
 class SparseQR : public SparseSolverBase<SparseQR<_MatrixType,_OrderingType> >
 {
   protected:
     typedef SparseSolverBase<SparseQR<_MatrixType,_OrderingType> > Base;
     using Base::m_isInitialized;
   public:
     using Base::_solve_impl;
     typedef _MatrixType MatrixType;
     typedef _OrderingType OrderingType;
     typedef typename MatrixType::Scalar Scalar;
     typedef typename MatrixType::RealScalar RealScalar;
     typedef typename MatrixType::StorageIndex StorageIndex;
     typedef SparseMatrix<Scalar,ColMajor,StorageIndex> QRMatrixType;
     typedef Matrix<StorageIndex, Dynamic, 1> IndexVector;
     typedef Matrix<Scalar, Dynamic, 1> ScalarVector;
     typedef PermutationMatrix<Dynamic, Dynamic, StorageIndex> PermutationType;
 
     enum {
       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
     };
     
   public:
     SparseQR () :  m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false)
     { }
     
     /** Construct a QR factorization of the matrix \a mat.
       * 
       * \warning The matrix \a mat must be in compressed mode (see SparseMatrix::makeCompressed()).
       * 
       * \sa compute()
       */
     explicit SparseQR(const MatrixType& mat) : m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false)
     {
       compute(mat);
     }
     
     /** Computes the QR factorization of the sparse matrix \a mat.
       * 
       * \warning The matrix \a mat must be in compressed mode (see SparseMatrix::makeCompressed()).
       * 
       * \sa analyzePattern(), factorize()
       */
     void compute(const MatrixType& mat)
     {
       analyzePattern(mat);
       factorize(mat);
     }
     void analyzePattern(const MatrixType& mat);
     void factorize(const MatrixType& mat);
     
     /** \returns the number of rows of the represented matrix. 
       */
     inline Index rows() const { return m_pmat.rows(); }
     
     /** \returns the number of columns of the represented matrix. 
       */
     inline Index cols() const { return m_pmat.cols();}
     
     /** \returns a const reference to the \b sparse upper triangular matrix R of the QR factorization.
       * \warning The entries of the returned matrix are not sorted. This means that using it in algorithms
       *          expecting sorted entries will fail. This include random coefficient accesses (SpaseMatrix::coeff()),
       *          and coefficient-wise operations. Matrix products and triangular solves are fine though.
       *
       * To sort the entries, you can assign it to a row-major matrix, and if a column-major matrix
       * is required, you can copy it again:
       * \code
       * SparseMatrix<double>          R  = qr.matrixR();  // column-major, not sorted!
       * SparseMatrix<double,RowMajor> Rr = qr.matrixR();  // row-major, sorted
       * SparseMatrix<double>          Rc = Rr;            // column-major, sorted
       * \endcode
       */
     const QRMatrixType& matrixR() const { return m_R; }
     
     /** \returns the number of non linearly dependent columns as determined by the pivoting threshold.
       *
       * \sa setPivotThreshold()
       */
     Index rank() const
     {
       eigen_assert(m_isInitialized && "The factorization should be called first, use compute()");
       return m_nonzeropivots; 
     }
     
     /** \returns an expression of the matrix Q as products of sparse Householder reflectors.
     * The common usage of this function is to apply it to a dense matrix or vector
     * \code
     * VectorXd B1, B2;
     * // Initialize B1
     * B2 = matrixQ() * B1;
     * \endcode
     *
     * To get a plain SparseMatrix representation of Q:
     * \code
     * SparseMatrix<double> Q;
     * Q = SparseQR<SparseMatrix<double> >(A).matrixQ();
     * \endcode
     * Internally, this call simply performs a sparse product between the matrix Q
     * and a sparse identity matrix. However, due to the fact that the sparse
     * reflectors are stored unsorted, two transpositions are needed to sort
     * them before performing the product.
     */
     SparseQRMatrixQReturnType<SparseQR> matrixQ() const 
     { return SparseQRMatrixQReturnType<SparseQR>(*this); }
     
     /** \returns a const reference to the column permutation P that was applied to A such that A*P = Q*R
       * It is the combination of the fill-in reducing permutation and numerical column pivoting.
       */
     const PermutationType& colsPermutation() const
     { 
       eigen_assert(m_isInitialized && "Decomposition is not initialized.");
       return m_outputPerm_c;
     }
     
     /** \returns A string describing the type of error.
       * This method is provided to ease debugging, not to handle errors.
       */
     std::string lastErrorMessage() const { return m_lastError; }
     
     /** \internal */
     template<typename Rhs, typename Dest>
     bool _solve_impl(const MatrixBase<Rhs> &B, MatrixBase<Dest> &dest) const
     {
       eigen_assert(m_isInitialized && "The factorization should be called first, use compute()");
       eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix");
 
       Index rank = this->rank();
       
       // Compute Q^* * b;
       typename Dest::PlainObject y, b;
       y = this->matrixQ().adjoint() * B;
       b = y;
       
       // Solve with the triangular matrix R
       y.resize((std::max<Index>)(cols(),y.rows()),y.cols());
       y.topRows(rank) = this->matrixR().topLeftCorner(rank, rank).template triangularView<Upper>().solve(b.topRows(rank));
       y.bottomRows(y.rows()-rank).setZero();
       
       // Apply the column permutation
       if (m_perm_c.size())  dest = colsPermutation() * y.topRows(cols());
       else                  dest = y.topRows(cols());
       
       m_info = Success;
       return true;
     }
 
     /** Sets the threshold that is used to determine linearly dependent columns during the factorization.
       *
       * In practice, if during the factorization the norm of the column that has to be eliminated is below
       * this threshold, then the entire column is treated as zero, and it is moved at the end.
       */
     void setPivotThreshold(const RealScalar& threshold)
     {
       m_useDefaultThreshold = false;
       m_threshold = threshold;
     }
     
     /** \returns the solution X of \f$ A X = B \f$ using the current decomposition of A.
       *
       * \sa compute()
       */
     template<typename Rhs>
     inline const Solve<SparseQR, Rhs> solve(const MatrixBase<Rhs>& B) const 
     {
       eigen_assert(m_isInitialized && "The factorization should be called first, use compute()");
       eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix");
       return Solve<SparseQR, Rhs>(*this, B.derived());
     }
     template<typename Rhs>
     inline const Solve<SparseQR, Rhs> solve(const SparseMatrixBase<Rhs>& B) const
     {
           eigen_assert(m_isInitialized && "The factorization should be called first, use compute()");
           eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix");
           return Solve<SparseQR, Rhs>(*this, B.derived());
     }
     
     /** \brief Reports whether previous computation was successful.
       *
       * \returns \c Success if computation was successful,
       *          \c NumericalIssue if the QR factorization reports a numerical problem
       *          \c InvalidInput if the input matrix is invalid
       *
       * \sa iparm()          
       */
     ComputationInfo info() const
     {
       eigen_assert(m_isInitialized && "Decomposition is not initialized.");
       return m_info;
     }
 
 
     /** \internal */
     inline void _sort_matrix_Q()
     {
       if(this->m_isQSorted) return;
       // The matrix Q is sorted during the transposition
       SparseMatrix<Scalar, RowMajor, Index> mQrm(this->m_Q);
       this->m_Q = mQrm;
       this->m_isQSorted = true;
     }
 
     
   protected:
     bool m_analysisIsok;
     bool m_factorizationIsok;
     mutable ComputationInfo m_info;
     std::string m_lastError;
     QRMatrixType m_pmat;            // Temporary matrix
     QRMatrixType m_R;               // The triangular factor matrix
     QRMatrixType m_Q;               // The orthogonal reflectors
     ScalarVector m_hcoeffs;         // The Householder coefficients
     PermutationType m_perm_c;       // Fill-reducing  Column  permutation
     PermutationType m_pivotperm;    // The permutation for rank revealing
     PermutationType m_outputPerm_c; // The final column permutation
     RealScalar m_threshold;         // Threshold to determine null Householder reflections
     bool m_useDefaultThreshold;     // Use default threshold
     Index m_nonzeropivots;          // Number of non zero pivots found
     IndexVector m_etree;            // Column elimination tree
     IndexVector m_firstRowElt;      // First element in each row
     bool m_isQSorted;               // whether Q is sorted or not
     bool m_isEtreeOk;               // whether the elimination tree match the initial input matrix
     
     template <typename, typename > friend struct SparseQR_QProduct;
     
 };
 
 /** \brief Preprocessing step of a QR factorization 
   * 
   * \warning The matrix \a mat must be in compressed mode (see SparseMatrix::makeCompressed()).
   * 
   * In this step, the fill-reducing permutation is computed and applied to the columns of A
   * and the column elimination tree is computed as well. Only the sparsity pattern of \a mat is exploited.
   * 
   * \note In this step it is assumed that there is no empty row in the matrix \a mat.
   */
 template <typename MatrixType, typename OrderingType>
 void SparseQR<MatrixType,OrderingType>::analyzePattern(const MatrixType& mat)
 {
   eigen_assert(mat.isCompressed() && "SparseQR requires a sparse matrix in compressed mode. Call .makeCompressed() before passing it to SparseQR");
   // Copy to a column major matrix if the input is rowmajor
   typename internal::conditional<MatrixType::IsRowMajor,QRMatrixType,const MatrixType&>::type matCpy(mat);
   // Compute the column fill reducing ordering
   OrderingType ord; 
   ord(matCpy, m_perm_c); 
   Index n = mat.cols();
   Index m = mat.rows();
   Index diagSize = (std::min)(m,n);
   
   if (!m_perm_c.size())
   {
     m_perm_c.resize(n);
     m_perm_c.indices().setLinSpaced(n, 0,StorageIndex(n-1));
   }
   
   // Compute the column elimination tree of the permuted matrix
   m_outputPerm_c = m_perm_c.inverse();
   internal::coletree(matCpy, m_etree, m_firstRowElt, m_outputPerm_c.indices().data());
   m_isEtreeOk = true;
   
   m_R.resize(m, n);
   m_Q.resize(m, diagSize);
   
   // Allocate space for nonzero elements: rough estimation
   m_R.reserve(2*mat.nonZeros()); //FIXME Get a more accurate estimation through symbolic factorization with the etree
   m_Q.reserve(2*mat.nonZeros());
   m_hcoeffs.resize(diagSize);
   m_analysisIsok = true;
 }
 
 /** \brief Performs the numerical QR factorization of the input matrix
   * 
   * The function SparseQR::analyzePattern(const MatrixType&) must have been called beforehand with
   * a matrix having the same sparsity pattern than \a mat.
   * 
   * \param mat The sparse column-major matrix
   */
 template <typename MatrixType, typename OrderingType>
 void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat)
 {
   using std::abs;
   
   eigen_assert(m_analysisIsok && "analyzePattern() should be called before this step");
   StorageIndex m = StorageIndex(mat.rows());
   StorageIndex n = StorageIndex(mat.cols());
   StorageIndex diagSize = (std::min)(m,n);
   IndexVector mark((std::max)(m,n)); mark.setConstant(-1);  // Record the visited nodes
   IndexVector Ridx(n), Qidx(m);                             // Store temporarily the row indexes for the current column of R and Q
   Index nzcolR, nzcolQ;                                     // Number of nonzero for the current column of R and Q
   ScalarVector tval(m);                                     // The dense vector used to compute the current column
   RealScalar pivotThreshold = m_threshold;
   
   m_R.setZero();
   m_Q.setZero();
   m_pmat = mat;
   if(!m_isEtreeOk)
   {
     m_outputPerm_c = m_perm_c.inverse();
     internal::coletree(m_pmat, m_etree, m_firstRowElt, m_outputPerm_c.indices().data());
     m_isEtreeOk = true;
   }
 
   m_pmat.uncompress(); // To have the innerNonZeroPtr allocated
   
   // Apply the fill-in reducing permutation lazily:
   {
     // If the input is row major, copy the original column indices,
     // otherwise directly use the input matrix
     // 
     IndexVector originalOuterIndicesCpy;
     const StorageIndex *originalOuterIndices = mat.outerIndexPtr();
     if(MatrixType::IsRowMajor)
     {
       originalOuterIndicesCpy = IndexVector::Map(m_pmat.outerIndexPtr(),n+1);
       originalOuterIndices = originalOuterIndicesCpy.data();
     }
     
     for (int i = 0; i < n; i++)
     {
       Index p = m_perm_c.size() ? m_perm_c.indices()(i) : i;
       m_pmat.outerIndexPtr()[p] = originalOuterIndices[i]; 
       m_pmat.innerNonZeroPtr()[p] = originalOuterIndices[i+1] - originalOuterIndices[i]; 
     }
   }
   
   /* Compute the default threshold as in MatLab, see:
    * Tim Davis, "Algorithm 915, SuiteSparseQR: Multifrontal Multithreaded Rank-Revealing
    * Sparse QR Factorization, ACM Trans. on Math. Soft. 38(1), 2011, Page 8:3 
    */
   if(m_useDefaultThreshold) 
   {
     RealScalar max2Norm = 0.0;
     for (int j = 0; j < n; j++) max2Norm = numext::maxi(max2Norm, m_pmat.col(j).norm());
     if(max2Norm==RealScalar(0))
       max2Norm = RealScalar(1);
     pivotThreshold = 20 * (m + n) * max2Norm * NumTraits<RealScalar>::epsilon();
   }
   
   // Initialize the numerical permutation
   m_pivotperm.setIdentity(n);
   
   StorageIndex nonzeroCol = 0; // Record the number of valid pivots
   m_Q.startVec(0);
 
   // Left looking rank-revealing QR factorization: compute a column of R and Q at a time
   for (StorageIndex col = 0; col < n; ++col)
   {
     mark.setConstant(-1);
     m_R.startVec(col);
     mark(nonzeroCol) = col;
     Qidx(0) = nonzeroCol;
     nzcolR = 0; nzcolQ = 1;
     bool found_diag = nonzeroCol>=m;
     tval.setZero(); 
     
     // Symbolic factorization: find the nonzero locations of the column k of the factors R and Q, i.e.,
     // all the nodes (with indexes lower than rank) reachable through the column elimination tree (etree) rooted at node k.
     // Note: if the diagonal entry does not exist, then its contribution must be explicitly added,
     // thus the trick with found_diag that permits to do one more iteration on the diagonal element if this one has not been found.
     for (typename QRMatrixType::InnerIterator itp(m_pmat, col); itp || !found_diag; ++itp)
     {
       StorageIndex curIdx = nonzeroCol;
       if(itp) curIdx = StorageIndex(itp.row());
       if(curIdx == nonzeroCol) found_diag = true;
       
       // Get the nonzeros indexes of the current column of R
       StorageIndex st = m_firstRowElt(curIdx); // The traversal of the etree starts here
       if (st < 0 )
       {
         m_lastError = "Empty row found during numerical factorization";
         m_info = InvalidInput;
         return;
       }
 
       // Traverse the etree 
       Index bi = nzcolR;
       for (; mark(st) != col; st = m_etree(st))
       {
         Ridx(nzcolR) = st;  // Add this row to the list,
         mark(st) = col;     // and mark this row as visited
         nzcolR++;
       }
 
       // Reverse the list to get the topological ordering
       Index nt = nzcolR-bi;
       for(Index i = 0; i < nt/2; i++) std::swap(Ridx(bi+i), Ridx(nzcolR-i-1));
        
       // Copy the current (curIdx,pcol) value of the input matrix
       if(itp) tval(curIdx) = itp.value();
       else    tval(curIdx) = Scalar(0);
       
       // Compute the pattern of Q(:,k)
       if(curIdx > nonzeroCol && mark(curIdx) != col ) 
       {
         Qidx(nzcolQ) = curIdx;  // Add this row to the pattern of Q,
         mark(curIdx) = col;     // and mark it as visited
         nzcolQ++;
       }
     }
 
     // Browse all the indexes of R(:,col) in reverse order
     for (Index i = nzcolR-1; i >= 0; i--)
     {
       Index curIdx = Ridx(i);
       
       // Apply the curIdx-th householder vector to the current column (temporarily stored into tval)
       Scalar tdot(0);
       
       // First compute q' * tval
       tdot = m_Q.col(curIdx).dot(tval);
 
       tdot *= m_hcoeffs(curIdx);
       
       // Then update tval = tval - q * tau
       // FIXME: tval -= tdot * m_Q.col(curIdx) should amount to the same (need to check/add support for efficient "dense ?= sparse")
       for (typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq)
         tval(itq.row()) -= itq.value() * tdot;
 
       // Detect fill-in for the current column of Q
       if(m_etree(Ridx(i)) == nonzeroCol)
       {
         for (typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq)
         {
           StorageIndex iQ = StorageIndex(itq.row());
           if (mark(iQ) != col)
           {
             Qidx(nzcolQ++) = iQ;  // Add this row to the pattern of Q,
             mark(iQ) = col;       // and mark it as visited
           }
         }
       }
     } // End update current column
     
     Scalar tau = RealScalar(0);
     RealScalar beta = 0;
     
     if(nonzeroCol < diagSize)
     {
       // Compute the Householder reflection that eliminate the current column
       // FIXME this step should call the Householder module.
       Scalar c0 = nzcolQ ? tval(Qidx(0)) : Scalar(0);
       
       // First, the squared norm of Q((col+1):m, col)
       RealScalar sqrNorm = 0.;
       for (Index itq = 1; itq < nzcolQ; ++itq) sqrNorm += numext::abs2(tval(Qidx(itq)));
       if(sqrNorm == RealScalar(0) && numext::imag(c0) == RealScalar(0))
       {
         beta = numext::real(c0);
         tval(Qidx(0)) = 1;
       }
       else
       {
         using std::sqrt;
         beta = sqrt(numext::abs2(c0) + sqrNorm);
         if(numext::real(c0) >= RealScalar(0))
           beta = -beta;
         tval(Qidx(0)) = 1;
         for (Index itq = 1; itq < nzcolQ; ++itq)
           tval(Qidx(itq)) /= (c0 - beta);
         tau = numext::conj((beta-c0) / beta);
           
       }
     }
 
     // Insert values in R
     for (Index  i = nzcolR-1; i >= 0; i--)
     {
       Index curIdx = Ridx(i);
       if(curIdx < nonzeroCol) 
       {
         m_R.insertBackByOuterInnerUnordered(col, curIdx) = tval(curIdx);
         tval(curIdx) = Scalar(0.);
       }
     }
 
     if(nonzeroCol < diagSize && abs(beta) >= pivotThreshold)
     {
       m_R.insertBackByOuterInner(col, nonzeroCol) = beta;
       // The householder coefficient
       m_hcoeffs(nonzeroCol) = tau;
       // Record the householder reflections
       for (Index itq = 0; itq < nzcolQ; ++itq)
       {
         Index iQ = Qidx(itq);
         m_Q.insertBackByOuterInnerUnordered(nonzeroCol,iQ) = tval(iQ);
         tval(iQ) = Scalar(0.);
       }
       nonzeroCol++;
       if(nonzeroCol<diagSize)
         m_Q.startVec(nonzeroCol);
     }
     else
     {
       // Zero pivot found: move implicitly this column to the end
       for (Index j = nonzeroCol; j < n-1; j++) 
         std::swap(m_pivotperm.indices()(j), m_pivotperm.indices()[j+1]);
       
       // Recompute the column elimination tree
       internal::coletree(m_pmat, m_etree, m_firstRowElt, m_pivotperm.indices().data());
       m_isEtreeOk = false;
     }
   }
   
   m_hcoeffs.tail(diagSize-nonzeroCol).setZero();
   
   // Finalize the column pointers of the sparse matrices R and Q
   m_Q.finalize();
   m_Q.makeCompressed();
   m_R.finalize();
   m_R.makeCompressed();
   m_isQSorted = false;
 
   m_nonzeropivots = nonzeroCol;
   
   if(nonzeroCol<n)
   {
     // Permute the triangular factor to put the 'dead' columns to the end
     QRMatrixType tempR(m_R);
     m_R = tempR * m_pivotperm;
     
     // Update the column permutation
     m_outputPerm_c = m_outputPerm_c * m_pivotperm;
   }
   
   m_isInitialized = true; 
   m_factorizationIsok = true;
   m_info = Success;
 }
 
 template <typename SparseQRType, typename Derived>
 struct SparseQR_QProduct : ReturnByValue<SparseQR_QProduct<SparseQRType, Derived> >
 {
   typedef typename SparseQRType::QRMatrixType MatrixType;
   typedef typename SparseQRType::Scalar Scalar;
   // Get the references 
   SparseQR_QProduct(const SparseQRType& qr, const Derived& other, bool transpose) : 
   m_qr(qr),m_other(other),m_transpose(transpose) {}
   inline Index rows() const { return m_qr.matrixQ().rows(); }
   inline Index cols() const { return m_other.cols(); }
   
   // Assign to a vector
   template<typename DesType>
   void evalTo(DesType& res) const
   {
     Index m = m_qr.rows();
     Index n = m_qr.cols();
     Index diagSize = (std::min)(m,n);
     res = m_other;
     if (m_transpose)
     {
       eigen_assert(m_qr.m_Q.rows() == m_other.rows() && "Non conforming object sizes");
       //Compute res = Q' * other column by column
       for(Index j = 0; j < res.cols(); j++){
         for (Index k = 0; k < diagSize; k++)
         {
           Scalar tau = Scalar(0);
           tau = m_qr.m_Q.col(k).dot(res.col(j));
           if(tau==Scalar(0)) continue;
           tau = tau * m_qr.m_hcoeffs(k);
           res.col(j) -= tau * m_qr.m_Q.col(k);
         }
       }
     }
     else
     {
       eigen_assert(m_qr.matrixQ().cols() == m_other.rows() && "Non conforming object sizes");
 
       res.conservativeResize(rows(), cols());
 
       // Compute res = Q * other column by column
       for(Index j = 0; j < res.cols(); j++)
       {
         Index start_k = internal::is_identity<Derived>::value ? numext::mini(j,diagSize-1) : diagSize-1;
         for (Index k = start_k; k >=0; k--)
         {
           Scalar tau = Scalar(0);
           tau = m_qr.m_Q.col(k).dot(res.col(j));
           if(tau==Scalar(0)) continue;
           tau = tau * numext::conj(m_qr.m_hcoeffs(k));
           res.col(j) -= tau * m_qr.m_Q.col(k);
         }
       }
     }
   }
   
   const SparseQRType& m_qr;
   const Derived& m_other;
   bool m_transpose; // TODO this actually means adjoint
 };
 
 template<typename SparseQRType>
 struct SparseQRMatrixQReturnType : public EigenBase<SparseQRMatrixQReturnType<SparseQRType> >
 {  
   typedef typename SparseQRType::Scalar Scalar;
   typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
   enum {
     RowsAtCompileTime = Dynamic,
     ColsAtCompileTime = Dynamic
   };
   explicit SparseQRMatrixQReturnType(const SparseQRType& qr) : m_qr(qr) {}
   template<typename Derived>
   SparseQR_QProduct<SparseQRType, Derived> operator*(const MatrixBase<Derived>& other)
   {
     return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(),false);
   }
   // To use for operations with the adjoint of Q
   SparseQRMatrixQTransposeReturnType<SparseQRType> adjoint() const
   {
     return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
   }
   inline Index rows() const { return m_qr.rows(); }
   inline Index cols() const { return m_qr.rows(); }
   // To use for operations with the transpose of Q FIXME this is the same as adjoint at the moment
   SparseQRMatrixQTransposeReturnType<SparseQRType> transpose() const
   {
     return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
   }
   const SparseQRType& m_qr;
 };
 
 // TODO this actually represents the adjoint of Q
 template<typename SparseQRType>
 struct SparseQRMatrixQTransposeReturnType
 {
   explicit SparseQRMatrixQTransposeReturnType(const SparseQRType& qr) : m_qr(qr) {}
   template<typename Derived>
   SparseQR_QProduct<SparseQRType,Derived> operator*(const MatrixBase<Derived>& other)
   {
     return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(), true);
   }
   const SparseQRType& m_qr;
 };
 
 namespace internal {
   
 template<typename SparseQRType>
 struct evaluator_traits<SparseQRMatrixQReturnType<SparseQRType> >
 {
   typedef typename SparseQRType::MatrixType MatrixType;
   typedef typename storage_kind_to_evaluator_kind<typename MatrixType::StorageKind>::Kind Kind;
   typedef SparseShape Shape;
 };
 
 template< typename DstXprType, typename SparseQRType>
 struct Assignment<DstXprType, SparseQRMatrixQReturnType<SparseQRType>, internal::assign_op<typename DstXprType::Scalar,typename DstXprType::Scalar>, Sparse2Sparse>
 {
   typedef SparseQRMatrixQReturnType<SparseQRType> SrcXprType;
   typedef typename DstXprType::Scalar Scalar;
   typedef typename DstXprType::StorageIndex StorageIndex;
   static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,Scalar> &/*func*/)
   {
     typename DstXprType::PlainObject idMat(src.rows(), src.cols());
     idMat.setIdentity();
     // Sort the sparse householder reflectors if needed
     const_cast<SparseQRType *>(&src.m_qr)->_sort_matrix_Q();
     dst = SparseQR_QProduct<SparseQRType, DstXprType>(src.m_qr, idMat, false);
   }
 };
 
 template< typename DstXprType, typename SparseQRType>
 struct Assignment<DstXprType, SparseQRMatrixQReturnType<SparseQRType>, internal::assign_op<typename DstXprType::Scalar,typename DstXprType::Scalar>, Sparse2Dense>
 {
   typedef SparseQRMatrixQReturnType<SparseQRType> SrcXprType;
   typedef typename DstXprType::Scalar Scalar;
   typedef typename DstXprType::StorageIndex StorageIndex;
   static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,Scalar> &/*func*/)
   {
     dst = src.m_qr.matrixQ() * DstXprType::Identity(src.m_qr.rows(), src.m_qr.rows());
   }
 };
 
 } // end namespace internal
 
 } // end namespace RivetEigen
 
 #endif

This page was automatically generated by the 2.3.7 LXR engine. The LXR team