summaryrefslogtreecommitdiff
path: root/runtimes/nn/depend/external/eigen/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h
diff options
context:
space:
mode:
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.h313
1 files changed, 313 insertions, 0 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
new file mode 100644
index 000000000..953d57c9d
--- /dev/null
+++ b/runtimes/nn/depend/external/eigen/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h
@@ -0,0 +1,313 @@
+// 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