diff options
Diffstat (limited to 'runtimes/nn/depend/external/eigen/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h')
| -rw-r--r-- | runtimes/nn/depend/external/eigen/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h | 313 |
1 files changed, 0 insertions, 313 deletions
diff --git a/runtimes/nn/depend/external/eigen/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h b/runtimes/nn/depend/external/eigen/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h deleted file mode 100644 index 953d57c9d..000000000 --- a/runtimes/nn/depend/external/eigen/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h +++ /dev/null @@ -1,313 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2012 Desire Nuentsa <desire.nuentsa_wakam@inria.fr> -// Copyright (C) 2014 Gael Guennebaud <gael.guennebaud@inria.fr> -// -// This Source Code Form is subject to the terms of the Mozilla -// Public License v. 2.0. If a copy of the MPL was not distributed -// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. - -#ifndef EIGEN_SUITESPARSEQRSUPPORT_H -#define EIGEN_SUITESPARSEQRSUPPORT_H - -namespace Eigen { - - template<typename MatrixType> class SPQR; - template<typename SPQRType> struct SPQRMatrixQReturnType; - template<typename SPQRType> struct SPQRMatrixQTransposeReturnType; - template <typename SPQRType, typename Derived> struct SPQR_QProduct; - namespace internal { - template <typename SPQRType> struct traits<SPQRMatrixQReturnType<SPQRType> > - { - typedef typename SPQRType::MatrixType ReturnType; - }; - template <typename SPQRType> struct traits<SPQRMatrixQTransposeReturnType<SPQRType> > - { - typedef typename SPQRType::MatrixType ReturnType; - }; - template <typename SPQRType, typename Derived> struct traits<SPQR_QProduct<SPQRType, Derived> > - { - 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 SuiteSparse_long. You can get it with SPQR::Index - * - * \tparam _MatrixType The type of the sparse matrix A, must be a column-major SparseMatrix<> - * - * \implsparsesolverconcept - * - * - */ -template<typename _MatrixType> -class SPQR : public SparseSolverBase<SPQR<_MatrixType> > -{ - protected: - typedef SparseSolverBase<SPQR<_MatrixType> > Base; - using Base::m_isInitialized; - public: - typedef typename _MatrixType::Scalar Scalar; - typedef typename _MatrixType::RealScalar RealScalar; - typedef SuiteSparse_long StorageIndex ; - typedef SparseMatrix<Scalar, ColMajor, StorageIndex> MatrixType; - typedef Map<PermutationMatrix<Dynamic, Dynamic, StorageIndex> > PermutationType; - enum { - ColsAtCompileTime = Dynamic, - MaxColsAtCompileTime = Dynamic - }; - public: - SPQR() - : m_ordering(SPQR_ORDERING_DEFAULT), m_allow_tol(SPQR_DEFAULT_TOL), m_tolerance (NumTraits<Scalar>::epsilon()), m_useDefaultThreshold(true) - { - cholmod_l_start(&m_cc); - } - - explicit SPQR(const _MatrixType& matrix) - : m_ordering(SPQR_ORDERING_DEFAULT), m_allow_tol(SPQR_DEFAULT_TOL), m_tolerance (NumTraits<Scalar>::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) - { - RealScalar max2Norm = 0.0; - for (int j = 0; j < mat.cols(); j++) max2Norm = numext::maxi(max2Norm, mat.col(j).norm()); - if(max2Norm==RealScalar(0)) - max2Norm = RealScalar(1); - pivotThreshold = 20 * (mat.rows() + mat.cols()) * max2Norm * NumTraits<RealScalar>::epsilon(); - } - cholmod_sparse A; - A = viewAsCholmod(mat); - m_rows = matrix.rows(); - Index col = matrix.cols(); - m_rank = SuiteSparseQR<Scalar>(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_rows; } - - /** - * Get the number of columns of the input matrix. - */ - inline Index cols() const { return m_cR->ncol; } - - template<typename Rhs, typename Dest> - void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &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<Upper>().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<Scalar,ColMajor, typename MatrixType::StorageIndex>(*m_cR); - m_isRUpToDate = true; - } - return m_R; - } - /// Get an expression of the matrix Q - SPQRMatrixQReturnType<SPQR> matrixQ() const - { - return SPQRMatrixQReturnType<SPQR>(*this); - } - /// Get the permutation that was applied to columns of A - PermutationType colsPermutation() const - { - eigen_assert(m_isInitialized && "Decomposition is not initialized."); - return PermutationType(m_E, m_cR->ncol); - } - /** - * Gets the rank of the matrix. - * It should be equal to matrixQR().cols if the matrix is full-rank - */ - Index rank() const - { - eigen_assert(m_isInitialized && "Decomposition is not initialized."); - return m_cc.SPQR_istat[4]; - } - /// Set the fill-reducing ordering method to be used - void setSPQROrdering(int ord) { m_ordering = ord;} - /// Set the tolerance tol to treat columns with 2-norm < =tol as zero - void setPivotThreshold(const RealScalar& tol) - { - m_useDefaultThreshold = false; - m_tolerance = tol; - } - - /** \returns a pointer to the SPQR workspace */ - cholmod_common *cholmodCommon() const { return &m_cc; } - - - /** \brief Reports whether previous computation was successful. - * - * \returns \c Success if computation was succesful, - * \c NumericalIssue if the sparse QR can not be computed - */ - ComputationInfo info() const - { - eigen_assert(m_isInitialized && "Decomposition is not initialized."); - return m_info; - } - protected: - bool m_analysisIsOk; - bool m_factorizationIsOk; - mutable bool m_isRUpToDate; - mutable ComputationInfo m_info; - int m_ordering; // Ordering method to use, see SPQR's manual - int m_allow_tol; // Allow to use some tolerance during numerical factorization. - RealScalar m_tolerance; // treat columns with 2-norm below this tolerance as zero - mutable cholmod_sparse *m_cR; // The sparse R factor in cholmod format - mutable MatrixType m_R; // The sparse matrix R in Eigen format - mutable StorageIndex *m_E; // The permutation applied to columns - mutable cholmod_sparse *m_H; //The householder vectors - mutable StorageIndex *m_HPinv; // The row permutation of H - mutable cholmod_dense *m_HTau; // The Householder coefficients - mutable Index m_rank; // The rank of the matrix - mutable cholmod_common m_cc; // Workspace and parameters - bool m_useDefaultThreshold; // Use default threshold - Index m_rows; - template<typename ,typename > friend struct SPQR_QProduct; -}; - -template <typename SPQRType, typename Derived> -struct SPQR_QProduct : ReturnByValue<SPQR_QProduct<SPQRType,Derived> > -{ - typedef typename SPQRType::Scalar Scalar; - typedef typename SPQRType::StorageIndex StorageIndex; - //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<typename ResType> - 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<Scalar>(method, m_spqr.m_H, m_spqr.m_HTau, m_spqr.m_HPinv, &y_cd, cc); - res = Matrix<Scalar,ResType::RowsAtCompileTime,ResType::ColsAtCompileTime>::Map(reinterpret_cast<Scalar*>(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<typename SPQRType> -struct SPQRMatrixQReturnType{ - - SPQRMatrixQReturnType(const SPQRType& spqr) : m_spqr(spqr) {} - template<typename Derived> - SPQR_QProduct<SPQRType, Derived> operator*(const MatrixBase<Derived>& other) - { - return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(),false); - } - SPQRMatrixQTransposeReturnType<SPQRType> adjoint() const - { - return SPQRMatrixQTransposeReturnType<SPQRType>(m_spqr); - } - // To use for operations with the transpose of Q - SPQRMatrixQTransposeReturnType<SPQRType> transpose() const - { - return SPQRMatrixQTransposeReturnType<SPQRType>(m_spqr); - } - const SPQRType& m_spqr; -}; - -template<typename SPQRType> -struct SPQRMatrixQTransposeReturnType{ - SPQRMatrixQTransposeReturnType(const SPQRType& spqr) : m_spqr(spqr) {} - template<typename Derived> - SPQR_QProduct<SPQRType,Derived> operator*(const MatrixBase<Derived>& other) - { - return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(), true); - } - const SPQRType& m_spqr; -}; - -}// End namespace Eigen -#endif |