// This file is part of Eigen, a lightweight C++ template library // for linear algebra. // // Copyright (C) 2012 Desire Nuentsa // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. #ifndef EIGEN_SUITESPARSEQRSUPPORT_H #define EIGEN_SUITESPARSEQRSUPPORT_H namespace Eigen { template class SPQR; template struct SPQRMatrixQReturnType; template struct SPQRMatrixQTransposeReturnType; template struct SPQR_QProduct; namespace internal { template struct traits > { typedef typename SPQRType::MatrixType ReturnType; }; template struct traits > { typedef typename SPQRType::MatrixType ReturnType; }; template struct traits > { typedef typename Derived::PlainObject ReturnType; }; } // End namespace internal /** * \ingroup SPQRSupport_Module * \class SPQR * \brief Sparse QR factorization based on SuiteSparseQR library * * This class is used to perform a multithreaded and multifrontal rank-revealing QR decomposition * of sparse matrices. The result is then used to solve linear leasts_square systems. * Clearly, a QR factorization is returned such that A*P = Q*R where : * * P is the column permutation. Use colsPermutation() to get it. * * Q is the orthogonal matrix represented as Householder reflectors. * Use matrixQ() to get an expression and matrixQ().transpose() to get the transpose. * You can then apply it to a vector. * * R is the sparse triangular factor. Use matrixQR() to get it as SparseMatrix. * NOTE : The Index type of R is always UF_long. You can get it with SPQR::Index * * \tparam _MatrixType The type of the sparse matrix A, must be a column-major SparseMatrix<> * NOTE * */ template class SPQR { public: typedef typename _MatrixType::Scalar Scalar; typedef typename _MatrixType::RealScalar RealScalar; typedef UF_long Index ; typedef SparseMatrix MatrixType; typedef PermutationMatrix PermutationType; public: SPQR() : m_isInitialized(false), m_ordering(SPQR_ORDERING_DEFAULT), m_allow_tol(SPQR_DEFAULT_TOL), m_tolerance (NumTraits::epsilon()), m_useDefaultThreshold(true) { cholmod_l_start(&m_cc); } SPQR(const _MatrixType& matrix) : m_isInitialized(false), m_ordering(SPQR_ORDERING_DEFAULT), m_allow_tol(SPQR_DEFAULT_TOL), m_tolerance (NumTraits::epsilon()), m_useDefaultThreshold(true) { cholmod_l_start(&m_cc); compute(matrix); } ~SPQR() { SPQR_free(); cholmod_l_finish(&m_cc); } void SPQR_free() { cholmod_l_free_sparse(&m_H, &m_cc); cholmod_l_free_sparse(&m_cR, &m_cc); cholmod_l_free_dense(&m_HTau, &m_cc); std::free(m_E); std::free(m_HPinv); } void compute(const _MatrixType& matrix) { if(m_isInitialized) SPQR_free(); MatrixType mat(matrix); /* Compute the default threshold as in MatLab, see: * Tim Davis, "Algorithm 915, SuiteSparseQR: Multifrontal Multithreaded Rank-Revealing * Sparse QR Factorization, ACM Trans. on Math. Soft. 38(1), 2011, Page 8:3 */ RealScalar pivotThreshold = m_tolerance; if(m_useDefaultThreshold) { using std::max; RealScalar max2Norm = 0.0; for (int j = 0; j < mat.cols(); j++) max2Norm = (max)(max2Norm, mat.col(j).norm()); if(max2Norm==RealScalar(0)) max2Norm = RealScalar(1); pivotThreshold = 20 * (mat.rows() + mat.cols()) * max2Norm * NumTraits::epsilon(); } cholmod_sparse A; A = viewAsCholmod(mat); Index col = matrix.cols(); m_rank = SuiteSparseQR(m_ordering, pivotThreshold, col, &A, &m_cR, &m_E, &m_H, &m_HPinv, &m_HTau, &m_cc); if (!m_cR) { m_info = NumericalIssue; m_isInitialized = false; return; } m_info = Success; m_isInitialized = true; m_isRUpToDate = false; } /** * Get the number of rows of the input matrix and the Q matrix */ inline Index rows() const {return m_cR->nrow; } /** * Get the number of columns of the input matrix. */ inline Index cols() const { return m_cR->ncol; } /** \returns the solution X of \f\$ A X = B \f\$ using the current decomposition of A. * * \sa compute() */ template inline const internal::solve_retval solve(const MatrixBase& B) const { eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()"); eigen_assert(this->rows()==B.rows() && "SPQR::solve(): invalid number of rows of the right hand side matrix B"); return internal::solve_retval(*this, B.derived()); } template void _solve(const MatrixBase &b, MatrixBase &dest) const { eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()"); eigen_assert(b.cols()==1 && "This method is for vectors only"); //Compute Q^T * b typename Dest::PlainObject y, y2; y = matrixQ().transpose() * b; // Solves with the triangular matrix R Index rk = this->rank(); y2 = y; y.resize((std::max)(cols(),Index(y.rows())),y.cols()); y.topRows(rk) = this->matrixR().topLeftCorner(rk, rk).template triangularView().solve(y2.topRows(rk)); // Apply the column permutation // colsPermutation() performs a copy of the permutation, // so let's apply it manually: for(Index i = 0; i < rk; ++i) dest.row(m_E[i]) = y.row(i); for(Index i = rk; i < cols(); ++i) dest.row(m_E[i]).setZero(); // y.bottomRows(y.rows()-rk).setZero(); // dest = colsPermutation() * y.topRows(cols()); m_info = Success; } /** \returns the sparse triangular factor R. It is a sparse matrix */ const MatrixType matrixR() const { eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()"); if(!m_isRUpToDate) { m_R = viewAsEigen(*m_cR); m_isRUpToDate = true; } return m_R; } /// Get an expression of the matrix Q SPQRMatrixQReturnType matrixQ() const { return SPQRMatrixQReturnType(*this); } /// Get the permutation that was applied to columns of A PermutationType colsPermutation() const { eigen_assert(m_isInitialized && "Decomposition is not initialized."); Index n = m_cR->ncol; PermutationType colsPerm(n); for(Index j = 0; j friend struct SPQR_QProduct; }; template struct SPQR_QProduct : ReturnByValue > { typedef typename SPQRType::Scalar Scalar; typedef typename SPQRType::Index Index; //Define the constructor to get reference to argument types SPQR_QProduct(const SPQRType& spqr, const Derived& other, bool transpose) : m_spqr(spqr),m_other(other),m_transpose(transpose) {} inline Index rows() const { return m_transpose ? m_spqr.rows() : m_spqr.cols(); } inline Index cols() const { return m_other.cols(); } // Assign to a vector template void evalTo(ResType& res) const { cholmod_dense y_cd; cholmod_dense *x_cd; int method = m_transpose ? SPQR_QTX : SPQR_QX; cholmod_common *cc = m_spqr.cholmodCommon(); y_cd = viewAsCholmod(m_other.const_cast_derived()); x_cd = SuiteSparseQR_qmult(method, m_spqr.m_H, m_spqr.m_HTau, m_spqr.m_HPinv, &y_cd, cc); res = Matrix::Map(reinterpret_cast(x_cd->x), x_cd->nrow, x_cd->ncol); cholmod_l_free_dense(&x_cd, cc); } const SPQRType& m_spqr; const Derived& m_other; bool m_transpose; }; template struct SPQRMatrixQReturnType{ SPQRMatrixQReturnType(const SPQRType& spqr) : m_spqr(spqr) {} template SPQR_QProduct operator*(const MatrixBase& other) { return SPQR_QProduct(m_spqr,other.derived(),false); } SPQRMatrixQTransposeReturnType adjoint() const { return SPQRMatrixQTransposeReturnType(m_spqr); } // To use for operations with the transpose of Q SPQRMatrixQTransposeReturnType transpose() const { return SPQRMatrixQTransposeReturnType(m_spqr); } const SPQRType& m_spqr; }; template struct SPQRMatrixQTransposeReturnType{ SPQRMatrixQTransposeReturnType(const SPQRType& spqr) : m_spqr(spqr) {} template SPQR_QProduct operator*(const MatrixBase& other) { return SPQR_QProduct(m_spqr,other.derived(), true); } const SPQRType& m_spqr; }; namespace internal { template struct solve_retval, Rhs> : solve_retval_base, Rhs> { typedef SPQR<_MatrixType> Dec; EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) template void evalTo(Dest& dst) const { dec()._solve(rhs(),dst); } }; } // end namespace internal }// End namespace Eigen #endif