SuiteSparseQRSupport.h revision 7faaa9f3f0df9d23790277834d426c3d992ac3ba
17faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// This file is part of Eigen, a lightweight C++ template library 27faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// for linear algebra. 37faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// 47faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Copyright (C) 2012 Desire Nuentsa <desire.nuentsa_wakam@inria.fr> 57faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// 67faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// This Source Code Form is subject to the terms of the Mozilla 77faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Public License v. 2.0. If a copy of the MPL was not distributed 87faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 97faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#ifndef EIGEN_SUITESPARSEQRSUPPORT_H 117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#define EIGEN_SUITESPARSEQRSUPPORT_H 127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeznamespace Eigen { 147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename MatrixType> class SPQR; 167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename SPQRType> struct SPQRMatrixQReturnType; 177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename SPQRType> struct SPQRMatrixQTransposeReturnType; 187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template <typename SPQRType, typename Derived> struct SPQR_QProduct; 197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez namespace internal { 207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template <typename SPQRType> struct traits<SPQRMatrixQReturnType<SPQRType> > 217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename SPQRType::MatrixType ReturnType; 237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez }; 247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template <typename SPQRType> struct traits<SPQRMatrixQTransposeReturnType<SPQRType> > 257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename SPQRType::MatrixType ReturnType; 277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez }; 287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template <typename SPQRType, typename Derived> struct traits<SPQR_QProduct<SPQRType, Derived> > 297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename Derived::PlainObject ReturnType; 317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez }; 327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } // End namespace internal 337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez/** 357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \ingroup SPQRSupport_Module 367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \class SPQR 377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \brief Sparse QR factorization based on SuiteSparseQR library 387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * This class is used to perform a multithreaded and multifrontal rank-revealing QR decomposition 407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * of sparse matrices. The result is then used to solve linear leasts_square systems. 417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Clearly, a QR factorization is returned such that A*P = Q*R where : 427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * P is the column permutation. Use colsPermutation() to get it. 447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Q is the orthogonal matrix represented as Householder reflectors. 467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Use matrixQ() to get an expression and matrixQ().transpose() to get the transpose. 477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * You can then apply it to a vector. 487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * R is the sparse triangular factor. Use matrixQR() to get it as SparseMatrix. 507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * NOTE : The Index type of R is always UF_long. You can get it with SPQR::Index 517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \tparam _MatrixType The type of the sparse matrix A, must be a column-major SparseMatrix<> 537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * NOTE 547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename _MatrixType> 577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezclass SPQR 587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez public: 607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename _MatrixType::Scalar Scalar; 617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename _MatrixType::RealScalar RealScalar; 627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef UF_long Index ; 637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef SparseMatrix<Scalar, ColMajor, Index> MatrixType; 647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef PermutationMatrix<Dynamic, Dynamic> PermutationType; 657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez public: 667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SPQR() 677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez : m_isInitialized(false), 687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_ordering(SPQR_ORDERING_DEFAULT), 697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_allow_tol(SPQR_DEFAULT_TOL), 707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_tolerance (NumTraits<Scalar>::epsilon()) 717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez cholmod_l_start(&m_cc); 737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SPQR(const _MatrixType& matrix) 767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez : m_isInitialized(false), 777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_ordering(SPQR_ORDERING_DEFAULT), 787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_allow_tol(SPQR_DEFAULT_TOL), 797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_tolerance (NumTraits<Scalar>::epsilon()) 807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez cholmod_l_start(&m_cc); 827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez compute(matrix); 837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ~SPQR() 867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SPQR_free(); 887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez cholmod_l_finish(&m_cc); 897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez void SPQR_free() 917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez cholmod_l_free_sparse(&m_H, &m_cc); 937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez cholmod_l_free_sparse(&m_cR, &m_cc); 947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez cholmod_l_free_dense(&m_HTau, &m_cc); 957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez std::free(m_E); 967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez std::free(m_HPinv); 977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez void compute(const _MatrixType& matrix) 1007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(m_isInitialized) SPQR_free(); 1027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixType mat(matrix); 1047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez cholmod_sparse A; 1057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez A = viewAsCholmod(mat); 1067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index col = matrix.cols(); 1077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_rank = SuiteSparseQR<Scalar>(m_ordering, m_tolerance, col, &A, 1087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez &m_cR, &m_E, &m_H, &m_HPinv, &m_HTau, &m_cc); 1097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (!m_cR) 1117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_info = NumericalIssue; 1137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_isInitialized = false; 1147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return; 1157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_info = Success; 1177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_isInitialized = true; 1187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_isRUpToDate = false; 1197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** 1217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Get the number of rows of the input matrix and the Q matrix 1227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 1237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez inline Index rows() const {return m_H->nrow; } 1247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** 1267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Get the number of columns of the input matrix. 1277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 1287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez inline Index cols() const { return m_cR->ncol; } 1297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \returns the solution X of \f$ A X = B \f$ using the current decomposition of A. 1317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \sa compute() 1337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 1347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename Rhs> 1357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez inline const internal::solve_retval<SPQR, Rhs> solve(const MatrixBase<Rhs>& B) const 1367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()"); 1387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(this->rows()==B.rows() 1397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez && "SPQR::solve(): invalid number of rows of the right hand side matrix B"); 1407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return internal::solve_retval<SPQR, Rhs>(*this, B.derived()); 1417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename Rhs, typename Dest> 1447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const 1457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()"); 1477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(b.cols()==1 && "This method is for vectors only"); 1487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez //Compute Q^T * b 1507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typename Dest::PlainObject y; 1517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez y = matrixQ().transpose() * b; 1527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Solves with the triangular matrix R 1537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index rk = this->rank(); 1547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez y.topRows(rk) = this->matrixR().topLeftCorner(rk, rk).template triangularView<Upper>().solve(y.topRows(rk)); 1557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez y.bottomRows(cols()-rk).setZero(); 1567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Apply the column permutation 1577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez dest.topRows(cols()) = colsPermutation() * y.topRows(cols()); 1587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_info = Success; 1607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \returns the sparse triangular factor R. It is a sparse matrix 1637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 1647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const MatrixType matrixR() const 1657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()"); 1677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(!m_isRUpToDate) { 1687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_R = viewAsEigen<Scalar,ColMajor, typename MatrixType::Index>(*m_cR); 1697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_isRUpToDate = true; 1707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return m_R; 1727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /// Get an expression of the matrix Q 1747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SPQRMatrixQReturnType<SPQR> matrixQ() const 1757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return SPQRMatrixQReturnType<SPQR>(*this); 1777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /// Get the permutation that was applied to columns of A 1797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez PermutationType colsPermutation() const 1807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(m_isInitialized && "Decomposition is not initialized."); 1827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index n = m_cR->ncol; 1837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez PermutationType colsPerm(n); 1847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(Index j = 0; j <n; j++) colsPerm.indices()(j) = m_E[j]; 1857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return colsPerm; 1867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** 1897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Gets the rank of the matrix. 1907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * It should be equal to matrixQR().cols if the matrix is full-rank 1917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 1927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index rank() const 1937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(m_isInitialized && "Decomposition is not initialized."); 1957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return m_cc.SPQR_istat[4]; 1967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /// Set the fill-reducing ordering method to be used 1987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez void setSPQROrdering(int ord) { m_ordering = ord;} 1997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /// Set the tolerance tol to treat columns with 2-norm < =tol as zero 2007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez void setPivotThreshold(const RealScalar& tol) { m_tolerance = tol; } 2017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \returns a pointer to the SPQR workspace */ 2037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez cholmod_common *cholmodCommon() const { return &m_cc; } 2047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Reports whether previous computation was successful. 2077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \returns \c Success if computation was succesful, 2097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \c NumericalIssue if the sparse QR can not be computed 2107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 2117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ComputationInfo info() const 2127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(m_isInitialized && "Decomposition is not initialized."); 2147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return m_info; 2157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez protected: 2177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez bool m_isInitialized; 2187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez bool m_analysisIsOk; 2197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez bool m_factorizationIsOk; 2207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez mutable bool m_isRUpToDate; 2217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez mutable ComputationInfo m_info; 2227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int m_ordering; // Ordering method to use, see SPQR's manual 2237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int m_allow_tol; // Allow to use some tolerance during numerical factorization. 2247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar m_tolerance; // treat columns with 2-norm below this tolerance as zero 2257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez mutable cholmod_sparse *m_cR; // The sparse R factor in cholmod format 2267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez mutable MatrixType m_R; // The sparse matrix R in Eigen format 2277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez mutable Index *m_E; // The permutation applied to columns 2287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez mutable cholmod_sparse *m_H; //The householder vectors 2297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez mutable Index *m_HPinv; // The row permutation of H 2307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez mutable cholmod_dense *m_HTau; // The Householder coefficients 2317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez mutable Index m_rank; // The rank of the matrix 2327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez mutable cholmod_common m_cc; // Workspace and parameters 2337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename ,typename > friend struct SPQR_QProduct; 2347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 2357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate <typename SPQRType, typename Derived> 2377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct SPQR_QProduct : ReturnByValue<SPQR_QProduct<SPQRType,Derived> > 2387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 2397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename SPQRType::Scalar Scalar; 2407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename SPQRType::Index Index; 2417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez //Define the constructor to get reference to argument types 2427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SPQR_QProduct(const SPQRType& spqr, const Derived& other, bool transpose) : m_spqr(spqr),m_other(other),m_transpose(transpose) {} 2437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez inline Index rows() const { return m_transpose ? m_spqr.rows() : m_spqr.cols(); } 2457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez inline Index cols() const { return m_other.cols(); } 2467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Assign to a vector 2477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename ResType> 2487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez void evalTo(ResType& res) const 2497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez cholmod_dense y_cd; 2517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez cholmod_dense *x_cd; 2527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int method = m_transpose ? SPQR_QTX : SPQR_QX; 2537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez cholmod_common *cc = m_spqr.cholmodCommon(); 2547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez y_cd = viewAsCholmod(m_other.const_cast_derived()); 2557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x_cd = SuiteSparseQR_qmult<Scalar>(method, m_spqr.m_H, m_spqr.m_HTau, m_spqr.m_HPinv, &y_cd, cc); 2567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez res = Matrix<Scalar,ResType::RowsAtCompileTime,ResType::ColsAtCompileTime>::Map(reinterpret_cast<Scalar*>(x_cd->x), x_cd->nrow, x_cd->ncol); 2577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez cholmod_l_free_dense(&x_cd, cc); 2587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const SPQRType& m_spqr; 2607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const Derived& m_other; 2617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez bool m_transpose; 2627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 2647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename SPQRType> 2657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct SPQRMatrixQReturnType{ 2667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SPQRMatrixQReturnType(const SPQRType& spqr) : m_spqr(spqr) {} 2687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename Derived> 2697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SPQR_QProduct<SPQRType, Derived> operator*(const MatrixBase<Derived>& other) 2707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(),false); 2727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SPQRMatrixQTransposeReturnType<SPQRType> adjoint() const 2747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return SPQRMatrixQTransposeReturnType<SPQRType>(m_spqr); 2767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // To use for operations with the transpose of Q 2787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SPQRMatrixQTransposeReturnType<SPQRType> transpose() const 2797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return SPQRMatrixQTransposeReturnType<SPQRType>(m_spqr); 2817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const SPQRType& m_spqr; 2837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 2847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename SPQRType> 2867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct SPQRMatrixQTransposeReturnType{ 2877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SPQRMatrixQTransposeReturnType(const SPQRType& spqr) : m_spqr(spqr) {} 2887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename Derived> 2897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SPQR_QProduct<SPQRType,Derived> operator*(const MatrixBase<Derived>& other) 2907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(), true); 2927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const SPQRType& m_spqr; 2947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 2957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeznamespace internal { 2977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename _MatrixType, typename Rhs> 2997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct solve_retval<SPQR<_MatrixType>, Rhs> 3007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez : solve_retval_base<SPQR<_MatrixType>, Rhs> 3017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 3027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef SPQR<_MatrixType> Dec; 3037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) 3047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename Dest> void evalTo(Dest& dst) const 3067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 3077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez dec()._solve(rhs(),dst); 3087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 3107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} // end namespace internal 3127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}// End namespace Eigen 3147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#endif 315