15c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles)// This file is part of Eigen, a lightweight C++ template library 25c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles)// for linear algebra. 35c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles)// 45c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles)// Copyright (C) 2012-2013 Desire Nuentsa <desire.nuentsa_wakam@inria.fr> 55c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles)// Copyright (C) 2012-2014 Gael Guennebaud <gael.guennebaud@inria.fr> 65c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles)// 75c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles)// This Source Code Form is subject to the terms of the Mozilla 85c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles)// Public License v. 2.0. If a copy of the MPL was not distributed 95c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles)// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 105c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 115c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles)#ifndef EIGEN_SPARSE_QR_H 125c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles)#define EIGEN_SPARSE_QR_H 135c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 145c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles)namespace Eigen { 155c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 165c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles)template<typename MatrixType, typename OrderingType> class SparseQR; 175c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles)template<typename SparseQRType> struct SparseQRMatrixQReturnType; 185c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles)template<typename SparseQRType> struct SparseQRMatrixQTransposeReturnType; 195c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles)template<typename SparseQRType, typename Derived> struct SparseQR_QProduct; 205c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles)namespace internal { 215c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) template <typename SparseQRType> struct traits<SparseQRMatrixQReturnType<SparseQRType> > 225c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) { 235c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) typedef typename SparseQRType::MatrixType ReturnType; 245c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) typedef typename ReturnType::Index Index; 255c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) typedef typename ReturnType::StorageKind StorageKind; 265c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) }; 275c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) template <typename SparseQRType> struct traits<SparseQRMatrixQTransposeReturnType<SparseQRType> > 285c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) { 295c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) typedef typename SparseQRType::MatrixType ReturnType; 3053e740f4a82e17f3ae59772501622dc354e42336Torne (Richard Coles) }; 315c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) template <typename SparseQRType, typename Derived> struct traits<SparseQR_QProduct<SparseQRType, Derived> > 325c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) { 335c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) typedef typename Derived::PlainObject ReturnType; 3493ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) }; 3551b2906e11752df6c18351cf520e30522d3b53a1Torne (Richard Coles)} // End namespace internal 36df95704c49daea886ddad70775bda23618d6274dBen Murdoch 3753e740f4a82e17f3ae59772501622dc354e42336Torne (Richard Coles)/** 3853e740f4a82e17f3ae59772501622dc354e42336Torne (Richard Coles) * \ingroup SparseQR_Module 3953e740f4a82e17f3ae59772501622dc354e42336Torne (Richard Coles) * \class SparseQR 4053e740f4a82e17f3ae59772501622dc354e42336Torne (Richard Coles) * \brief Sparse left-looking rank-revealing QR factorization 4153e740f4a82e17f3ae59772501622dc354e42336Torne (Richard Coles) * 42e52495584422c5edb5b2944981473a2e208da323Torne (Richard Coles) * This class implements a left-looking rank-revealing QR decomposition 43e52495584422c5edb5b2944981473a2e208da323Torne (Richard Coles) * of sparse matrices. When a column has a norm less than a given tolerance 44e52495584422c5edb5b2944981473a2e208da323Torne (Richard Coles) * it is implicitly permuted to the end. The QR factorization thus obtained is 4553e740f4a82e17f3ae59772501622dc354e42336Torne (Richard Coles) * given by A*P = Q*R where R is upper triangular or trapezoidal. 4651b2906e11752df6c18351cf520e30522d3b53a1Torne (Richard Coles) * 471e202183a5dc46166763171984b285173f8585e5Torne (Richard Coles) * P is the column permutation which is the product of the fill-reducing and the 481e202183a5dc46166763171984b285173f8585e5Torne (Richard Coles) * rank-revealing permutations. Use colsPermutation() to get it. 491e202183a5dc46166763171984b285173f8585e5Torne (Richard Coles) * 501e202183a5dc46166763171984b285173f8585e5Torne (Richard Coles) * Q is the orthogonal matrix represented as products of Householder reflectors. 511e202183a5dc46166763171984b285173f8585e5Torne (Richard Coles) * Use matrixQ() to get an expression and matrixQ().transpose() to get the transpose. 5253e740f4a82e17f3ae59772501622dc354e42336Torne (Richard Coles) * You can then apply it to a vector. 5309380295ba73501a205346becac22c6978e4671dTorne (Richard Coles) * 54f79f16f17ddc4f842d7b7a38603e280e94be826aTorne (Richard Coles) * R is the sparse triangular or trapezoidal matrix. The later occurs when A is rank-deficient. 55d5428f32f5d1719f774f62e19147104ca245a3abTorne (Richard Coles) * matrixR().topLeftCorner(rank(), rank()) always returns a triangular factor of full rank. 5609380295ba73501a205346becac22c6978e4671dTorne (Richard Coles) * 5753e740f4a82e17f3ae59772501622dc354e42336Torne (Richard Coles) * \tparam _MatrixType The type of the sparse matrix A, must be a column-major SparseMatrix<> 5853e740f4a82e17f3ae59772501622dc354e42336Torne (Richard Coles) * \tparam _OrderingType The fill-reducing ordering method. See the \link OrderingMethods_Module 5953e740f4a82e17f3ae59772501622dc354e42336Torne (Richard Coles) * OrderingMethods \endlink module for the list of built-in and external ordering methods. 6053e740f4a82e17f3ae59772501622dc354e42336Torne (Richard Coles) * 6153e740f4a82e17f3ae59772501622dc354e42336Torne (Richard Coles) * \warning The input sparse matrix A must be in compressed mode (see SparseMatrix::makeCompressed()). 6206f816c7c76bc45a15e452ade8a34e8af077693eTorne (Richard Coles) * 6306f816c7c76bc45a15e452ade8a34e8af077693eTorne (Richard Coles) */ 6406f816c7c76bc45a15e452ade8a34e8af077693eTorne (Richard Coles)template<typename _MatrixType, typename _OrderingType> 6506f816c7c76bc45a15e452ade8a34e8af077693eTorne (Richard Coles)class SparseQR 6606f816c7c76bc45a15e452ade8a34e8af077693eTorne (Richard Coles){ 6753e740f4a82e17f3ae59772501622dc354e42336Torne (Richard Coles) public: 68f79f16f17ddc4f842d7b7a38603e280e94be826aTorne (Richard Coles) typedef _MatrixType MatrixType; 6951b2906e11752df6c18351cf520e30522d3b53a1Torne (Richard Coles) typedef _OrderingType OrderingType; 7051b2906e11752df6c18351cf520e30522d3b53a1Torne (Richard Coles) typedef typename MatrixType::Scalar Scalar; 7153e740f4a82e17f3ae59772501622dc354e42336Torne (Richard Coles) typedef typename MatrixType::RealScalar RealScalar; 7253e740f4a82e17f3ae59772501622dc354e42336Torne (Richard Coles) typedef typename MatrixType::Index Index; 7319cde67944066db31e633d9e386f2aa9bf9fadb3Torne (Richard Coles) typedef SparseMatrix<Scalar,ColMajor,Index> QRMatrixType; 741e202183a5dc46166763171984b285173f8585e5Torne (Richard Coles) typedef Matrix<Index, Dynamic, 1> IndexVector; 751e202183a5dc46166763171984b285173f8585e5Torne (Richard Coles) typedef Matrix<Scalar, Dynamic, 1> ScalarVector; 761e202183a5dc46166763171984b285173f8585e5Torne (Richard Coles) typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType; 77591b958dee2cf159d33a0b931e6231072eaf38d5Ben Murdoch public: 785c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) SparseQR () : m_isInitialized(false), m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false) 795c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) { } 805c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 815c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) /** Construct a QR factorization of the matrix \a mat. 825c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) * 835c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) * \warning The matrix \a mat must be in compressed mode (see SparseMatrix::makeCompressed()). 845c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) * 85323480423219ecd77329f8326dc5e0e3b50926d4Torne (Richard Coles) * \sa compute() 86f523d2789ac2f83c4eca0ee4d5161bfdb5f2d052Torne (Richard Coles) */ 875c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) SparseQR(const MatrixType& mat) : m_isInitialized(false), m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false) 88f523d2789ac2f83c4eca0ee4d5161bfdb5f2d052Torne (Richard Coles) { 89f523d2789ac2f83c4eca0ee4d5161bfdb5f2d052Torne (Richard Coles) compute(mat); 905c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) } 915c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 925c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) /** Computes the QR factorization of the sparse matrix \a mat. 935c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) * 945c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) * \warning The matrix \a mat must be in compressed mode (see SparseMatrix::makeCompressed()). 95f523d2789ac2f83c4eca0ee4d5161bfdb5f2d052Torne (Richard Coles) * 965c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) * \sa analyzePattern(), factorize() 975c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) */ 985c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) void compute(const MatrixType& mat) 995c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) { 1005c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) analyzePattern(mat); 1015c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) factorize(mat); 102f79f16f17ddc4f842d7b7a38603e280e94be826aTorne (Richard Coles) } 1035c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) void analyzePattern(const MatrixType& mat); 1045c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) void factorize(const MatrixType& mat); 1055c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 10651b2906e11752df6c18351cf520e30522d3b53a1Torne (Richard Coles) /** \returns the number of rows of the represented matrix. 10751b2906e11752df6c18351cf520e30522d3b53a1Torne (Richard Coles) */ 1085c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) inline Index rows() const { return m_pmat.rows(); } 1095c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 1105c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) /** \returns the number of columns of the represented matrix. 1115c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) */ 1125c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) inline Index cols() const { return m_pmat.cols();} 1135c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 1145c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) /** \returns a const reference to the \b sparse upper triangular matrix R of the QR factorization. 1155c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) */ 1165c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) const QRMatrixType& matrixR() const { return m_R; } 1175c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 1185c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) /** \returns the number of non linearly dependent columns as determined by the pivoting threshold. 1195c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) * 1205c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) * \sa setPivotThreshold() 1215c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) */ 12209380295ba73501a205346becac22c6978e4671dTorne (Richard Coles) Index rank() const 123a9984bf9ddc3cf73fdae3f29134a2bab379e7029Ben Murdoch { 124bfe3590b1806e3ff18f46ee3af5d4b83078f305aTorne (Richard Coles) eigen_assert(m_isInitialized && "The factorization should be called first, use compute()"); 125c0e19a689c8ac22cdc96b291a8d33a5d3b0b34a4Torne (Richard Coles) return m_nonzeropivots; 1265c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) } 127926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) 128926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) /** \returns an expression of the matrix Q as products of sparse Householder reflectors. 129926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) * The common usage of this function is to apply it to a dense matrix or vector 13053e740f4a82e17f3ae59772501622dc354e42336Torne (Richard Coles) * \code 1315c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) * VectorXd B1, B2; 1325c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) * // Initialize B1 133323480423219ecd77329f8326dc5e0e3b50926d4Torne (Richard Coles) * B2 = matrixQ() * B1; 1345c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) * \endcode 135323480423219ecd77329f8326dc5e0e3b50926d4Torne (Richard Coles) * 136926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) * To get a plain SparseMatrix representation of Q: 1375c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) * \code 1385c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) * SparseMatrix<double> Q; 1395c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) * Q = SparseQR<SparseMatrix<double> >(A).matrixQ(); 140323480423219ecd77329f8326dc5e0e3b50926d4Torne (Richard Coles) * \endcode 141323480423219ecd77329f8326dc5e0e3b50926d4Torne (Richard Coles) * Internally, this call simply performs a sparse product between the matrix Q 142323480423219ecd77329f8326dc5e0e3b50926d4Torne (Richard Coles) * and a sparse identity matrix. However, due to the fact that the sparse 143323480423219ecd77329f8326dc5e0e3b50926d4Torne (Richard Coles) * reflectors are stored unsorted, two transpositions are needed to sort 144f523d2789ac2f83c4eca0ee4d5161bfdb5f2d052Torne (Richard Coles) * them before performing the product. 145323480423219ecd77329f8326dc5e0e3b50926d4Torne (Richard Coles) */ 146323480423219ecd77329f8326dc5e0e3b50926d4Torne (Richard Coles) SparseQRMatrixQReturnType<SparseQR> matrixQ() const 147323480423219ecd77329f8326dc5e0e3b50926d4Torne (Richard Coles) { return SparseQRMatrixQReturnType<SparseQR>(*this); } 148926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) 1495c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) /** \returns a const reference to the column permutation P that was applied to A such that A*P = Q*R 150926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) * It is the combination of the fill-in reducing permutation and numerical column pivoting. 151926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) */ 152926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) const PermutationType& colsPermutation() const 153926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) { 1545c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) eigen_assert(m_isInitialized && "Decomposition is not initialized."); 15551b2906e11752df6c18351cf520e30522d3b53a1Torne (Richard Coles) return m_outputPerm_c; 156926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) } 15751b2906e11752df6c18351cf520e30522d3b53a1Torne (Richard Coles) 1585c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) /** \returns A string describing the type of error. 1595c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) * This method is provided to ease debugging, not to handle errors. 1606f543c786fc42989f552b4daa774ca5ff32fa697Ben Murdoch */ 1618abfc5808a4e34d6e03867af8bc440dee641886fTorne (Richard Coles) std::string lastErrorMessage() const { return m_lastError; } 16251b2906e11752df6c18351cf520e30522d3b53a1Torne (Richard Coles) 163bfe3590b1806e3ff18f46ee3af5d4b83078f305aTorne (Richard Coles) /** \internal */ 16451b2906e11752df6c18351cf520e30522d3b53a1Torne (Richard Coles) template<typename Rhs, typename Dest> 1658abfc5808a4e34d6e03867af8bc440dee641886fTorne (Richard Coles) bool _solve(const MatrixBase<Rhs> &B, MatrixBase<Dest> &dest) const 1668abfc5808a4e34d6e03867af8bc440dee641886fTorne (Richard Coles) { 1675c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) eigen_assert(m_isInitialized && "The factorization should be called first, use compute()"); 1685c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix"); 169f91f5fa1608c2cdd9af1842fb5dadbe78275be2aBo Liu 1705c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) Index rank = this->rank(); 1715c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 1725c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) // Compute Q^T * b; 1735c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) typename Dest::PlainObject y, b; 1745c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) y = this->matrixQ().transpose() * B; 1755c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) b = y; 176d5428f32f5d1719f774f62e19147104ca245a3abTorne (Richard Coles) 1775c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) // Solve with the triangular matrix R 17809380295ba73501a205346becac22c6978e4671dTorne (Richard Coles) y.resize((std::max)(cols(),Index(y.rows())),y.cols()); 179f91f5fa1608c2cdd9af1842fb5dadbe78275be2aBo Liu y.topRows(rank) = this->matrixR().topLeftCorner(rank, rank).template triangularView<Upper>().solve(b.topRows(rank)); 1805c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) y.bottomRows(y.rows()-rank).setZero(); 1815c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 1825c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) // Apply the column permutation 1835c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) if (m_perm_c.size()) dest = colsPermutation() * y.topRows(cols()); 1845c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) else dest = y.topRows(cols()); 1855c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 1865c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) m_info = Success; 1875c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) return true; 1885c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) } 1895c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 1905c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 1915c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) /** Sets the threshold that is used to determine linearly dependent columns during the factorization. 192e69819bd8e388ea4ad1636a19aa6b2eed4952191Ben Murdoch * 193e69819bd8e388ea4ad1636a19aa6b2eed4952191Ben Murdoch * In practice, if during the factorization the norm of the column that has to be eliminated is below 194f79f16f17ddc4f842d7b7a38603e280e94be826aTorne (Richard Coles) * this threshold, then the entire column is treated as zero, and it is moved at the end. 195e69819bd8e388ea4ad1636a19aa6b2eed4952191Ben Murdoch */ 196e69819bd8e388ea4ad1636a19aa6b2eed4952191Ben Murdoch void setPivotThreshold(const RealScalar& threshold) 1975c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) { 1985c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) m_useDefaultThreshold = false; 1995c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) m_threshold = threshold; 2005c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) } 2015c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 2025c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) /** \returns the solution X of \f$ A X = B \f$ using the current decomposition of A. 2035c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) * 2045c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) * \sa compute() 2055c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) */ 2065c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) template<typename Rhs> 2075c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) inline const internal::solve_retval<SparseQR, Rhs> solve(const MatrixBase<Rhs>& B) const 2085c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) { 2095c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) eigen_assert(m_isInitialized && "The factorization should be called first, use compute()"); 2105c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix"); 2115c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) return internal::solve_retval<SparseQR, Rhs>(*this, B.derived()); 2125c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) } 2135c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) template<typename Rhs> 2145c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) inline const internal::sparse_solve_retval<SparseQR, Rhs> solve(const SparseMatrixBase<Rhs>& B) const 2155c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) { 2165c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) eigen_assert(m_isInitialized && "The factorization should be called first, use compute()"); 2175c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix"); 2185c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) return internal::sparse_solve_retval<SparseQR, Rhs>(*this, B.derived()); 2195c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) } 2205c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 2215c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) /** \brief Reports whether previous computation was successful. 2225c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) * 2235c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) * \returns \c Success if computation was successful, 2245c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) * \c NumericalIssue if the QR factorization reports a numerical problem 2255c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) * \c InvalidInput if the input matrix is invalid 2265c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) * 2275c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) * \sa iparm() 2285c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) */ 2295c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) ComputationInfo info() const 2305c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) { 2315c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) eigen_assert(m_isInitialized && "Decomposition is not initialized."); 2325c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) return m_info; 2335c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) } 234926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) 235926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) protected: 236926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) inline void sort_matrix_Q() 237926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) { 238926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) if(this->m_isQSorted) return; 2395c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) // The matrix Q is sorted during the transposition 2405c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) SparseMatrix<Scalar, RowMajor, Index> mQrm(this->m_Q); 2415c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) this->m_Q = mQrm; 2425c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) this->m_isQSorted = true; 2435c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) } 2445c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 2455c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 2465c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) protected: 2475c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) bool m_isInitialized; 2485c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) bool m_analysisIsok; 2495c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) bool m_factorizationIsok; 2505c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) mutable ComputationInfo m_info; 2515c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) std::string m_lastError; 2525c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) QRMatrixType m_pmat; // Temporary matrix 2535c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) QRMatrixType m_R; // The triangular factor matrix 2545c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) QRMatrixType m_Q; // The orthogonal reflectors 2555c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) ScalarVector m_hcoeffs; // The Householder coefficients 256a854de003a23bf3c7f95ec0f8154ada64092ff5cTorne (Richard Coles) PermutationType m_perm_c; // Fill-reducing Column permutation 2575c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) PermutationType m_pivotperm; // The permutation for rank revealing 2585c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) PermutationType m_outputPerm_c; // The final column permutation 25906f816c7c76bc45a15e452ade8a34e8af077693eTorne (Richard Coles) RealScalar m_threshold; // Threshold to determine null Householder reflections 2605c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) bool m_useDefaultThreshold; // Use default threshold 2615c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) Index m_nonzeropivots; // Number of non zero pivots found 2625c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) IndexVector m_etree; // Column elimination tree 2635c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) IndexVector m_firstRowElt; // First element in each row 2645c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) bool m_isQSorted; // whether Q is sorted or not 2655c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 2665c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) template <typename, typename > friend struct SparseQR_QProduct; 2675c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) template <typename > friend struct SparseQRMatrixQReturnType; 2685c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 2695c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles)}; 2705c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 2715c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles)/** \brief Preprocessing step of a QR factorization 2725c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) * 2735c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) * \warning The matrix \a mat must be in compressed mode (see SparseMatrix::makeCompressed()). 2745c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) * 2755c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) * In this step, the fill-reducing permutation is computed and applied to the columns of A 2765c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) * and the column elimination tree is computed as well. Only the sparsity pattern of \a mat is exploited. 2775c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) * 2785c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) * \note In this step it is assumed that there is no empty row in the matrix \a mat. 2795c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) */ 2805c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles)template <typename MatrixType, typename OrderingType> 2815c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles)void SparseQR<MatrixType,OrderingType>::analyzePattern(const MatrixType& mat) 2825c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles){ 2835c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) eigen_assert(mat.isCompressed() && "SparseQR requires a sparse matrix in compressed mode. Call .makeCompressed() before passing it to SparseQR"); 2845c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) // Compute the column fill reducing ordering 2855c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) OrderingType ord; 2865c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) ord(mat, m_perm_c); 2875c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) Index n = mat.cols(); 2885c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) Index m = mat.rows(); 2895c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) Index diagSize = (std::min)(m,n); 2905c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 2915c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) if (!m_perm_c.size()) 2925c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) { 2935c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) m_perm_c.resize(n); 2945c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) m_perm_c.indices().setLinSpaced(n, 0,n-1); 2955c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) } 2965c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 2975c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) // Compute the column elimination tree of the permuted matrix 2985c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) m_outputPerm_c = m_perm_c.inverse(); 2995c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) internal::coletree(mat, m_etree, m_firstRowElt, m_outputPerm_c.indices().data()); 3005c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 3015c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) m_R.resize(m, n); 3025c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) m_Q.resize(m, diagSize); 3035c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 3045c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) // Allocate space for nonzero elements : rough estimation 3055c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) m_R.reserve(2*mat.nonZeros()); //FIXME Get a more accurate estimation through symbolic factorization with the etree 3065c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) m_Q.reserve(2*mat.nonZeros()); 3075c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) m_hcoeffs.resize(diagSize); 3085c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) m_analysisIsok = true; 3095c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles)} 3105c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 3115c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles)/** \brief Performs the numerical QR factorization of the input matrix 3125c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) * 31351b2906e11752df6c18351cf520e30522d3b53a1Torne (Richard Coles) * The function SparseQR::analyzePattern(const MatrixType&) must have been called beforehand with 3145c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) * a matrix having the same sparsity pattern than \a mat. 31551b2906e11752df6c18351cf520e30522d3b53a1Torne (Richard Coles) * 3165c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) * \param mat The sparse column-major matrix 3175c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) */ 31851b2906e11752df6c18351cf520e30522d3b53a1Torne (Richard Coles)template <typename MatrixType, typename OrderingType> 3195c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles)void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat) 32051b2906e11752df6c18351cf520e30522d3b53a1Torne (Richard Coles){ 3215c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) using std::abs; 3225c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) using std::max; 3235c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 3245c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) eigen_assert(m_analysisIsok && "analyzePattern() should be called before this step"); 325c0e19a689c8ac22cdc96b291a8d33a5d3b0b34a4Torne (Richard Coles) Index m = mat.rows(); 3265c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) Index n = mat.cols(); 3275c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) Index diagSize = (std::min)(m,n); 3285c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) IndexVector mark((std::max)(m,n)); mark.setConstant(-1); // Record the visited nodes 3295c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) IndexVector Ridx(n), Qidx(m); // Store temporarily the row indexes for the current column of R and Q 3305c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) Index nzcolR, nzcolQ; // Number of nonzero for the current column of R and Q 3315c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) ScalarVector tval(m); // The dense vector used to compute the current column 3325c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) RealScalar pivotThreshold = m_threshold; 3335c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 3345c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) m_pmat = mat; 335c0e19a689c8ac22cdc96b291a8d33a5d3b0b34a4Torne (Richard Coles) m_pmat.uncompress(); // To have the innerNonZeroPtr allocated 3365c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) // Apply the fill-in reducing permutation lazily: 3375c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) for (int i = 0; i < n; i++) 338f5e4ad553afbc08dd2e729bb77e937a9a94d5827Torne (Richard Coles) { 3395c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) Index p = m_perm_c.size() ? m_perm_c.indices()(i) : i; 340f5e4ad553afbc08dd2e729bb77e937a9a94d5827Torne (Richard Coles) m_pmat.outerIndexPtr()[p] = mat.outerIndexPtr()[i]; 3415c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) m_pmat.innerNonZeroPtr()[p] = mat.outerIndexPtr()[i+1] - mat.outerIndexPtr()[i]; 3425c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) } 343591b958dee2cf159d33a0b931e6231072eaf38d5Ben Murdoch 3445c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) /* Compute the default threshold as in MatLab, see: 345591b958dee2cf159d33a0b931e6231072eaf38d5Ben Murdoch * Tim Davis, "Algorithm 915, SuiteSparseQR: Multifrontal Multithreaded Rank-Revealing 3465c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) * Sparse QR Factorization, ACM Trans. on Math. Soft. 38(1), 2011, Page 8:3 3475c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) */ 3485c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) if(m_useDefaultThreshold) 3495c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) { 3505c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) RealScalar max2Norm = 0.0; 3516f543c786fc42989f552b4daa774ca5ff32fa697Ben Murdoch for (int j = 0; j < n; j++) max2Norm = (max)(max2Norm, m_pmat.col(j).norm()); 3525c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) pivotThreshold = 20 * (m + n) * max2Norm * NumTraits<RealScalar>::epsilon(); 3535c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) } 3545c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 3558abfc5808a4e34d6e03867af8bc440dee641886fTorne (Richard Coles) // Initialize the numerical permutation 3568abfc5808a4e34d6e03867af8bc440dee641886fTorne (Richard Coles) m_pivotperm.setIdentity(n); 3575c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 3585c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) Index nonzeroCol = 0; // Record the number of valid pivots 3595c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) m_Q.startVec(0); 3605c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 361f5e4ad553afbc08dd2e729bb77e937a9a94d5827Torne (Richard Coles) // Left looking rank-revealing QR factorization: compute a column of R and Q at a time 362f5e4ad553afbc08dd2e729bb77e937a9a94d5827Torne (Richard Coles) for (Index col = 0; col < n; ++col) 36309380295ba73501a205346becac22c6978e4671dTorne (Richard Coles) { 36409380295ba73501a205346becac22c6978e4671dTorne (Richard Coles) mark.setConstant(-1); 36509380295ba73501a205346becac22c6978e4671dTorne (Richard Coles) m_R.startVec(col); 36609380295ba73501a205346becac22c6978e4671dTorne (Richard Coles) mark(nonzeroCol) = col; 367f5e4ad553afbc08dd2e729bb77e937a9a94d5827Torne (Richard Coles) Qidx(0) = nonzeroCol; 368f5e4ad553afbc08dd2e729bb77e937a9a94d5827Torne (Richard Coles) nzcolR = 0; nzcolQ = 1; 369f5e4ad553afbc08dd2e729bb77e937a9a94d5827Torne (Richard Coles) bool found_diag = nonzeroCol>=m; 37009380295ba73501a205346becac22c6978e4671dTorne (Richard Coles) tval.setZero(); 371f5e4ad553afbc08dd2e729bb77e937a9a94d5827Torne (Richard Coles) 372f5e4ad553afbc08dd2e729bb77e937a9a94d5827Torne (Richard Coles) // Symbolic factorization: find the nonzero locations of the column k of the factors R and Q, i.e., 373926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) // all the nodes (with indexes lower than rank) reachable through the column elimination tree (etree) rooted at node k. 3745c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) // Note: if the diagonal entry does not exist, then its contribution must be explicitly added, 37509380295ba73501a205346becac22c6978e4671dTorne (Richard Coles) // thus the trick with found_diag that permits to do one more iteration on the diagonal element if this one has not been found. 37609380295ba73501a205346becac22c6978e4671dTorne (Richard Coles) for (typename MatrixType::InnerIterator itp(m_pmat, col); itp || !found_diag; ++itp) 37709380295ba73501a205346becac22c6978e4671dTorne (Richard Coles) { 37809380295ba73501a205346becac22c6978e4671dTorne (Richard Coles) Index curIdx = nonzeroCol; 379926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) if(itp) curIdx = itp.row(); 3805c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) if(curIdx == nonzeroCol) found_diag = true; 3815c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 382d5428f32f5d1719f774f62e19147104ca245a3abTorne (Richard Coles) // Get the nonzeros indexes of the current column of R 38309380295ba73501a205346becac22c6978e4671dTorne (Richard Coles) Index st = m_firstRowElt(curIdx); // The traversal of the etree starts here 38409380295ba73501a205346becac22c6978e4671dTorne (Richard Coles) if (st < 0 ) 3855c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) { 3865c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) m_lastError = "Empty row found during numerical factorization"; 3875c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) m_info = InvalidInput; 3885c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) return; 3895c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) } 3905c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 3915c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) // Traverse the etree 39209380295ba73501a205346becac22c6978e4671dTorne (Richard Coles) Index bi = nzcolR; 3935c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) for (; mark(st) != col; st = m_etree(st)) 39409380295ba73501a205346becac22c6978e4671dTorne (Richard Coles) { 395e1f1df5f01594c0e62e751e4b46e779b85c2faa5Torne (Richard Coles) Ridx(nzcolR) = st; // Add this row to the list, 3965c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) mark(st) = col; // and mark this row as visited 3975c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) nzcolR++; 3985c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) } 3995c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 400e1f1df5f01594c0e62e751e4b46e779b85c2faa5Torne (Richard Coles) // Reverse the list to get the topological ordering 401c0e19a689c8ac22cdc96b291a8d33a5d3b0b34a4Torne (Richard Coles) Index nt = nzcolR-bi; 4025c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) for(Index i = 0; i < nt/2; i++) std::swap(Ridx(bi+i), Ridx(nzcolR-i-1)); 4035c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 404a854de003a23bf3c7f95ec0f8154ada64092ff5cTorne (Richard Coles) // Copy the current (curIdx,pcol) value of the input matrix 4055c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) if(itp) tval(curIdx) = itp.value(); 40609380295ba73501a205346becac22c6978e4671dTorne (Richard Coles) else tval(curIdx) = Scalar(0); 4075c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 4085c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) // Compute the pattern of Q(:,k) 4095c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) if(curIdx > nonzeroCol && mark(curIdx) != col ) 4105c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) { 4118abfc5808a4e34d6e03867af8bc440dee641886fTorne (Richard Coles) Qidx(nzcolQ) = curIdx; // Add this row to the pattern of Q, 4125c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) mark(curIdx) = col; // and mark it as visited 413a9984bf9ddc3cf73fdae3f29134a2bab379e7029Ben Murdoch nzcolQ++; 4145c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) } 4155c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) } 416323480423219ecd77329f8326dc5e0e3b50926d4Torne (Richard Coles) 4175c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) // Browse all the indexes of R(:,col) in reverse order 4185c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) for (Index i = nzcolR-1; i >= 0; i--) 4195c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) { 4205c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) Index curIdx = Ridx(i); 4215c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 42251b2906e11752df6c18351cf520e30522d3b53a1Torne (Richard Coles) // Apply the curIdx-th householder vector to the current column (temporarily stored into tval) 423f79f16f17ddc4f842d7b7a38603e280e94be826aTorne (Richard Coles) Scalar tdot(0); 424926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) 4255c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) // First compute q' * tval 4268abfc5808a4e34d6e03867af8bc440dee641886fTorne (Richard Coles) tdot = m_Q.col(curIdx).dot(tval); 427bfe3590b1806e3ff18f46ee3af5d4b83078f305aTorne (Richard Coles) 4288abfc5808a4e34d6e03867af8bc440dee641886fTorne (Richard Coles) tdot *= m_hcoeffs(curIdx); 4298abfc5808a4e34d6e03867af8bc440dee641886fTorne (Richard Coles) 43051b2906e11752df6c18351cf520e30522d3b53a1Torne (Richard Coles) // Then update tval = tval - q * tau 4315c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) // FIXME: tval -= tdot * m_Q.col(curIdx) should amount to the same (need to check/add support for efficient "dense ?= sparse") 432c0e19a689c8ac22cdc96b291a8d33a5d3b0b34a4Torne (Richard Coles) for (typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq) 4335c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) tval(itq.row()) -= itq.value() * tdot; 434926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) 4351e202183a5dc46166763171984b285173f8585e5Torne (Richard Coles) // Detect fill-in for the current column of Q 436926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) if(m_etree(Ridx(i)) == nonzeroCol) 4371e202183a5dc46166763171984b285173f8585e5Torne (Richard Coles) { 438926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) for (typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq) 4395c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) { 4405c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) Index iQ = itq.row(); 4415c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) if (mark(iQ) != col) 4425c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) { 4435c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) Qidx(nzcolQ++) = iQ; // Add this row to the pattern of Q, 4445c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) mark(iQ) = col; // and mark it as visited 4455c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) } 446a854de003a23bf3c7f95ec0f8154ada64092ff5cTorne (Richard Coles) } 4475c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) } 4485c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) } // End update current column 4495c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 4505c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) Scalar tau; 4515c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) RealScalar beta = 0; 4525c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 4535c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) if(nonzeroCol < diagSize) 4545c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) { 455a9984bf9ddc3cf73fdae3f29134a2bab379e7029Ben Murdoch // Compute the Householder reflection that eliminate the current column 45651b2906e11752df6c18351cf520e30522d3b53a1Torne (Richard Coles) // FIXME this step should call the Householder module. 4575c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) Scalar c0 = nzcolQ ? tval(Qidx(0)) : Scalar(0); 4585c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 459926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) // First, the squared norm of Q((col+1):m, col) 460926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) RealScalar sqrNorm = 0.; 4615c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) for (Index itq = 1; itq < nzcolQ; ++itq) sqrNorm += numext::abs2(tval(Qidx(itq))); 462926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) if(sqrNorm == RealScalar(0) && numext::imag(c0) == RealScalar(0)) 4635c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) { 464926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) tau = RealScalar(0); 4655c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) beta = numext::real(c0); 4665c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) tval(Qidx(0)) = 1; 4675c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) } 468f79f16f17ddc4f842d7b7a38603e280e94be826aTorne (Richard Coles) else 469f79f16f17ddc4f842d7b7a38603e280e94be826aTorne (Richard Coles) { 4705c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) using std::sqrt; 4715c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) beta = sqrt(numext::abs2(c0) + sqrNorm); 4725c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) if(numext::real(c0) >= RealScalar(0)) 4735c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) beta = -beta; 4745c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) tval(Qidx(0)) = 1; 4755c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) for (Index itq = 1; itq < nzcolQ; ++itq) 4765c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) tval(Qidx(itq)) /= (c0 - beta); 4775c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) tau = numext::conj((beta-c0) / beta); 4785c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 4795c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) } 4805c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) } 481c0e19a689c8ac22cdc96b291a8d33a5d3b0b34a4Torne (Richard Coles) 4825c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) // Insert values in R 483926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) for (Index i = nzcolR-1; i >= 0; i--) 4845c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) { 4855c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) Index curIdx = Ridx(i); 4865c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) if(curIdx < nonzeroCol) 4875c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) { 4885c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) m_R.insertBackByOuterInnerUnordered(col, curIdx) = tval(curIdx); 4895c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) tval(curIdx) = Scalar(0.); 4905c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) } 4915c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) } 4925c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 4935c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) if(nonzeroCol < diagSize && abs(beta) >= pivotThreshold) 4945c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) { 4955c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) m_R.insertBackByOuterInner(col, nonzeroCol) = beta; 4965c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) // The householder coefficient 4975c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) m_hcoeffs(nonzeroCol) = tau; 4985c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) // Record the householder reflections 4995c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) for (Index itq = 0; itq < nzcolQ; ++itq) 5005c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) { 5015c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) Index iQ = Qidx(itq); 5025c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) m_Q.insertBackByOuterInnerUnordered(nonzeroCol,iQ) = tval(iQ); 5035c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) tval(iQ) = Scalar(0.); 5045c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) } 5055c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) nonzeroCol++; 5065c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) if(nonzeroCol<diagSize) 5075c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) m_Q.startVec(nonzeroCol); 5085c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) } 5095c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) else 5105c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) { 5115c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) // Zero pivot found: move implicitly this column to the end 5125c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) for (Index j = nonzeroCol; j < n-1; j++) 5135c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) std::swap(m_pivotperm.indices()(j), m_pivotperm.indices()[j+1]); 5145c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 5155c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) // Recompute the column elimination tree 51651b2906e11752df6c18351cf520e30522d3b53a1Torne (Richard Coles) internal::coletree(m_pmat, m_etree, m_firstRowElt, m_pivotperm.indices().data()); 51793ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) } 51809380295ba73501a205346becac22c6978e4671dTorne (Richard Coles) } 519a854de003a23bf3c7f95ec0f8154ada64092ff5cTorne (Richard Coles) 52093ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) m_hcoeffs.tail(diagSize-nonzeroCol).setZero(); 52193ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) 52293ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) // Finalize the column pointers of the sparse matrices R and Q 52393ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) m_Q.finalize(); 52493ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) m_Q.makeCompressed(); 52551b2906e11752df6c18351cf520e30522d3b53a1Torne (Richard Coles) m_R.finalize(); 52693ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) m_R.makeCompressed(); 52709380295ba73501a205346becac22c6978e4671dTorne (Richard Coles) m_isQSorted = false; 52809380295ba73501a205346becac22c6978e4671dTorne (Richard Coles) 52993ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) m_nonzeropivots = nonzeroCol; 53093ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) 53193ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) if(nonzeroCol<n) 53293ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) { 53393ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) // Permute the triangular factor to put the 'dead' columns to the end 53451b2906e11752df6c18351cf520e30522d3b53a1Torne (Richard Coles) MatrixType tempR(m_R); 53593ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) m_R = tempR * m_pivotperm; 53609380295ba73501a205346becac22c6978e4671dTorne (Richard Coles) 537a854de003a23bf3c7f95ec0f8154ada64092ff5cTorne (Richard Coles) // Update the column permutation 53893ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) m_outputPerm_c = m_outputPerm_c * m_pivotperm; 53993ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) } 54093ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) 54193ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) m_isInitialized = true; 54293ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) m_factorizationIsok = true; 54351b2906e11752df6c18351cf520e30522d3b53a1Torne (Richard Coles) m_info = Success; 54493ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles)} 54509380295ba73501a205346becac22c6978e4671dTorne (Richard Coles) 546a854de003a23bf3c7f95ec0f8154ada64092ff5cTorne (Richard Coles)namespace internal { 54793ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) 54893ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles)template<typename _MatrixType, typename OrderingType, typename Rhs> 54993ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles)struct solve_retval<SparseQR<_MatrixType,OrderingType>, Rhs> 55093ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) : solve_retval_base<SparseQR<_MatrixType,OrderingType>, Rhs> 55193ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles){ 55251b2906e11752df6c18351cf520e30522d3b53a1Torne (Richard Coles) typedef SparseQR<_MatrixType,OrderingType> Dec; 55393ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) 55409380295ba73501a205346becac22c6978e4671dTorne (Richard Coles) 555a854de003a23bf3c7f95ec0f8154ada64092ff5cTorne (Richard Coles) template<typename Dest> void evalTo(Dest& dst) const 55693ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) { 55793ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) dec()._solve(rhs(),dst); 55893ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) } 55993ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles)}; 56093ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles)template<typename _MatrixType, typename OrderingType, typename Rhs> 56151b2906e11752df6c18351cf520e30522d3b53a1Torne (Richard Coles)struct sparse_solve_retval<SparseQR<_MatrixType, OrderingType>, Rhs> 56293ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) : sparse_solve_retval_base<SparseQR<_MatrixType, OrderingType>, Rhs> 56309380295ba73501a205346becac22c6978e4671dTorne (Richard Coles){ 564a854de003a23bf3c7f95ec0f8154ada64092ff5cTorne (Richard Coles) typedef SparseQR<_MatrixType, OrderingType> Dec; 56593ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec, Rhs) 56693ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) 56793ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) template<typename Dest> void evalTo(Dest& dst) const 56893ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) { 56993ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) this->defaultEvalTo(dst); 57051b2906e11752df6c18351cf520e30522d3b53a1Torne (Richard Coles) } 57193ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles)}; 57209380295ba73501a205346becac22c6978e4671dTorne (Richard Coles)} // end namespace internal 573a854de003a23bf3c7f95ec0f8154ada64092ff5cTorne (Richard Coles) 57493ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles)template <typename SparseQRType, typename Derived> 57593ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles)struct SparseQR_QProduct : ReturnByValue<SparseQR_QProduct<SparseQRType, Derived> > 57693ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles){ 57793ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) typedef typename SparseQRType::QRMatrixType MatrixType; 57893ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) typedef typename SparseQRType::Scalar Scalar; 57951b2906e11752df6c18351cf520e30522d3b53a1Torne (Richard Coles) typedef typename SparseQRType::Index Index; 58093ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) // Get the references 58109380295ba73501a205346becac22c6978e4671dTorne (Richard Coles) SparseQR_QProduct(const SparseQRType& qr, const Derived& other, bool transpose) : 582a854de003a23bf3c7f95ec0f8154ada64092ff5cTorne (Richard Coles) m_qr(qr),m_other(other),m_transpose(transpose) {} 58393ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) inline Index rows() const { return m_transpose ? m_qr.rows() : m_qr.cols(); } 58493ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) inline Index cols() const { return m_other.cols(); } 58593ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) 58693ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) // Assign to a vector 58793ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles) template<typename DesType> 5885c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) void evalTo(DesType& res) const 5895c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) { 5905c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) Index m = m_qr.rows(); 5915c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) Index n = m_qr.cols(); 5925c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) Index diagSize = (std::min)(m,n); 5935c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) res = m_other; 5945c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) if (m_transpose) 5955c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) { 5965c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) eigen_assert(m_qr.m_Q.rows() == m_other.rows() && "Non conforming object sizes"); 5975c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) //Compute res = Q' * other column by column 5985c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) for(Index j = 0; j < res.cols(); j++){ 5995c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) for (Index k = 0; k < diagSize; k++) 600926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) { 6015c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) Scalar tau = Scalar(0); 602926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) tau = m_qr.m_Q.col(k).dot(res.col(j)); 603926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) if(tau==Scalar(0)) continue; 604926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) tau = tau * m_qr.m_hcoeffs(k); 605926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) res.col(j) -= tau * m_qr.m_Q.col(k); 606926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) } 607926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) } 608926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) } 6095c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) else 610926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) { 611926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) eigen_assert(m_qr.m_Q.rows() == m_other.rows() && "Non conforming object sizes"); 6125c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) // Compute res = Q * other column by column 613926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) for(Index j = 0; j < res.cols(); j++) 614926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) { 6155c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) for (Index k = diagSize-1; k >=0; k--) 616926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) { 617926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) Scalar tau = Scalar(0); 618926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) tau = m_qr.m_Q.col(k).dot(res.col(j)); 6195c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) if(tau==Scalar(0)) continue; 620926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) tau = tau * m_qr.m_hcoeffs(k); 6215c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) res.col(j) -= tau * m_qr.m_Q.col(k); 6225c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) } 623926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) } 6245c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) } 625926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) } 6265c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 627926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) const SparseQRType& m_qr; 6285c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) const Derived& m_other; 629926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) bool m_transpose; 630926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles)}; 63153e740f4a82e17f3ae59772501622dc354e42336Torne (Richard Coles) 6325c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles)template<typename SparseQRType> 63353e740f4a82e17f3ae59772501622dc354e42336Torne (Richard Coles)struct SparseQRMatrixQReturnType : public EigenBase<SparseQRMatrixQReturnType<SparseQRType> > 634926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles){ 6355c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) typedef typename SparseQRType::Index Index; 6365c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) typedef typename SparseQRType::Scalar Scalar; 6375c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; 6385c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) SparseQRMatrixQReturnType(const SparseQRType& qr) : m_qr(qr) {} 639926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) template<typename Derived> 6405c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) SparseQR_QProduct<SparseQRType, Derived> operator*(const MatrixBase<Derived>& other) 641926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) { 6425c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(),false); 6435c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) } 6445c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) SparseQRMatrixQTransposeReturnType<SparseQRType> adjoint() const 64509380295ba73501a205346becac22c6978e4671dTorne (Richard Coles) { 6465c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr); 647a9984bf9ddc3cf73fdae3f29134a2bab379e7029Ben Murdoch } 6485c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) inline Index rows() const { return m_qr.rows(); } 6495c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) inline Index cols() const { return (std::min)(m_qr.rows(),m_qr.cols()); } 650c0e19a689c8ac22cdc96b291a8d33a5d3b0b34a4Torne (Richard Coles) // To use for operations with the transpose of Q 651926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) SparseQRMatrixQTransposeReturnType<SparseQRType> transpose() const 6525c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) { 6535c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr); 6545c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) } 6555c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) template<typename Dest> void evalTo(MatrixBase<Dest>& dest) const 6565c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) { 657926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) dest.derived() = m_qr.matrixQ() * Dest::Identity(m_qr.rows(), m_qr.rows()); 6585c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) } 6595c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) template<typename Dest> void evalTo(SparseMatrixBase<Dest>& dest) const 660926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) { 661926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) Dest idMat(m_qr.rows(), m_qr.rows()); 662926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) idMat.setIdentity(); 6635c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) // Sort the sparse householder reflectors if needed 664926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) const_cast<SparseQRType *>(&m_qr)->sort_matrix_Q(); 665926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) dest.derived() = SparseQR_QProduct<SparseQRType, Dest>(m_qr, idMat, false); 6665c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) } 6675c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 668926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) const SparseQRType& m_qr; 669c0e19a689c8ac22cdc96b291a8d33a5d3b0b34a4Torne (Richard Coles)}; 670926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) 671c0e19a689c8ac22cdc96b291a8d33a5d3b0b34a4Torne (Richard Coles)template<typename SparseQRType> 672926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles)struct SparseQRMatrixQTransposeReturnType 6735c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles){ 674926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) SparseQRMatrixQTransposeReturnType(const SparseQRType& qr) : m_qr(qr) {} 6755c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) template<typename Derived> 6761e202183a5dc46166763171984b285173f8585e5Torne (Richard Coles) SparseQR_QProduct<SparseQRType,Derived> operator*(const MatrixBase<Derived>& other) 677926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) { 6785c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(), true); 679926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles) } 6805c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) const SparseQRType& m_qr; 6815c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles)}; 6825c87bf8b86a7c82ef50fb7a89697d8e02e2553beTorne (Richard Coles) 68393ac45cfc74041c8ae536ce58a9534d46db2024eTorne (Richard Coles)} // end namespace Eigen 68409380295ba73501a205346becac22c6978e4671dTorne (Richard Coles) 6851e202183a5dc46166763171984b285173f8585e5Torne (Richard Coles)#endif 686926b001d589ce2f10facb93dd4b87578ea35a855Torne (Richard Coles)